# Why do I get unexpected fft result of filter output ?

2 views (last 30 days)
K_S_ on 29 Jul 2022
Commented: Paul on 1 Aug 2022
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
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)
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 on 1 Aug 2022
Hi KS
Let's open the figures
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
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?

### Categories

Find more on Array and Matrix Mathematics in Help Center and File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!