# The frequency response function (FRF) fails to detect the antiresonance of a system

I am trying to identify a vibrational systems by computing the frequency response function (FRF) of the system when a chirp signal is applied to its input. After comparing the FRF computed and the bode plot of the transfer function I found that the resonance has been well detected but the anti-resonance has not.

The system has the following form:

$$\textrm{sys} = \frac{s^2+\omega_{n_{z_1}}^2}{s^2+\omega_{n_{p_1}}^2}, \quad \textrm{where } \omega_{n_z} = 10\textrm{ Hz} \textrm{ and } \omega_{n_p} = 100\textrm{ Hz}.\\$$

I have a hypothesis:

The steady state response of a linear system to a sinusoidal signal is another sinusoidal signal with the same frequency with an amplitude and a phase that depend on the system. When I applied a chirp signal, the system never reached the steady state. So while the output of the system to $\sin(\omega_{n_{z_1}}\cdot 2\pi)$ should be zero, I am still obtaining some sinusoidal signal of frequency $100\textrm{ Hz}$ as a consequence of stimulating (weakly) that part of the spectrum with the transient response.

In my opinion the correct computation should be measuring only the components of the outputs that correspond to the frequency of the input at certain particular moment.

• Do you know how to correct this problem?
• Is it possible to apply a moving filter or window to correct this problem?

Here is my code:

% Identification of the Frequency Response of a Transfer Function with
% Resonance and Antiresonance peaks
%% Cleaning the house
clc;
clear;
close all;
%% Definition of the System
wn_z1 = 10*2*pi;           % Anti-Resonance Natural frequency 10 Hz
wn_p1 = 80*2*pi;           % Resonance Natural frequency 80 Hz
chi_z = 0;                 % Damping factor of the zeros
chi_p = 0;                % Damping factor of the poles
s     = tf('s');
sys_1 = (s^2+2*chi_z*wn_z1*s+ wn_z1^2)/(s^2+2*chi_p*wn_p1*s+wn_p1^2);
zpk(sys_1)
%% Definition of the Chirp Input
h                    = 0.001;
time                 = 0:0.001:60;
freq_ini             = 0.1*2*pi; % Initial frequency considered (0.1 Hz)
freq_final           = 240*2*pi; % Final frequency considered (240 Hz)
chirp_input          = 1.5*chirp(time,freq_ini,time(end),freq_final,'logarithmic');

frequency_evolution  = logspace(log10(freq_ini),log10(freq_final),length(time));
figure(1);
subplot(2,1,1);plot(time,chirp_input); xlabel('time (s)'); ylabel('Chirp    amplitude')
subplot(2,1,2);plot(time,frequency_evolution/(2*pi)); xlabel('time (s)'); ylabel('Frequency (Hz)');
%% Execution of the simulation
[sys_output,~]       = lsim(sys_1,chirp_input,time);
%% Computation of the FRF
fft_input   = fft(chirp_input,2^12);                                           %% Considering zero padding
fft_output  = fft(sys_output ,2^12);                                           %% Considering zero padding
fft_input   = fft_input(1:length(fft_input)/2);
fft_output  = fft_output(1:length(fft_output)/2);
fft_output  = fft_output(:);
fft_input   = fft_input(:);
FRF              = abs(fft_output)./abs(fft_input);
frequency_vector = linspace(0,1/(2*h),length(FRF));

%% Comparing the Bode Plot (Matlab) with the FRF (own implementation)
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
close all;
bodemag(sys_1,opts,'r'); hold on;
semilogx(frequency_vector, 20*log10(FRF ),'b')
xlim([0.1,frequency_vector(end)])
legend('Matlab Bode','FRF Identification')

• Can you please post plots of the input and output chirps? Is the system a physical system or a simulated one? Have you tried to plot on a semilogy? According to the chirp response, the system doesn't exhibit that zero at 10Hz, BUT! your chirp might be too fast. There could be a slight "kink" at the curve around 10Hz but the Y range is too wide for it to be seen. – A_A Jun 23 '16 at 15:51

I agree with A_A that the chirp is too fast.

I adjusted a few things in the code

• increase time to $600$ (seconds)
• start chirp at $0.5\textrm{ Hz}$
• end chirp at $2\textrm{ Hz}$
• removed zero padding (not required here)

And the result is...

The full code is

% Identification of the Frequency Response of a Transfer Function with
% Resonance and Antiresonance peaks
%% Cleaning the house
clc;
clear;
close all;
%% Definition of the System
wn_z1 = 10*2*pi;           % Anti-Resonance Natural frequency 10 Hz
wn_p1 = 80*2*pi;           % Resonance Natural frequency 80 Hz
chi_z = 0;                 % Damping factor of the zeros
chi_p = 0;                % Damping factor of the poles
s     = tf('s');
sys_1 = (s^2+2*chi_z*wn_z1*s+ wn_z1^2)/(s^2+2*chi_p*wn_p1*s+wn_p1^2);
zpk(sys_1)
%% Definition of the Chirp Input
h                    = 0.001;
time                 = 0:0.001:600;
freq_ini             = 0.5*2*pi; % Initial frequency considered (0.1 Hz)
freq_final           = 2*2*pi; % Final frequency considered (240 Hz)
chirp_input          = 1.5*chirp(time,freq_ini,time(end),freq_final,'logarithmic');

frequency_evolution  = logspace(log10(freq_ini),log10(freq_final),length(time));
figure(1);
subplot(2,1,1);plot(time,chirp_input); xlabel('time (s)'); ylabel('Chirp    amplitude')
subplot(2,1,2);plot(time,frequency_evolution/(2*pi)); xlabel('time (s)'); ylabel('Frequency (Hz)');
%% Execution of the simulation
[sys_output,~]       = lsim(sys_1,chirp_input,time);
%% Computation of the FRF
fft_input   = fft(chirp_input);                                           %% Considering zero padding
fft_output  = fft(sys_output );                                           %% Considering zero padding
fft_input   = fft_input(1:length(fft_input)/2);
fft_output  = fft_output(1:length(fft_output)/2);
fft_output  = fft_output(:);
fft_input   = fft_input(:);
FRF              = abs(fft_output)./abs(fft_input);
frequency_vector = linspace(0,1/(2*h),length(FRF));

%% Comparing the Bode Plot (Matlab) with the FRF (own implementation)
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
close all;
bodemag(sys_1,opts,'r'); hold on;
semilogx(frequency_vector, 20*log10(FRF ),'b')
xlim([0.1,frequency_vector(end)])
legend('Matlab Bode','FRF Identification')

• Thanks for your comment. I still have two question: Is it unnecessary to excite the system with a chirp signal with frequency components in those part of the spectrum that I want to identify? Why here zero padding is unnecesary? Is it not always advantageous practice? Regards, JC – Juan Carlos Ibáñez Jun 23 '16 at 17:11
• In general you can identify your system (reliably) only in the frequency region that is properly excited. In your simulated example this may seem different, as you have no noise. In reality, when you identify based on measured (=noisy) data, things can look different. About zero padding you find more in this link. My intention was to not bother with fft-lengths with the increased signal length (to be honest - I was to lazy to think about the length of the "new" signals) – snowflake Jun 24 '16 at 5:35

True that a chirp signal helps to get the FRF, but every time we change the frequency we can't reach the steady state, so this will cause bias in the estimation. As an advice try to use the multisine excitation, they are more suitable for such cases.