MATLAB: Deviation between tfestimate+lsim and bode on same transfer function

bodeCommunications ToolboxControl System ToolboxlsimSignal Processing Toolboxtfestimatetransfer 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

Best Answer

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.
Related Question