Transfer function estimation from experimental data

64 views (last 30 days)
I have a particular system's frequency response data from experimental measurements and I would like to estimate a transfer function. Based on what I've read here, I have the magnitude and phase response response stored as an idfrd and I have tried playing around with the system identification toolbox (transfer function and state space models) without much luck. Nothing I've come up with really seems to fit it well. Am I missing something here?
Here's what the Bode plot looks like. Actual data is also attached.
  2 Comments
Star Strider
Star Strider on 1 Jul 2018
It may help to have your original, time-domain data. There doesn’t appear to be any relevant information in the idfrd data you posted. There are no poles, and I cannot even identify any zeros.
david
david on 2 Jul 2018
Don't have the original time domain data that corresponds to this set of idfrd data, but I took another data set using a chirped sinusoid. Sampled at 100 kHz. File has the system input and output response in time domain.

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 2 Jul 2018
This is the best I can do:
d = load('td_data.mat');
Fs = 100E+3; % Hz
tdu = d.input;
tdy = d.output;
L = numel(tdu);
t = linspace(0, L, L)/Fs;
figure
plot(t, tdu, t, tdy)
grid
Fn = Fs/2;
FTy = fft(tdy)/L;
FTu = fft(tdu)/L;
FTH = FTy./FTu; % Empirical Transfer Function
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
semilogy(Fv, abs(FTy(Iv))*2, Fv, abs(FTu(Iv))*2, Fv, abs(FTH(Iv))*2)
grid
xlim([0 5E+3])
legend('Output', 'Input', '‘Transfer Function’', 'Location','E')
figure
plot(Fv, imag(FTH(Iv))) % Poles, Zeros On Imaginary Axis
grid
xlim([0 5E+3])
figure
plot(Fv, abs(FTH(Iv))*2)
grid
xlim([0 5E+3])
D = iddata(tdy, tdu, 1/Fs);
sys = tfest(D, 5, 5);
figure
bode(sys)
An easy way to estimate the number of poles and zeros your system has is to plot the imaginary part of the transfer function as a function of frequency. The poles and zeros will be readily visible. I counted 5 poles before your system degenerated into noise. This is in figure(3).
I would incrementally add more poles and the equivalent number of zeros until it shows the desired high-pass characteristic. It has a zero at the origin, a pole at infinity, and a maximum between about 250-2000 Hz (in figure(4)), so I would experiment to get those characteristics in ‘sys’. I began to see that with 10 poles and 10 zeros.

Community Treasure Hunt

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

Start Hunting!