MatLab log analysis

[ztan]

I've been using MatLab to crunch my logfiles for a while now. Please post script examples that you find useful on this thread if you use MatLab.

Benefits: Fast, I mean properly fast. I can combine multiple scripts for analysis all at once (just loaded 112MB in 20 log files over 24.5 seconds).

Drawbacks: Expensive, unless you work/study in an institution that you access the software through.

Basic information including .csv loading and data filtering and binning scripts shared by nil5 on the Romraider forum. The scripts that I use are based on these and you'll need these scripts in the same folder that you run other scripts from if you want to use examples that I post:

http://www.romraider.com/forum/viewtopic.php?t=7233

Tips:

1. Please note that you will need to make sure that the scripts match your logging parameter names which are case sensitive.

2. Use a revision number at the end of each logfile - it will make it easy to batch analyze a lot of logfiles together.

Starting off:

The load_csv.m script that I use (I think I found it somewhere on the Romraider forums, I am not smart enough to have written it). Put in the appropriate revision number for the files that you want to analyze before running it:

clear 
close all
filepath = '.\';
filenames = ls([filepath '\*rev_x.xx.csv']);
s = size(filenames);
tic
disp('Loading data...');
for f = 1:s(1)
    fprintf(1,'(%d) %s\n',f,filenames(f,:));
    data_tmp = data_load([ filepath filenames(f,:) ]); 
    
    if (f == 1)
        data = data_tmp;
    else
        data = data_cat(data,data_tmp);
    end
end
toc
data0 = data;

This loads up all the .csv files in a single matrix with all the log parameters in line.

I'll post examples later on in the thread.

[diodelphi]

please continue
I love matlab, i think will be so helpful for tuners

[ztan]

Simple example: closed loop trims for scaling MAF

Load data using load_csv.m script

Data needed for analysis:
time, STFT, LTFT, MAF_V

Data needed for filtering:
ECT(Coolant temp), CL_OL status

data = data0;

%generate total trim
data.tottrim=data.STFT+data.LTFT;

%generate dMAF/dt
data.dMAF_V=diff(data.MAF_V)./diff(data.time);
data.dMAF_V=vertcat(data.MAF_V(1),data.dMAF_V);
data.dMAF_V=abs(data.dMAF_V);
fn = fieldnames(data);
N = length(fn);

%make sure engine is up to temp
ect_idx=find(data.ECT>80);
data=data_filter(data,ect_idx);

%discard wide variations in MAF
dMAF_V_idx=find(data.dMAF_V<0.2);
data=data_filter(data,dMAF_V_idx);

%find closed loop operation
CL_OL_idx=find(data.CL_OL==2);
data=data_filter(data,CL_OL_idx);

%MAF_V scale bins
MAF_V_Scale = [0.898,0.938,0.977,1.016,1.055,1.094,1.133,1.172,1.211,1.250,1.289,1.328,1.367,1.406,1.445,1.484,1.523,1.563,1.602,1.641,1.680,1.719,1.758,1.797,1.836,1.875,1.914,1.953,1.992,2.031,2.070,2.109,2.148,2.188,2.227,2.266,2.305,2.344,2.383,2.422,2.461,2.500,2.578,2.773,2.969,3.203,3.438,3.711,3.906,4.063,4.297,4.492,4.727,5.000];

MAF_V_Scale = reshape(MAF_V_Scale,54,1);

%put trim data in MAF_V bins

MAF_bin = data_bin(data,data.MAF_V, MAF_V_Scale);

MAF_corr_mean = zeros(length(MAF_bin),1);
MAF_corr_mode = zeros(length(MAF_bin),1);
MAF_corr_median = zeros(length(MAF_bin),1);
MAF_corr_std = zeros(length(MAF_bin),1);

%define minimum count number
MIN_CNT = 10;

for m= 1:length(MAF_bin)
       tottrim = MAF_bin(m).data.tottrim;
            
           MAF_corr_mean(m) = mean(tottrim);
           MAF_corr_mode(m) = mode(tottrim);
           MAF_corr_median(m) = median(tottrim);
           MAF_corr_std(m) = std(tottrim);
           MAF_corr_cnt(m) = length(tottrim);
       if (MAF_corr_cnt(m) < MIN_CNT)
           MAF_corr_mean(m) = 0;
           MAF_corr_mode(m) = 0;
           MAF_corr_median(m) = 0;
           MAF_corr_std(m) = 0;
       end
