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