FFT doesn't give desired frequency response

Hi there,

I am currently trying to evaluate the frequency response of a system, defined as: TargetForce --> System --> ActualForce. I have a timetable (TimetableData.txt) containing both channels, evenly resampled. My function [FFT_struct] = FFT_from_tt(tt) is supposed to give the single-sided and padded fft of both channels. I have to problems/misunderstandings with the result:

• The amplitudes don't match up with the time-domain amplitudes (see "timedomain.png" vs. "targetforce_frequencydomain_amplitude.png"). I've tried the function with one of Matlab's code examples and it does work there, which gives me a hint that something is wrong with my data - but I have no clue what that could be.
• The phase-response looks strange (see "frequencyResponse_Phase.png"). I obtained it by the code below (see end of code snippet). I would have expected it to be the "inverse" and something between 0 and -pi.

As I said, I suggest that I treat my specific data wrong - but I really need help with finding out. Thanks in advance!

function [FFT_struct] = FFT_from_tt(tt)
% FFT_FROM_TT Calculates FFT columnwise from evenly resampled timetable.
% Usage of padding 

% Tolerance for fft analysis - UNUSED HERE
tol = 0; 

sample_time = seconds(mean(diff(tt.Time)));
sample_rate = 1/sample_time;
L = size(tt, 1);
NFFT=2^nextpow2(L);
chns = size(tt, 2);
freq = sample_rate/2*linspace(0,1,NFFT/2+1);

%% FFT 
for ch = 1:chns
    Y = fft(fillmissing(table2array(tt(:, ch)), 'linear'), NFFT)/L;  % Double-sided transformation
    Y = 2*Y(1:NFFT/2+1);  % Single-sided transformation
    Y(abs(Y)<tol) = 0; % Tolerance for frequency response analysis

    FFT_struct.(tt.Properties.VariableNames{ch}).freq = freq;
    FFT_struct.(tt.Properties.VariableNames{ch}).Y = Y;
    FFT_struct.(tt.Properties.VariableNames{ch}).ampl = abs(Y); 
    FFT_struct.(tt.Properties.VariableNames{ch}).phase = unwrap(angle(Y));  % [rad]   

end



end




%% NOT IN THIS FUNCTION - Calculation of phase response
YTF = YActual./YTarget;
plot(freq, unwrap(angle(YTF)))

Answer 1

[Mathieu NOE]

hello

there are different estimators used to compute a transfer function

using just the ratio of two FFT is doable but does not provide the best result

here (below) , we are using the cross spectrum divided by the input auto spectrum method (called H1 estimator if I recall correctly)

You can notice that the validity of your data does not go much beyong 3 Hz , as the coherence drops very rapidely above this limit. So acquiring data at 400 Hz is quite fast for a system where the effective bandwithh seems to be only 5 Hz.

I did'nt do it here , but you can either decimate the data before doing the fft , or simply on the acquisition system, reduce the sampling rate ( 100 hz would largely suffice)

Your txt file I modified it (remove the text "sec" in the first column) so I could easily load the data with regular readmatrix

here's your transfer function (Bode plot) - the lower plot is your coherence between the two channels

you see the mag and phase plots indicates the actual force follows the target force between 0 and 3 Hz , then you have a low pass roll of characteristic (normal !) , and above 3 hz the quality of the signals are not very good (noise and noise )

unzip("TimetableData.zip")

data = readmatrix('TimetableData.txt',"NumHeaderLines",10); % Time,TargetForce,ActualForce
t = data(:,1);
x = data(:,2); %input  : TargetForce
y = data(:,3); %output : ActualForce
dt = mean(diff(t));
fs = 1/dt;   % sampling frequency

NFFT = 2048*2;
NOVERLAP = round(0.75*NFFT);  % 75 percent overlap

%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);

%% alternative with supplied function 
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP); 
% Txy = transfer function (complex), Cxy = coherence, F = freq vector

% Bode plots
figure(1),
subplot(3,1,1),semilogx(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),semilogx(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),semilogx(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');

%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_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]

as mentionned above , the sampling rate is a bit fast for the required bandwith

here I decimated the signals by a factor 20 to reduce the freq range and FFT computation burden by the same time

unzip("TimetableData.zip")
data = readmatrix('TimetableData.txt',"NumHeaderLines",10); % Time,TargetForce,ActualForce
t = data(:,1);
x = data(:,2); %input  : TargetForce
y = data(:,3); %output : ActualForce

% rmove NaN values 
ind = ~isnan(t) & ~isnan(x) & ~isnan(y);
t = t(ind);
x = x(ind); %input  : TargetForce
y = y(ind); %output : ActualForce

dt = mean(diff(t));
fs = 1/dt;   % sampling frequency

% decimate as we don't need the full frequency range up to Nyquist (200 Hz)
r = 20;
x = decimate(x,r);
y = decimate(y,r);
fs = fs/r;

NFFT = 256;
NOVERLAP = round(0.75*NFFT);  % 75 percent overlap
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied function 
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP); 
% Txy = transfer function (complex), Cxy = coherence, F = freq vector
% Bode plots
figure(1),
subplot(3,1,1),semilogx(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),semilogx(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),semilogx(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');

%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_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

problem solved ?

[Klafo]

Sorry for the late answer, works! Thank you very much sir!

[Mathieu NOE]

as always, my pleasure !