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

40 Ansichten (letzte 30 Tage)
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('°')

Akzeptierte Antwort

Mathieu NOE
Mathieu NOE am 17 Jan. 2022
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
  5 Kommentare
Mathieu NOE
Mathieu NOE am 20 Jan. 2022
hello Ricky
hope you're doing well
do you mind accepting my answer ? tx

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by