Designed filter shows half the cut-off frequency that it should
2 次查看(过去 30 天)
显示 更早的评论
hi everyone,
I tried to implement the designed filter (with fdesign), using its coefficients and test it with a chirp signal. here is my code:
clear screen
clearvars
format long
SamplingFrequency = 2000;
SamplingPeriod = 1/SamplingFrequency;
CenterFrequency = 450;
TransitionWidth = 100;
PassbandRipple = 1;
StopbandAttenuation = 80;
DesignSpec = fdesign.lowpass('Fp,Fst,Ap,Ast',...
CenterFrequency-TransitionWidth/2, ...
CenterFrequency+TransitionWidth/2, ...
PassbandRipple,StopbandAttenuation, ...
SamplingFrequency);
LPF = design(DesignSpec,'equiripple','SystemObject',true);
fvtool(LPF)
Coefficients = LPF.Numerator;
%%% here we define different input scenarios to be tested on the FIR filter
%%% Impulse Response
FIR_In = cat(2, 1, zeros(1,size(Coefficients, 2)));
FIR_Out = myFIR(FIR_In, Coefficients);
%%% Chirp Response
clear FIR_Out FIR_In
TotalTime = 100; %seconds
StartFrequency = 5; %Herz
StopFrequency = 1000; %Herz
t = 0:SamplingPeriod:TotalTime;
InstantanousFrequency = linspace(StartFrequency,StopFrequency,size(t,2));
disp (size(InstantanousFrequency))
x = zeros(1, size(t,2));
for i = 1:size(t,2)
x(i) = cos (2*pi*InstantanousFrequency(i)*t(i));
end
FIR_In = x;
FIR_In = cat(2, FIR_In, zeros(1,size(Coefficients, 2)));
FIR_Out = myFIR(FIR_In, Coefficients);
figure
plot (InstantanousFrequency,FIR_In(1,1:size(InstantanousFrequency,2)),InstantanousFrequency,FIR_Out(1,1:size(InstantanousFrequency,2)))
function FIR_Out = myFIR(FIR_In, Coefficients)
FIR_Out = zeros(1, size(FIR_In, 2));
x = zeros(1, size(Coefficients, 2));
for CurrentInputIndex = 1:1:size(FIR_In, 2)
%%% delay line
for DelayLineIndex = size(Coefficients, 2):-1:2
x(DelayLineIndex) = x(DelayLineIndex-1);
end
x(1) = FIR_In(CurrentInputIndex);
%%%
Acc = 0;
for TapIndex = 1:1:size(Coefficients, 2)
Acc = Acc + x(TapIndex) * Coefficients(TapIndex);
end
FIR_Out(CurrentInputIndex) = Acc;
end
end
as it is shown in the result, I get the transistion band on half the value that it should be (should be 400-500) and also the double passband. I think I am missing something about the nyquist criteria but I don't know what. Fs is 2000 and I am only chirping to 1000Hz. what am I doing wrong here?![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/909980/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/909980/image.jpeg)
0 个评论
回答(1 个)
Abhimenyu
2024-2-10
Hi Keyvan,
From the information shared by you, I could infer that you are trying to implement the designed filter using the "fdesign" function of MATLAB. The "Nyquist criterion" is correctly defined. The generation of the "chirp" signal is affecting the performance of the filter because the "InstantaneousFrequency" is not integrated to get the phase of the signal. The "cumtrapz" function of MATLAB can be used to integrate the "InstantaneousFrequency" to get the phase of the "chirp" signal, which is then used to generate the signal itself.
Please follow the example MATLAB R2023b code below to understand more:
%% NO CHANGE HAS BEEN MADE TO THE CODE IN THIS SECTION
clear screen
clearvars
format long
% Filter specifications
SamplingFrequency = 2000;
CenterFrequency = 450;
TransitionWidth = 100;
PassbandRipple = 1;
StopbandAttenuation = 80;
% Design the filter
DesignSpec = fdesign.lowpass('Fp,Fst,Ap,Ast',...
CenterFrequency - TransitionWidth/2, ...
CenterFrequency + TransitionWidth/2, ...
PassbandRipple, StopbandAttenuation, ...
SamplingFrequency);
LPF = design(DesignSpec, 'equiripple', 'SystemObject', true);
Coefficients = LPF.Numerator;
% Generate a chirp signal
TotalTime = 5; % seconds
StartFrequency = 5; % Hz
StopFrequency = 1000; % Hz
t = 0:1/SamplingFrequency:TotalTime;
InstantaneousFrequency = linspace(StartFrequency, StopFrequency, numel(t));
The "x" variable which is the "chirp" signal has been computed below using the MATLAB function "cumtrapz".
%% CHANGE MADE : GENERATION USING "cumtrapz" function
x = cos(2 * pi * cumtrapz(t, InstantaneousFrequency));
FIR_In = x;
% Apply the filter using the custom myFIR function
FIR_Out = myFIR(FIR_In, Coefficients);
figure
plot (InstantaneousFrequency,FIR_In(1,1:size(InstantaneousFrequency,2)),InstantaneousFrequency,FIR_Out(1,1:size(InstantaneousFrequency,2)))
%% NO CHANGE HAS BEEN MADE TO THE CODE IN THIS SECTION
function FIR_Out = myFIR(FIR_In, Coefficients)
FIR_Out = zeros(1, size(FIR_In, 2)); % Initialize the output signal
x = zeros(1, size(Coefficients, 2)); % Buffer for the filter's delay line
% Process each sample of the input signal
for CurrentInputIndex = 1:size(FIR_In, 2)
% Shift the delay line to make room for the new sample
for DelayLineIndex = size(Coefficients, 2):-1:2
x(DelayLineIndex) = x(DelayLineIndex-1);
end
x(1) = FIR_In(CurrentInputIndex); % Insert the new sample at the beginning
% Perform the dot product of the delay line and the filter coefficients
Acc = 0;
for TapIndex = 1:size(Coefficients, 2)
Acc = Acc + x(TapIndex) * Coefficients(TapIndex);
end
FIR_Out(CurrentInputIndex) = Acc; % Store the result in the output signal
end
end
For more information on the "cumtrapz" function, follow this MATLAB R2023b documentation link: https://www.mathworks.com/help/matlab/ref/cumtrapz.html
Regards,
Abhimenyu
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Filter Design 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!