Determine the Bode diagram or a transfer function with input output data without tfest

In a partially unknown system, I measured a square-wave signal as input and the associated system response. Now I would like to use fft() to determine the transfer function or the Bode diagram. I have attached the measurements. I have the following code (also from this forum) but it results in this very bad bode plot. where is my mistake?

input = x3bzeit((23000:39000),2); 
output = x3bzeit((23000:39000),3);
time = x3bzeit((23000:39000),1)*10^(-6);

Fs = 1/mean(diff(time));            % Sampling Frequency
Fn = Fs/2;                          % Nyquist Frequency
L = numel(time);
FTinp = fft(input)/L;
FTout = fft(output)/L;
TF = FTout ./ FTinp;                % Transfer Function
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv);                   % Index Vector
figure
subplot(2,1,1)
plot(Fv, 20*log10(abs(TF(Iv))))
set(gca, 'XScale', 'log')

title('Amplitude')
ylabel('dB')
subplot(2,1,2)
plot(Fv, angle(TF(Iv))*180/pi)
title('Phase')
ylabel('°')

Answer 1

[Mathieu NOE]

hello

this is a better code ; a better transfer function estimate is based on the auto power spectra and cross spectrum between input and output;

see the plot inludes also the coherence , which must be close to one for the valid portion of the TF ; when coherence drops , this is mostly here because there is no energy in the signals in a frequency range (no excitation , only uncorrelated noise)

we can see here that the ringing of the time signal matches with a resonant second order low pass transfer funtion ; the resonance is around 50 kHz

I took the full signal length to make more averages , the result is better than just select a fraction of it

clc
clearvars

load('x3bzeit.mat');

% input = x3bzeit((23000:39000),2); 
% output = x3bzeit((23000:39000),3);
% time = x3bzeit((23000:39000),1)*10^(-6);

input = x3bzeit(:,2); 
output = x3bzeit(:,3);
time = x3bzeit(:,1)*10^(-6);


Fs = 1/mean(diff(time));            % Sampling Frequency
nfft = 10000;
window = hanning(nfft);
noverlap = round(0.75*nfft);
[Txy,Cxy,Fv] = tfe_and_coh(input,output,nfft,Fs,window,noverlap);
figure
subplot(3,1,1)
semilogx(Fv, 20*log10(abs(Txy)))
title('Amplitude')
ylabel('dB')
subplot(3,1,2)
semilogx(Fv, angle(Txy)*180/pi)
title('Phase')
ylabel('°')
subplot(3,1,3)
semilogx(Fv, Cxy)
title('Coherence')



%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = tfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function  and Coherence Estimate 
    
% compute PSD and CSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x , y if length is less than the window length
    x(nwind)=0;
    y(nwind)=0;  
    n=nwind;
end
x = x(:);		% Make sure x is a column vector
y = y(:);		% Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;

Pxx = zeros(nfft,1); 
Pyy = zeros(nfft,1); 
Pxy = zeros(nfft,1); 

for i=1:k
    xw = window.*x(index);
    yw = window.*y(index);
    index = index + (nwind - noverlap);
    
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft);
    Xx2 = abs(Xx).^2;
    Yy2 = abs(Yy).^2;
    Xy2 = Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pyy = Pyy + Yy2;
    Pxy = Pxy + Xy2;
    
end

% Select first half
if ~any(any(imag([x y])~=0))   % if x and y are not complex
    if rem(nfft,2)    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end

    Pxx = Pxx(select);
    Pyy = Pyy(select);
    Pxy = Pxy(select);
    
else
    select = 1:nfft;
end
Txy = Pxy ./ Pxx;                   % transfer function estimate 
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy);    % coherence function estimate 

f = (select - 1)'*Fs/nfft;

end

[Ricky Schulze]

Perfect, it works mostly for the other data too thanks a lot. I just have a question about the process what exactly do you charge for

Xy2 = Yy.*conj(Xx);

Is that the cross spectrum?

[Mathieu NOE]

Glad it works

yes Xy2 is the cross spectrum for one buffer of data

then Pxy is the sum for all buffers

[Ricky Schulze]

And why do you use Xy2 to determine the transfer function instead of Yy2?

I only know the process with Yy2

[Mathieu NOE]