end

%%Plot data

figure(1);
scatter(data.MAF_V,data.tottrim,2);
title('MAF Trims');
hold on;
plot(MAF_V_Scale,MAF_corr_mean,'-*r','MarkerSize',4);
plot(MAF_V_Scale,MAF_corr_median,'-*g','MarkerSize',4);
hold off;

In the above example, the mean trim data for each MAF_V bin ends up in a 54:1 element matrix MAF_corr_mean which can be used to rescale your MAF.

Example trim plot:

[steve99]

Have you seen VGI,s maf scaling utility ?

[ztan]

Script for high level log analysis - display IAM, FLKC, and FBKC in one place over however many logfiles you've loaded:

figure(1);

% engine load
subplot(3,1,1);
plot(data.IAM,'.-');
title('IAM');
subplot(3,1,2);
plot(data.FBKC,'.-');
title('FBKC');
subplot(3,1,3);
plot(data.FLKC,'.-');
title('FLKC');

Example below from 7 daily drive logs; >110k log entries.
Time to load files: 5.2 sec, plot processed in under a second:

[ztan]

To isolate a part of a log for more detailed analysis, you can use the data filter script that nil5 provided.

For example, to look at the FLKC data in the above example, you could set an index and then filter the data appropriately:

%set index to the frames that you want
%example set index to select samples 20000-30000

idx=[2e4:3e4];

%filter data
%full stack of logfiles in data0

data=data_filter(data0,idx);

[ztan]

Example code is for analyzing FLKC decrements.
FLKC as well as FLKC_Offset cell need to be logged with other parameters.

The script identifies areas where FLKC has dropped without the offset cell changing, parses the data stack for samples pre and post FLKC decrement and plots whatever data that you want related to the knock event.

Average FLKC at each decrement is returned in a matrix that is load and RPM binned appropriately (FLKC_mean in this example).

Scripts for analyzing FBKC and IAM are similar but not as complex.

% find discrete decreases in FLKC

data.flkc_dec=zeros(length(data.time),1);

for m=2:1:length(data.time)
%note difference if FLKC value is lower than last FLKC
    if data.FLKC(m)<data.FLKC(m-1) & data.FLKC(m)<0;
            data.flkc_dec(m)=data.FLKC(m);
%discard data if FLKC offset cell has changed
            if  data.FLKC_Offset(m)~=data.FLKC_Offset(m-1)
                data.flkc_dec(m)=0;
            else
            end
    else
    end
end

idx_flkc = find(data.flkc_dec ~= 0);

      
%% Get indices that occurred immediately before and after the knock was detected and plot

SAMPLES_PRE = 20;
SAMPLES_POST = 19;

if (isempty(idx_flkc))
    display('FLKC: None.');
else

idx_flkc_pre = (idx_flkc*ones(1,(SAMPLES_PRE+1)))-(ones(length(idx_flkc),1)*(0:SAMPLES_PRE));
idx_flkc_post = (idx_flkc*ones(1,(SAMPLES_POST)))+(ones(length(idx_flkc),1)*(1:SAMPLES_POST));

idx_flkc2 = unique(cat(2,idx_flkc_pre,idx_flkc_post));

% generate data sets with pre and post samples

data_flkc2 = data_filter(data,idx_flkc2);

% Generate some plots

figure(1);

