How to obtain transfer function coefficients from data plot?

Hi all,

I am trying to model the air fuel ratio response to different NOx inputs, and I would like to obtain my system of ODE's from this data.

NOx_i = readtable('NOxIN.csv');
AFR_o = readtable('afrOUT.csv');

AFR_time=AFR_o{:,1};
AFR_val=AFR_o{:,2};

save('AFRout.mat','AFR_time','AFR_val');

load('AFRout.mat','AFR_val');
[AFR_num,AFR_den] = tfdata(AFR_val);

I have extracted the data points and imported them to matlab using readtable, currently as two seperate tables containing time and NOx, and time and AFR. I am currently trying this code to get the transfer function coefficients, but I don't think I am giving the tfdata function the proper input type. It needs to recieve a dynamic system input, but I'm unsure how to get my two tables and data sets into a system that can be read by tfdata.

Any help or suggestions would be greatly appreciated!

Answer 1

[Mathieu NOE]

hello

you need to use tfest

sys = tfest(data,np) estimates a continuous-time transfer function sys using the time-domain or frequency-domain data data and containing np poles. The number of zeros in sys is max(np-1,0).

just looking at the curves , it seems the relationship between the two is a simple as a first order polynomial (because the two curves are very similar only mirrored)

so if it does not work with tfest , I would try to fit a polynomial

[Derek Johnson]

I will try this out and let you know. Thank you!

[Derek Johnson]

I believe my data sets must be continous in order to make them a system. Right now my two time sets have different sampling times, but the same start and end. Meaning the AFR has data points at 0,21,72,143,178 seconds and the NOx has points at 0,14,49,90,178 (these are just random numbers to try to communicate the issue), I dont think tfest is working because my system isn't properly defined as a result of these different time inputs. Is there any way I can get these two sets of data on the same continous time domain?

Answer 2

[Star Strider]

‘Is there any way I can get these two sets of data on the same continous time domain?’

Use the Signal Processing Toolbox resample function on both signals using either the same time vector, or the same sampling frequency (these are both options). After that, you should be able to use the System Identification Toolbox funcitons with the resampled signals.

[Derek Johnson]

I must not be understanding the resample function properly because it is saying something is wrong with my arguments, I have tried a few different things and I suspect the issue is with the input signal, but the signal is a matrix with a column of time values and a column of AFR values, why would this not work as an input to resample it at 1800 Hz?

NOx_o = readmatrix('NOxOUT.csv');
AFR_i = readmatrix('afrIN.csv');

% AFR_time=AFR_i{:,1};
% AFR_val=AFR_i{:,2};
% AFR_time_flip=(AFR_i{:,1}).';
% AFR_val_flip=(AFR_i{:,2}).';
% 
% NOx_time=NOx_o{:,1};
% NOx_val=NOx_o{:,2};
% NOx_time_flip=(NOx_o{:,1}).';
% NOx_val_flip=(NOx_o{:,2}).';
% 
% save('AFRout.mat','AFR_time','AFR_val');
% 
% load('AFRout.mat');
% %[AFR_num,AFR_den] = tfdata(AFR_val);

% figure;
% plot(AFR_time, AFR_val);
% title('AFR Input Data');
% xlabel('Time (sec)');
% ylabel('Air Fuel Ratio (AFR)');
% 
% figure;
% plot(NOx_time, NOx_val);
% title('NOx Output Data');
% xlabel('Time (sec)');
% ylabel('NOx (ppm)');

t = (0:1:1800);
fs = 1;
AFR = resample(AFR_i,1800,1);
figure;
plot(t,AFR);

This is what the command window puts back:

Check for missing argument or incorrect argument data type in call to function 'resample'.

Error in projtest360 (line 33)

AFR = resample(AFR_i,1800,1);

[Star Strider]

I’m not certain what you’re doing. Since you specified a time vector ‘t’, see if using:

• y = resample(x,tx) resamples the values, x, of a signal sampled at the instants specified in vector tx. The function interpolates x linearly onto a vector of uniformly spaced instants with the same endpoints and number of samples as tx. NaNs are treated as missing data and are ignored.

will do what you want. (There is no specific link to that resample section.) Noite that there are other options.

I usually specify a sampling frequency, and let the function calculate everything else, using:

• [y,ty] = resample(x,tx,___) returns in ty the instants that correspond to the resampled signal.

with ‘tx’ being the original vector of sampling times. Use the same ‘t’ vector to resample both vectors. That should work. See the function documentation for details. Experiment to understand what the function is doing.

[Derek Johnson]

I finally got the function to return a figure with no errors. But the result is far from the original sample.

The result shown is with 1800 Hz sampling frequency, I have tried changing the frequency, the interpolation method, and I can't get this graph to change at all. Any suggestions on how to come closer to my original air fuel ratio curve, which I'll add below.

NOx_o = readtable('NOxOUT.csv');
AFR_i = readtable('afrIN.csv');

AFR_time=AFR_i{:,1};
AFR_val=AFR_i{:,2};
AFR_time_flip=(AFR_i{:,1}).';
AFR_val_flip=(AFR_i{:,2}).';

NOx_time=NOx_o{:,1};
NOx_val=NOx_o{:,2};
NOx_time_flip=(NOx_o{:,1}).';
NOx_val_flip=(NOx_o{:,2}).';

AFR_i = AFR_i{:,:};
NOx_o = NOx_o{:,:};

figure;
plot(AFR_time, AFR_val);
title('AFR Input Data');
xlabel('Time (sec)');
ylabel('Air Fuel Ratio (AFR)');
% 
% figure;
% plot(NOx_time, NOx_val);
% title('NOx Output Data');
% xlabel('Time (sec)');
% ylabel('NOx (ppm)');

%tx = (0:1:1800);
%fs = 1;
[y,ty] = resample(AFR_val,AFR_time,1800);
figure;
plot(ty,y);

And here is my workspace.

[Star Strider]

The code appears to be correct.

The problem may be the sampling frequency of 1800. You have 75 samples going from 0 to 1800 seconds.

I would begin with:

Fs = 1/mean(diff(AFR_time));

then round that up to an acceptable value (so if that is about 75/1800=0.042 samples/second round it up to either 0.05 or 0.1), and use that instead of 1800.

[Derek Johnson]

Is there a way I can evaluate the accuracy of the tfest returned system? I am not sure how to plot the system against various inputs, step inputs, etc... Or if there is another way I could see the accuracy of the systems. I right now have two transfer functions, one with 1 pole and one with 2 poles, and I don't know which is best to represent my system.

Note: The signals ended up coming out good with resample and tfest seems to have worked as of now. Thanks for all the help so far!

[Star Strider]

To determine the number of poles and zeros (and their approximate locations), calculate the Fourier transform of the time-domain isgnal, then plot only the imaginary part of the signal as a function of frequency. The pole locations will be the peaks (or at least the positive extrema), and the zeros will be the zero-crossings. This will at least tell the number of poles and zeros. Note that it could also reveal if there are poles or zeros at 0 or at +Inf, so it may be necessary to count those as well.

This should work, although in certain situations, it can be a bit ambiguous and require some interpretation in the context of what you know about the system.

This is the most reliable way I know of to estimate the number and approximate locations of the poles and zeros with any precision.

[Derek Johnson]

I think I am confused I how to get the Fourier transform of the time-domain signal. The tf gives me the num and den in laplace, then I am trying to get the tf into the time domain, and then plug in t values till 1800 seconds, then get the fft of that response. What am I misunderstanding here? I believe the error is that I cant plug t values into the t_sig because it is 1x1 sym, but I am not sure what else to try.

Here is my relevant code, workspace, and what t_sig is. Thanks for any help!

%% transfer function estimation

sys = tfest(data,2);
%sys2 = tfest(data,1);
[num,den] = tfdata(sys);
syms s
num_poly = poly2sym(cell2mat(num),s);
den_poly = poly2sym(cell2mat(den),s);
G_poly = num_poly/den_poly;
t_sig = ilaplace(G_poly);
t = (0:1:1800);
time_response = t_sig(t);
fft = fft(time_response);

[Star Strider]

I suggested plotting the imaginary component of the Fourier transform in order to estimate with reasonable accuracy the number of poles and zeros if the system you are identifyting. It likely has no other utility than that, since the Bode plot of the identified system will simulate it for you and plot the estimated transfer function.

Use the compare function to see how well the identified system matches the original data.