Deviation between tfestimate+lsim and bode on same transfer function

Goal:

Create an input signal to find the frequency response of a physical system.

Process:

Since I have an idea/hypothesis of what the frequency response of this physical system looks like I made a dummy transfer function (double mass spring damper) on which I could test the desired input signal. From this transfer function I made two bode plots: one using the "bode" function and one using the "tfestimate" function (in combination with the "lsim" function and my input signal).

Problem:

The course of both plots is extremely similar in terms of magnitude and phase values however there is quite a big difference in the frequencies at which similar values occur.

I tried using white noise instead of the input signal I desired, which gave the same results.

Code:

%% Parameters for Transfer Function
m0 = 0.0;
k0 = 0;
m1 = 3.6;  % Apperent stem mass at current robot position (not plant mass!) [kg]
m2 = 3.6;   % Felt mass of the branch (+ the remaining stem influence) [kg]
b1 = 10;     % Damping stem at current robot position [Ns/m]
b2 = 10;     % Damping of the branch (+ the remaining stem influence) [Ns/m]
k1= 400;    % Stiffness stem at current robot position [N/m]
k2 = 400;   % Stiffness of the branch (+ the remaining stem influence) [N/m]
x1 = 0; %0.05;   %Preloading of stem spring [m]

%% Analytical TF (double mass spring damper)
num = [m2, b2, k2];
den = [(m0+m1)*m2, (m0+m1)*b2+m2*b1+m2*b2, m2*k2+(m0+m1)*k2+m2*(k0+k1)+b1*b2 b1*k2+b2*(k0+k1), (k0+k1)*k2];
transfer = tf(num,den);
figure()
bode(transfer)

%% Input signal TF's
Fs = 100;
Ts = 1/Fs;
t = 60;
time = 0:Ts:t-Ts;

% white noise
in = wgn(length(time),1,10);

% % Frequency sweep (desired signal)
% bias = 0.0;
% f0 = 0.0;                                 %% Initial frequency (Hz)
% f1 = 10.0;                                %% End frequency (Hz)
% Amp0 =25.0;                              %% Initial amplitude (mm)
% Amp1 = 5.0;                               %% End amplitude (mm)
% target_time = 60.0;                       %% Time in which to get from initial to final frequency (sec)
% k = (f1-f0)/target_time;                  %% The frequency change rate (Hz/s)
% l = (Amp0-Amp1)/target_time;              %% The amplitude change rate (mm/s)
% for t=1:length(time) 
%     in(t) = (Amp0-l*time(t))*sin(2*pi*((k*time(t))/2+f0)*time(t)) + bias; 
% end

%% TF simulation
out = lsim(transfer, in, time);

% Enhancing spectral estimation for bode
res = 0.1;
nfft = Fs/res;
window = hann(nfft);
noverlap = nfft/2;
[H,hz] = tfestimate(in, out,window,noverlap,nfft,Fs);

% Bode
figure()
subplot(2,1,1)
semilogx(hz,mag2db(abs(H)))
title('Bode plot from white noise + lsim + tfestimate')
ylabel('Amplitude [dB]')
grid
subplot(2,1,2)
semilogx(hz,rad2deg(angle(H)))
ylabel('Phase [deg]')
xlabel('Frequency [rad/s]');
grid

[Anish Mitra]

Check the units of frequency in both cases.

[Lars Janssen]

Good thinking! I thought I already ruled this out, but in hindsight I should have seen the factor of 2*pi...

For anybody else who finds this thread, the following has happened: when using "tfestimate" the unit in which the frequency is outputted can change based on the input that is given. Here you can find that if you use [txy,w] = tfestimate(___) the "w" variable will be in rad/sample (when I saw "rad" I mistakenly assumed everything was the same as in using the "bode" function), however when you use [txy,f] = tfestimate(___,fs) (which I did) then the "f" variable will be in cycles/unit time, where fs specifies the unit time.

I solved this problem by simply multiplying "f" with 2*pi, but you can also get rid of the "fs" parameter in the "tfestimate" function and multiply "w" by "fs" to get the same result.

Thank you for the help @Anish Mitra. What do I do now? Since there is no answers I cannot accept an answer. Just close the question?

[Ameer Hamza]

Lars, there is no need to close the question. I think that the discussion you added in the comment is quite useful and will be helpful for others searching for a similar question. I recommend that you add this discussion in an answer and accept your own answer.

[Lars Janssen]

Will do! Thanks again for the help, it's truly appreciated!

Answer 1

[Lars Janssen]

For anybody else who finds this thread, the following has happened: when using "tfestimate" the unit in which the frequency is outputted can change based on the input that is given. Here you can find that if you use [txy,w] = tfestimate(___) the "w" variable will be in rad/sample, however when you use [txy,f] = tfestimate(___,fs) then the "f" variable will be in cycles/unit time, where fs specifies the unit time.

I solved this problem by simply multiplying "f" with 2*pi, but you can also get rid of the "fs" parameter in the "tfestimate" function and multiply "w" by "fs" to get the same result.