% RPM
subplot(8,1,1);
plot(data_flkc2.RPM,'.-');
title('RPM');
grid;
% Load
subplot(8,1,2);
plot(data_flkc2.Load,'.-r');
title('Load');
grid;
% timing
subplot(8,1,3);
plot(data_flkc2.TotTiming,'.-');
hold all;
plot(data_flkc2.BaseTiming,'.-r');
title('Timing');
grid;
% afr 
subplot(8,1,4);
plot(data_flkc2.AFR,'.-');
hold all;
plot(data_flkc2.AFR_command,'.-r');
title('AFR');
grid;
% Throttle
subplot(8,1,5);
plot(data_flkc2.Throttle,'.-');
title('Throttle');
grid;
% Trims
subplot(8,1,6);
plot(data_flkc2.LTFT,'.-');
hold all;
plot(data_flkc2.STFT,'.-r');
title('Trims');
grid;
% cl_ol
subplot(8,1,7);
plot(data_flkc2.CL_OL,'.-');
title('CL/OL Status');
grid;
% Knock
subplot(8,1,8);
plot(data_flkc2.FBKC,'.-');
hold all;
plot(data_flkc2.FLKC,'.-r');
hold all;
plot(data_flkc2.IAM,'.-g');
title('Knock')
grid;

end

%% data binning
data_flkc=data_filter(data,idx_flkc);

load_bins = [0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4];
rpm_bins = [800,1000,1200,1600,2000,2200,2400,2600,2800,3200,3600,4000,4400,4600,4800,5200,5600,6000,6400,6800,7000,7200,7400];

% Bin data by Load
data_flkc_load_binned = data_bin(data_flkc, data_flkc.Load, load_bins);

% Bin the data by RPM
flkc_binned = repmat(struct(),1,length(data_flkc_load_binned));

for m=1:length(data_flkc_load_binned)
    data_tmp = data_flkc_load_binned(m).data;
    flkc_binned(m).data = data_bin(data_tmp,data_tmp.RPM,rpm_bins);
end

% Now we have all the data binned!!!
% Construct a matrix of the parameter to observe
FLKC_mean = zeros(length(rpm_bins),length(load_bins));

% Put parameters in matrix
for m= 1:length(load_bins)
   for n = 1:length(rpm_bins)
        flkc2 = flkc_binned(m).data(n).data.flkc_dec;

        if (isempty(flkc2))
            FLKC_mean(n,m) = 0;
        else
            FLKC_mean(n,m) = mean(flkc2);
        end
   end
end

%Average FLKC for each FLKC decrement incident returned in FLKC_mean

%% Histogram plot

figure(100)

surf(load_bins,rpm_bins,FLKC_mean);
xlabel('load');
ylabel('rpm');
zlabel('FLKC decrements');
title('FLKC histogram');

[Joe-G]

This makes me want to learn matlab some more

[Joe-G]

Originally Posted by ztan

I've been using MatLab to crunch my logfiles for a while now. Please post script examples that you find useful on this thread if you use MatLab.

Benefits: Fast, I mean properly fast. I can combine multiple scripts for analysis all at once (just loaded 112MB in 20 log files over 24.5 seconds).

Drawbacks: Expensive, unless you work/study in an institution that you access the software through.

Basic information including .csv loading and data filtering and binning scripts shared by nil5 on the Romraider forum. The scripts that I use are based on these and you'll need these scripts in the same folder that you run other scripts from if you want to use examples that I post:

http://www.romraider.com/forum/viewtopic.php?t=7233

Tips:

1. Please note that you will need to make sure that the scripts match your logging parameter names which are case sensitive.

2. Use a revision number at the end of each logfile - it will make it easy to batch analyze a lot of logfiles together.

Starting off:

The load_csv.m script that I use (I think I found it somewhere on the Romraider forums, I am not smart enough to have written it). Put in the appropriate revision number for the files that you want to analyze before running it:

clear 
close all
filepath = '.\';
filenames = ls([filepath '\*rev_x.xx.csv']);
s = size(filenames);
tic
disp('Loading data...');
for f = 1:s(1)
    fprintf(1,'(%d) %s\n',f,filenames(f,:));
    data_tmp = data_load([ filepath filenames(f,:) ]); 
    
    if (f == 1)
        data = data_tmp;
    else
        data = data_cat(data,data_tmp);
    end
end
toc
data0 = data;

This loads up all the .csv files in a single matrix with all the log parameters in line.

I'll post examples later on in the thread.

Do you have any log files that work with these scripts? I'd like to give Matlab a shot. This would be a great time for me to learn how to use it.

[ztan]

Originally Posted by Joe-G

Do you have any log files that work with these scripts? I'd like to give Matlab a shot. This would be a great time for me to learn how to use it.