there are different estimators - see the details here

What is a Frequency Response Function (FRF)? (siemens.com)

frequency response - Understanding the H1 and H2 estimators - Signal Processing Stack Exchange

On the use of the Hs estimator for the experimental assessment of transmissibility matrices (archives-ouvertes.fr)

[Mathieu NOE]

hello Ricky

hope you're doing well

do you mind accepting my answer ? tx

Answer 2

[Kumaresh Kumaresh]

Hello Mr. Mathieu,

Thank you for your code. I have used your code to estimate the transfer function based on power spectra for my input and output data.

However, coherence is maintained one, is that mean energy is always higher in the signal ?

Can you please share some insights ? Thank you

Here is the same code:

clc; clear all;

clc

clearvars

%importdata('velocityW.csv');

load('velocityW.mat');

input = ans(:,2);

output = ans(:,3);

time = ans(:,1); %*10^(-6);

Fs = 1/mean(diff(time)); % Sampling Frequency

nfft = 10000;

window = hanning(nfft);

noverlap = round(0.75*nfft);

[Txy,Cxy,Fv] = tfe_and_coh(input,output,nfft,Fs,window,noverlap);

figure

subplot(3,1,1)

semilogx(Fv, 20*log10(abs(Txy)))

title('Amplitude')

ylabel('dB')

subplot(3,1,2)

semilogx(Fv, angle(Txy)*180/pi)

title('Phase')

ylabel('°')

subplot(3,1,3)

semilogx(Fv, Cxy)

title('Coherence')

%%%%%%%%%%%%%%%%%%%%%%%

function [Txy,Cxy,f] = tfe_and_coh(x,y,nfft,Fs,window,noverlap)

% Transfer Function and Coherence Estimate

% compute PSD and CSD

window = window(:);

n = length(x); % Number of data points

nwind = length(window); % length of window

if n < nwind % zero-pad x , y if length is less than the window length

x(nwind)=0;

y(nwind)=0;

n=nwind;

end

x = x(:); % Make sure x is a column vector

y = y(:); % Make sure y is a column vector

k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows

% (k = fix(n/nwind) for noverlap=0)

index = 1:nwind;

Pxx = zeros(nfft,1);

Pyy = zeros(nfft,1);

Pxy = zeros(nfft,1);

for i=1:k

xw = window.*x(index);

yw = window.*y(index);

index = index + (nwind - noverlap);

Xx = fft(xw,nfft);

Yy = fft(yw,nfft);

Xx2 = abs(Xx).^2;

Yy2 = abs(Yy).^2;

Xy2 = Yy.*conj(Xx);

Pxx = Pxx + Xx2;

Pyy = Pyy + Yy2;

Pxy = Pxy + Xy2;

end

% Select first half

if ~any(any(imag([x y])~=0)) % if x and y are not complex

if rem(nfft,2) % nfft odd

select = [1:(nfft+1)/2];

else

select = [1:nfft/2+1]; % include DC AND Nyquist

end

Pxx = Pxx(select);

Pyy = Pyy(select);

Pxy = Pxy(select);

else

select = 1:nfft;

end

Txy = Pxy ./ Pxx; % transfer function estimate

Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate

f = (select - 1)'*Fs/nfft;

end

[Mathieu NOE]

hello @Kumaresh Kumaresh

coherence express the ratio of power in the output signal vs the input signal . It's a number between 0 and 1 .

Coherence (signal processing) - Wikipedia

if you make only one average (computation) because you choose nfft as long or even longer than the signals length, then mathematically coherece can only be one.

to see the "true" coherence you must have more than one average in the TFE computation.

here for example with nfft = 512 (also take a power of two for nfft to have faster fft computation)

Notice than now your coherenece is no more 1 all accross the frequency range. I would even say that'is is very low , so not very good identification IMO - how were those signas generated and acquired ?

load('velocityW.mat');
input = ans(:,2); 
output = ans(:,3);
time = ans(:,1);    %*10^(-6);
Fs = 1/mean(diff(time));            % Sampling Frequency

