find transfer function of input and output signals and then perform an FFT.

45 views (last 30 days)
Hello,
I have collected two signals as input x(t) and output y(t) to the system. I want to find the transfer function(which is basically Output/Input as per my knowledge). But some of the values in the input signal are zero so I think performing this division would give 'divide by zero' error.
I tried using 'tf', 'tfestimate' and 'iddata' on these signals but I am not getting the right plot of the transfer function.
Ultimately after finding the transfer function, I want to perform fft on it and plot a single-sided amplitude spectrum(TF vs freq plot).
Can someone please help me with this? I am new to the Control systems and signal processing concepts so it will be really useful for me to understand the basics.
Thanks.

Accepted Answer

Star Strider
Star Strider on 20 Nov 2019
This should get you started:
[D,S] = xlsread('Test_data_1.xlsx');
torg = D(:,1); % Time Vector
inp = D(:,2); % Input Vector
out = D(:,3); % Output Vector
[Tu,~,idx] = unique(torg); % Unique Times & Index
inpm = accumarray(idx, inp, [], @mean); % Input Vector (Means Of Duplicates)
outm = accumarray(idx, out, [], @mean); % Output Vector (Means Of Duplicates)
treg = linspace(min(Tu), max(Tu), numel(Tu)).' - min(Tu); % Uniform Time Vector (Corrected To Begin At Zero)
Ts = mean(diff(treg)); % Sampling Interval
Fs = 1/Ts; % Sampling Frequency
inpr = interp1(Tu-min(Tu), inpm, treg); % Interpolated Input Data (To ‘treg’)
outr = interp1(Tu-min(Tu), outm, treg); % Interpolated Output Data (To ‘treg’)
L = numel(treg);
figure
plot(treg, inpr, treg, outr)
grid
xlabel('Time')
ylabel('Amplitude')
legend('Input', 'Output')
Fn = Fs/2;
inprm = inpr-mean(inpr);
outrm = outr-mean(outr);
FTinpr = fft(inprm)/L; % Fourier Transform Of Mean-Corrected Data
FToutr = fft(outrm)/L; % Fourier Transform Of Mean-Corrected Data
Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector (One-Sided)
Iv = 1:numel(Fv); % Index Vector
figure
subplot(3,1,1)
plot(Fv, abs(FTinpr(Iv))*2)
grid
title('Input Spectrum')
subplot(3,1,2)
plot(Fv, abs(FToutr(Iv))*2)
grid
title('Output Spectrum')
subplot(3,1,3)
plot(Fv, abs(FToutr(Iv)./FTinpr(Iv)))
smf = sgolayfilt(abs(FToutr(Iv)./FTinpr(Iv)), 3, 75);
hold on
plot(Fv, smf, 'LineWidth',2)
hold off
grid
title('Transfer Function Spectrum')
legend('Original','Savitzsky-Golay Filtered')
dm = iddata(outm, inpm, Ts); % Create ‘iddata’ Object
sys = tfest(dm, 3);
figure
bodeplot(sys);
grid
Your data had a significant number of duplicate values, so I used accumarray to take the mean of the duplicates, creating the ‘inpm’ and ‘outm’ vectors. I then used these in the rest of the analysis.
You will need to experiment with the numbers of poles and zeros to use with tfest, to correcttly approximate the transfer function.
  4 Comments
Giulia Lattanzi
Giulia Lattanzi on 18 Aug 2022
Hi, can I ask please what's the best approach to determine a right number of poles?
Star Strider
Star Strider on 18 Aug 2022
@Giulia Lattanzi — The way I generally determine them it is to take the fft of the transfer function and then plot only the imaginary part as a function of frequency. The poles (and their frequencies) as well as the zeros (and their frequencies) should readily reveal themselves. Be sure to note whether there are any of either at zero or infinity. This should give you a general idea.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!