PM me with your email address and I'll send you some logs.
My logfiles are quite big and if you string a few drives together, you might need to be crunching many, many MB all together.

[ztan]

This example looks at adjustment of DI fuelling using the GDI flow table. Other tables can be used including GDI pressure target and GDI pressure multiplier.

When looking at the ROM disassimbly, a lot is indirectly referenced and quite difficult to pick apart.
For the GDI flow table in the A01G ROM:
X-axis is Fuel Pressure, FFF8936C is the RAM address for data used in the GDI Flow table
Y-axis looks like a Direct Injector Pulse Width value, FFF887E8 is the RAM address for data used in the GDI Flow table (I log this as Direct_IPW_0).
The return from the GDI flow lookup gets stored in FFF887EC and logs very closely with FFF887F0 (I believe this is Direct IPW).

This script can be run in 2 states, I've only shown the first:
1. Combined PI : DI operation
2. Full DI

Similar to MAF scaling with trims, this script only looks at closed loop operation.

An error between the two can be calculated to see if PI activation skews anything, but I have not found that particularly useful.

%reset dataset
data = data0;

%generate total trim
data.tottrim=data.STFT+data.LTFT;

%generate dMAF/dt
data.dMAF_V=diff(data.MAF_V)./diff(data.time);
%add first dMAF/dt value - expand dMAF_V by one
data.dMAF_V=vertcat(data.MAF_V(1),data.dMAF_V);
%make absolute value
data.dMAF_V=abs(data.dMAF_V);

%engine up to temp
ect_idx=find(data.ECT>80);
data=data_filter(data,ect_idx);

%dMAF_V/dt<0.2
dMAF_V_idx=find(data.dMAF_V<0.2);
data=data_filter(data,dMAF_V_idx);

%CL operation
CL_OL_idx=find(data.CL_OL==2);
data=data_filter(data,CL_OL_idx);

%set minimum count for matrix
MIN_CNT = 30;
%set plot limit
ZLIM = [-10,10  ];

%% Data binning to GDI flow table

%set fuel pressure and direct_ipw bins
fuelpress_bins = [2.0,4.0,8.0,12.0,16.0,20.0];
direct_ipw_bins = [0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0,3.0,4.0,5.0];

data_comb_fuelpress_binned = data_bin(data, data.FuelPress, fuelpress_bins);

data_comb_binned = repmat(struct(),1,length(data_comb_fuelpress_binned));

for m = 1:length(data_comb_fuelpress_binned)
    data_comb_tmp = data_comb_fuelpress_binned(m).data;
    data_comb_binned(m).data = data_bin(data_comb_tmp,data_comb_tmp.Direct_IPW_0,direct_ipw_bins); 
end

% construct a matrix of the total trim
comb_trim_mean = zeros(length(direct_ipw_bins),length(fuelpress_bins));
comb_trim_median = zeros(length(direct_ipw_bins),length(fuelpress_bins));
comb_trim_cnt = zeros(length(direct_ipw_bins),length(fuelpress_bins));

%corr_matrix_std = corr_matrix_mean;

for m= 1:length(fuelpress_bins)
   for n = 1:length(direct_ipw_bins)
       comb_trim = data_comb_binned(m).data(n).data.tottrim;
       if (isempty(comb_trim))
           comb_trim_mean(n,m) = 0;
           comb_trim_median(n,m) = 0;
           comb_trim_cnt(n,m) = 0;

       else
           comb_trim_mean(n,m) = mean(comb_trim);
           comb_trim_median(n,m) = median (comb_trim);
           comb_trim_cnt(n,m) = length (comb_trim);
       end
       if  (comb_trim_cnt(n,m) < MIN_CNT)
           comb_trim_mean(n,m) = 0;
           comb_trim_median(n,m) = 0;
       end    
   end
end

%%plot data
figure(1);

surf(fuelpress_bins,direct_ipw_bins,comb_trim_mean);
xlabel('FuelPress');
ylabel('Direct_IPW');
zlabel('TotalTrim');
zlim(ZLIM);
title('GDI Flow Rate Combined Trims');

The fuelling trims for each cell in the GDI flow rate table where there are more than 30 samples (MIN_CNT) is returned in the matrix comb_trim_mean. Coding for DI pressure target and multiplier are similar when you find the parameters that the ECU uses to input into the table.