nfft = 512;
window = hanning(nfft);
noverlap = round(0.75*nfft);
[Txy,Cxy,Fv] = tfe_and_coh(input,output,nfft,Fs,window,noverlap);
figure
subplot(3,1,1)
semilogx(Fv, 20*log10(abs(Txy)))
title('Amplitude')
ylabel('dB')
subplot(3,1,2)
semilogx(Fv, angle(Txy)*180/pi)
title('Phase')
ylabel('°')
subplot(3,1,3)
semilogx(Fv, Cxy)
title('Coherence')
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = tfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function  and Coherence Estimate 
    
% compute PSD and CSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x , y if length is less than the window length
    x(nwind)=0;
    y(nwind)=0;  
    n=nwind;
end
x = x(:);		% Make sure x is a column vector
y = y(:);		% Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1); 
Pyy = zeros(nfft,1); 
Pxy = zeros(nfft,1); 
for i=1:k
    xw = window.*x(index);
    yw = window.*y(index);
    index = index + (nwind - noverlap);
    
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft);
    Xx2 = abs(Xx).^2;
    Yy2 = abs(Yy).^2;
    Xy2 = Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pyy = Pyy + Yy2;
    Pxy = Pxy + Xy2;
    
end
% Select first half
if ~any(any(imag([x y])~=0))   % if x and y are not complex
    if rem(nfft,2)    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end
    Pxx = Pxx(select);
    Pyy = Pyy(select);
    Pxy = Pxy(select);
    
else
    select = 1:nfft;
end
Txy = Pxy ./ Pxx;                   % transfer function estimate 
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy);    % coherence function estimate 
f = (select - 1)'*Fs/nfft;
end

[Mathieu NOE]

hello again

I was curious and so I plotted the time signals , input above and output below. the sinputs seems to be a noisy sinus with a sudden change in frequency , the output looks to me very strange , mostly noise (the input signal shape can be barely recognized).

this is not a good input signal to make a good quality FRF estimate on a large frequency range. you need to use a clean signal which energy is pread on a wide frequency range (random, bursts, sweeps)

[Kumaresh Kumaresh]

Hello Mr. Mathieu,

Thank you for your response.

Could you please elaborate me more about nfft. How nfft = 512 ? What is the optimum fit value ?

Answer for second comment:

During my simulation, first half was calculated very slowly (good sinusoidal behavior in flow pattern), whereas second half was calculated with faster time step. The reason why you observe disorderness and strange behavior in the second half profile. So, if I could able to extract the time data in a uniform manner (clear signal), I believe the data will be reasonable to interpret.

Thank you once again. I will update my results here in some days (after extracting the clean signal) to discuss with you once again.

[Mathieu NOE]

hello again

to the 1st comment :

• nfft must be less or equal to your signal length (I will not dig here into the topic of zero padding signals for fft enhanced frequency resolution - this is another debate)
• nfft is choosen depending of your data (quality) and frequency resolution :
• if the signals are very clean, you don't need many averages so nfft can be set to a larger value - and then your FRF estimation will have a refined frequency resolution (frequency increment : df = Fs/nfft)
• if your siganls are noisy, you may want to sacrifice a bit the frequency resolution and do more averages , so take a lower nfft value , and as a result the tfe_and_coh function will do more averages
• increasing the overlap factr allows you to do more averages for the same signal length (here we picked 75 %)

to the 2nd comment :

• if you use a signal that basically contains only 2 frequencies , you can only expect to perform a valid FRF estimation at those 2 frequencies , the rest of the FRF will be garbage (with coherence very low)

[Kumaresh Kumaresh]

Thank you for your explanation once again

[Mathieu NOE]

my pleasure !

[Kumaresh Kumaresh]

Hello Mr. Mathieu,

Hope you are doing good. Based on our previous conversations, I would like to share with you my plot with some questions.

here is the Bode plot for pressure variable

.

here is the Bode plot for velocity

Pressure is not the dominating parameter is my study, as we could see at higher frequencies, amplitude and phase becomes zero, whereas the velocity dominates at higher frequency, for input frequencies of 100 and 200Hz, when the amplitude becomes zero, the phase shift is higher oscillating between +180 to -180, proving unstability.

Kindly correct me if I am wrong and tell me if I missed something in the plot.

In addition, I would like to plot Nyquist or Nichols plot. Because in Nyquist or Nichols plots, the information will be in single chart + easier to determine stability, however frequency is not explicitly called out.

