Why do I get unexpected fft result of filter output ?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Paul
am 1 Aug. 2022
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)
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)
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?
Antworten (0)
Siehe auch
Kategorien
Mehr zu Spectral Measurements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!