[ztan]

This is a script to calculate your volumetric efficiency. Can be used to tune AVCS if you have a lot of time on your hands to maximize the amount of air running through your engine at a given rpm/load point.

Based on the romraider derivation at http://www.romraider.com/forum/viewtopic.php?t=8054 but using metric units that I log with.

data = data0;

%calculate VE
%VE=100*((MAF g/s * R L.kPa/K.mol * IAT K)/(MAP kPa * MMA g/mol))/(Disp cycle/L * RPM rev/min / (60 s/min * 2 rev/cycle))
%R = 1.205912 (MAP in psi), R = 0.08314472 (MAP in bar), R= 8.314472 (MAP in kPa), R = 62.363669 (MAP in mmHg).
%MMA = Molar Mass of Air Constant =28.97 g/mol

data.VE=100.*((data.MAF).*8.314472.*(data.IAT+273.15))./((data.MAP).*28.97)./(1.998.*(data.RPM)./(60*2));

MIN_CNT = 30;
ZLIM = [0,120];

%%
load_bins = [0.15,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.85,0.90,0.95,1.00,1.10,1.20,1.30,1.40];
rpm_bins = [800,1000,1200,1600,2000,2200,2400,2800,3200,3600,3800,4000,4200,4400,4600,4800,5200,5600,6000,6400,6800,7000,7200,7400];
%%
data_VE_load_binned = data_bin(data, data.Load, load_bins);

data_VE_binned = repmat(struct(),1,length(data_VE_load_binned));

for m = 1:length(data_VE_load_binned)
    data_VE_tmp = data_VE_load_binned(m).data;
    data_VE_binned(m).data = data_bin(data_VE_tmp,data_VE_tmp.RPM,rpm_bins); 
end

%%
% now we have all the data binned!!!
% construct a matrix of the parameter to observe
VE_mean = zeros(length(rpm_bins),length(load_bins));
VE_median = zeros(length(rpm_bins),length(load_bins));
VE_cnt = zeros(length(rpm_bins),length(load_bins));

%%
%corr_matrix_std = corr_matrix_mean;

for m = 1:length(load_bins)
   for n = 1:length(rpm_bins)
       VE_tmp = data_VE_binned(m).data(n).data.VE;
       if (isempty(VE_tmp))
           VE_mean(n,m) = 0;
           VE_median(n,m) = 0;
           VE_cnt(n,m) = 0;

       else
           VE_mean(n,m) = mean(VE_tmp);
           VE_median(n,m) = median (VE_tmp);
           VE_cnt(n,m) = length (VE_tmp);
       end
       if  (VE_cnt(n,m) < MIN_CNT)
           VE_mean(n,m) = 0;
           VE_median(n,m) = 0;
       end    
   end
end

%%

figure(1);

surf(load_bins,rpm_bins,VE_mean);
xlabel('Load');
ylabel('RPM');
zlabel('VE');
zlim(ZLIM);
title('Volumetric Efficiency');

Return for VE is given in the matrix VE_mean

Attached is stock intake and exhaust with Shiv's AVCS timing, you can see the torque dip in the VE plot.

You can also see at low RPM when the throttle is cracked open, the MAF sensor reads air coming into the intake, but the engine speed hasn't responded by the time the line is logged leading to a falsely high VE.

[ztan]

There are a few assumptions made here.

See Kodename47's PI : DI thread for the assumptions I have made and formula derivations.

Port injector duty cycle:

%calculate Port IDC
%IDC = 100 % * Port_IPW msec / ((6000 msec/min * 2 rev/cycle)/ RPM rev/min)

data.Port_IDC=(data.Port_IPW).*(data.RPM)./1200;

Direct injector duty cycle:

%calculate DI IDC
%max DI_IPW = 5ms, max FuelPressure = 20MPa

data.DI_IDC=(data.Direct_IPW).*(data.FuelPress);

Data binning follows code from Volumetric Efficiency script.
Figures are Port and DI duty using stock Port injectors and stock PI : DI map (23% port duty, 65% DI duty at redline and 1.2-1.3 g/rev load):