And I read that Nyquist plot easily handles time delay, which I am curious to plot and reading more about it.

Sorry for too many questions.

Share your suggestions please.

Thank you

[Mathieu NOE]

hello again and welcome back

I still have some issue interpreting your test protocol because - as far as I have understood so far , your test signal is more or less a sinusoid (f1 frequency) so normally most of the usefull signal should be at that frequency ... unless you have some noise (correlated or not) on top of your signal (can be input or output noises)

what I still don't understand is how you can draw conclusions in frequency ranges that are far away from f1 , it basically then only a transfer function estimate between low amplitude noisy signals so I don't know how you can trust the results of tfe in those frequency range...

[Mathieu NOE]

I may have to look again at your time signals (as done previously - see my former comments on that matter)

[Kumaresh Kumaresh]

Hello Mr. Mathieu,

Thank you for your comments. I am hereby attaching test.mat file.

1 column - time

2 column - input

all other columns - output, can be compared with input singal separately.

Kindly check it when you are free

Thank you

[Mathieu NOE]

hello again

I just looke d at your data : the input is a clean 200 Hz tone (not very long though), but the output(s) - I don't know - it looks just like random noise , you can't barelly distinguish a 200 Hz response in there , and that's confirmed by the fft spectrum .

so for me there is absolutely NO possibility to make any transfer function estimate with reasonnable quality as no information from the input is observed in the outputs (IMHO)

again - I believe I already spoke about that - to compute a valid transfer function in a given frequency range, you need the input signal has ebergy in that frequency range. So plotting a FRF estimate over a broad frequency range while the test signal is a single frequency sinewave is not good at all ! (sorry to be so direct !)

 load('test.mat')

dt = mean(diff(ans(:,1)));
Fs = 1/dt;

time = ans(:,1);
input = ans(:,2);
outputs = ans(:,3:end);

% input plot 
figure, 
subplot(211),plot(time,input);grid on
[f1,fft_spectrum1] = do_fft(time,input);
subplot(212),semilogx(f1,fft_spectrum1,'-*');grid on
title('Spectrum')
ylabel('FFT amp')
xlabel('Frequency[hz]')

% output plot 
figure, 
subplot(211),plot(time,outputs);grid on
for k = 1:size(outputs,2)
[f2,fft_spectrum2(:,k)] = do_fft(time,outputs(:,k));
end
subplot(212),semilogx(f2,fft_spectrum2,'-*');grid on
title('Spectrum')
ylabel('FFT amp')
xlabel('Frequency[hz]')

●

   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)

time = time(:);
data = data(:);

dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length


%% use windowing or not at your conveniance
% no window
fft_spectrum = abs(fft(data))*2/nfft;

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;

end

[Kumaresh Kumaresh]

Hello,

Thank you for your response. My work is related to computational simulation of aircraft gas turbine combustor. Based on the input conditions (in the form of sine wave, which is clean) applied, output distorted sine wave is obtained. Based on the physics, pressure and velocity are two critical parameters, which are examined.

I didn't apply numerator or denominator specifically like other traditional problems. No feedback loop analysis done for this open-loop data. As you mentioned, there might be no possibility to compute a valid transfer function in a given frequency range (as no information from the input is observed in the outputs)

With the given problem I have, I approached Bode Plot to understand the behavior in the frequency domain and examine the stability.

Kindly share your valuable insights

Thank you

[Mathieu NOE]

hello again

are you limited to sine inputs ? is this why we see f1 frequency tests carried at f1 = 0 / 100 / 200 Hz ?

you can build a Bde plot from discrete sine tests , but what worries me (again) is that in the example above, the 200 Hz sinewave (input) could not be seen in the spectrum of the output (as if the system has no gain in this frequency range)

can you try other frequencies below 200 Hz ?

all the best

[Kumaresh Kumaresh]

Sorry for replying late. Thank you for your comments.

Transfer function for pressure = p'_combustor / p'_inlet; where p' is pressure fluctuation.

Transfer function for velocity = u'_combustor / u'_inlet; where v' is velocity fluctuation.

So my gain is dimensionless.

Yes I am limited to sine inputs.

"200 Hz sinewave (input) could not be seen in the spectrum of the output" -- because my system (gas turbine combustor) is complex.