Why do I get unexpected fft result of filter output ?

3 次查看（过去 30 天）
K_S_
K_S_ 2022-7-29
评论： Paul 2022-8-1
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)|')
  3 个评论
显示 1更早的评论
Paul
Paul 2022-8-1
Ran in:
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by