Bode plot obtained from measured data is different from that obtained from transfer function

Hi All. I want obatain a bode diagram of a system without using bode function from Matlab, so I tried to get it by using time domain input and output data of that system. I implemented code in Matlab R2020a as we could see at below. However, the bode plot obtained from measured data is different from that obtained from transfer function as we could see the figure at below. Could any one tell me where I went wrong?

%% pre-operation
clear
close all
clc

%% define transfer function
num = [1 1 1000];
den = [1 0.1 20 10];
gs = tf(num,den);

%% generate excitation signal(chirp type)
ts = 0.001;             % sampling time step
tMax = 2;               % max time
t = 0:ts:(tMax-ts);     % time sequence

a = 1;                  % signal amplitude     
beta = 250;             % signal change vel
fIni = 0;               % initial freq
fEnd = beta*t(end);     % target freq
in = a*cos(2*pi*(beta*t+fIni).*t);  % excitation signal(chirp type)

freq = beta*t;          % freq sequence

%% plot excitation signal
figure
hold on
plot(t,in)
hold off
legend
grid on
box on
xlabel('time/s')
ylabel('counts')
title('Excitation signal in time domain')

%% get output data
out = lsim(gs,in,t)';

%% plot time domain input and output signal
figure
hold on
plot(t,in)
plot(t,out)
hold off
legend
grid on
box on
xlabel('time/s')
ylabel('counts')
title('Input and output in time domain from simulink')

%% plot bode diagram from transfer function
figure
bodeOp = bodeoptions;
bodeOp.FreqUnits = 'Hz';
bodemag(gs,bodeOp)
hold on

%% plot bode diagram from measured data
xSimInputfft = fft(in);
ySimOutputfft = fft(out);
gjw = ySimOutputfft./xSimInputfft;
mag = 20*log10(abs(gjw));
semilogx(freq,mag,'-o')
grid on
legend('bode from tf','bode form measured data')
xlim([0.1,1000])

Result

[Yang Philip]

Hi All. After get some hints from Jon and other references. I adjusted my code like below and it seems works better now.

However, I still have some confusions:

1. Why should I reduce sampling time step ts and signal frequency change vel beta, expand max time tMax to get better results (It seems the smaller ts and beta, bigger tMax could lead better result)?
2. How could I get bode plot for an unstable system through experiment method which means I didn't know transfer function of the system at first?

Maybe any one could give me some hints.

Code:

%% pre-operation
clear
close all
clc

%% define transfer function
num = [1 1];
den = [1 1 10];
gs = tf(num,den);

if ~isstable(gs)
    disp('Err: system unstable!')
    return
end

%% generate excitation signal(chirp type)
ts = 0.00025;           % sampling time step
tMax = 5;               % max time
t = 0:ts:(tMax-ts);     % time sequence

a = 1;                  % signal amplitude     
beta = 50;              % signal change vel
fIni = 0;               % initial freq
fEnd = beta*t(end);     % target freq
in = a*cos(2*pi*(beta*t+fIni).*t);  % excitation signal(chirp type)

freq = beta*t;          % freq sequence

%% plot excitation signal
figure
hold on
plot(t,in)
hold off
legend
grid on
box on
xlabel('time/s')
ylabel('counts')
title('Excitation signal in time domain')

%% get output data
out = lsim(gs,in,t)';

%% plot time domain input and output signal
figure
hold on
plot(t,in)
plot(t,out)
hold off
legend
grid on
box on
xlabel('time/s')
ylabel('counts')
title('Input and output in time domain from simulink')

%% bode diagram data calculation
input = in;
output = out;
time = t;

Fs = 1/ts;                          % 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

%% plot bode from tf and measured data
figure
hold on
bodeOp = bodeoptions;
bodeOp.FreqUnits = 'Hz';
bodemag(gs,bodeOp)
semilogx(Fv, 20*log10(abs(TF(Iv))),'-.','linewidth',2)
legend('bode plot from tf','bode plot from measured data')
grid on
xlim([0.1,1000])

Result:

Reference:

1. https://in.mathworks.com/matlabcentral/answers/375211-calculate-frequency-response-from-two-signal-on-domain-time-in-matlab;
2. https://in.mathworks.com/matlabcentral/answers/496878-is-it-possible-that-ploting-bode-diagram-without-transfer-fuction.

Answer 1

[Jon]

There are a number of problems with your problem formulation.

I think the primary one may be that your system, as defined, is unstable.

You can calculate roots(den) and see that it has poles in the right half plane.

Did you intend for the system to be unstable? I would have to dig into it more but my intuition is that identifying frequency responses of unstable systems using a chirp may be problematic.

Also, if you are looking for a chirp with a linearly increasing frequency I think you should define it as:

in = a*cos(2*pi*(beta*t+fIni));  % excitation signal(chirp type)

rather than

in = a*cos(2*pi*(beta*t+fIni).*t);  % excitation signal(chirp type)

which results in the time value being squared.

[Yang Philip]

Thanks for your advice Jon! I changed my plant to a stable system then I get better result. But, I think the chirp signal expression in my code does not wrong after I do some research.

[Jon]

Hi,

I'm glad that you were able to get better results once you used a stable system.

I think in general to perform an experimental system identification of an unstable system you would have to include it in an overal closed loop that would stabilize it. There are many subtelties to performing this type of analysis due the resulting correlation between the input and output signals.

I was thinking that you wanted a simple swept cosine with linearly increasing frequency. I think maybe what you have is some form of quadratic chirp, but I am not familiar with the details. You should however be careful that the sampling frequency (inverse of your sample period) is at least twice the highest frequency in your chirp input signal (Nyquist Criteria) avoid aliasing.

If you think I have answered your question it would be good to mark it as answered when you have a chance so that others will know an answer is available.