Why do I get unexpected fft result of filter output ?

To make sure the notch filter is working, I input a sine wave signal and performed frequency analysis of the output result to confirm the amount of attenuation.

I ran the attached code below and the simlink model.

As a result, since the gain of the filter was increased by 0.1 times at 50Hz, I thought the output amplitude would also increase by 0.1 times, but it did not.

Please tell me the reason.

I am a beginner in frequency analysis, so it may be due to lack of knowledge of frequency analysis, not Matlab, but I want to know.

%% Sampling setting
fs = 1000;                        % Sampling frequency [Hz]
Ts = 1/fs;
ntrans = 1000;                    % Transient response 
nsteady = 1000;                   % Steady-state response
nn = ntrans + nsteady;            
Tsim = nn*Ts;
t = (0:nn)'*Ts;
%% Input setting %%
ampli = 100;                       % sin wave amplitude
fn =50;                            % sin wave frequency
%% Notch filter setting %%
wn = 2*pi*fn;                       
zeta = 0.1;
d = 0.1;
b = [1 2*d*zeta*wn wn^2];
a = [1 2*zeta*wn wn^2];
H = tf(b,a);

figure(1)
title('Bode Plot of Notch Filter')
bode(H)
filename = 'test_simlink.slx'
open_system(filename)
out = sim(filename)
figure(2)
x = out.x;
plot(t,x)
title('Input Signal x(t)')
xlabel('t (sec)')
ylabel('x')
figure(3) 
y = out.y;
plot(t,y)
title('Output Signal y(t)')
xlabel('t (sec)')
ylabel('y')

%% FFT %%
x = x(ntrans+2:end,1);
X = fft(x);
L = length(x);
X = abs(X/L);
X = X(1:L/2+1);
X(2:end-1) = 2*X(2:end-1);
f = fs*(0:(L/2))/L;
figure(4)
plot(f,X) 
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('f (Hz)')
ylabel('|X(f)|')
y = y(ntrans+2:end,1);
Y = fft(y);
L = length(y);
Y = abs(Y/L);
Y = Y(1:L/2+1);
Y(2:end-1) = 2*Y(2:end-1);
f = fs*(0:(L/2))/L;
figure(5)
plot(f,Y) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('f (Hz)')
ylabel('|Y(f)|')

[Paul]

Hi K_S_,

How exactly is H implemented in the Simulink model?

Please edit the question to either attach or insert figures 4 and 5 into the question, via either the paperclip or image icon in the "Insert" section of the toolbar.

[K_S_]

Hi Paul,

I attached figures 4 and 5.

[Paul]

Hi KS

Let's open the figures

openfig(websave('fig4.fig','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1084040/fig4.fig'));

openfig(websave('fig5.fig','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1084045/fig5.fig'));

That does look odd for the stated filter.

Testing the code for the frequency domain analysis using a sine wave input, not the Simulink model

%% Sampling setting
fs = 1000;                        % Sampling frequency [Hz]
Ts = 1/fs;
ntrans = 1000;                    % Transient response 
nsteady = 1000;                   % Steady-state response
nn = ntrans + nsteady;            
Tsim = nn*Ts;
t = (0:nn)'*Ts;
%% Input setting %%
ampli = 100;                       % sin wave amplitude
fn =50;                            % sin wave frequency
%% Notch filter setting %%
wn = 2*pi*fn;                       
zeta = 0.1;
d = 0.1;
b = [1 2*d*zeta*wn wn^2];
a = [1 2*zeta*wn wn^2];
H = tf(b,a);

The filter gain at the sine wave frequency is:

m = bode(H,2*pi*fn)

m = 0.1000

So we'd expect the output amplitude (in steady state as done in the code) to be 1/10 of the input amplitude.

% figure(1)
% title('Bode Plot of Notch Filter')
% bode(H)
% filename = 'test_simlink.slx'
% open_system(filename)
% out = sim(filename)
% figure(2)
% x = out.x;
x = ampli*sin(2*pi*fn*t);
% plot(t,x)
% title('Input Signal x(t)')
% xlabel('t (sec)')
% ylabel('x')
% figure(3) 
% y = out.y;
y = lsim(H,x,t);
% plot(t,y)
% title('Output Signal y(t)')
% xlabel('t (sec)')
% ylabel('y')
%% FFT %%
x = x(ntrans+2:end,1);
X = fft(x);
L = length(x);
X = abs(X/L);
X = X(1:L/2+1);
X(2:end-1) = 2*X(2:end-1);
f = fs*(0:(L/2))/L;
figure(4)
plot(f,X) 
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('f (Hz)')
ylabel('|X(f)|')

y = y(ntrans+2:end,1);
Y = fft(y);
L = length(y);
Y = abs(Y/L);
Y = Y(1:L/2+1);
Y(2:end-1) = 2*Y(2:end-1);
f = fs*(0:(L/2))/L;
figure(5)
plot(f,Y) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('f (Hz)')
ylabel('|Y(f)|')

We see that the amplitude of the output is a bit larger than expected. But recall that lsim actually converts the continous-time system into a discrete-time approximation, so the actual expected gain is

m = bode(c2d(H,Ts,'foh'),2*pi*fn)

m = 0.1074

which is basically the gain we observe (100*0.0174 = 10.74).

The analysis code seems to be correct, so the problem must be somwhere else. I don't open Simulink models online and so can't comment on the implementation.

Is H implemented in a Transfer Function block?

What are the solver settings? Fixed step with a step size of 0.001?

How is the y(t) collected for output?