Designing a fir filter and how to apply it to a signal

52 ビュー (過去 30 日間)
Nediljko Ivisic
Nediljko Ivisic 2020 年 7 月 2 日
コメント済み: Star Strider 2020 年 7 月 2 日
Hello, I have to design a bandstop fir filter for a college project using hamming window with passband frequencys of 0.15 and 0.3 and apply it to a signal x = sin(2*pi*0.2*t) + sin(2*pi*0.1*t), 0<t<128. So far I wrote this:
fs = 10;
fp1 = 0.15;
fp2 = 0.3;
f1 = 0.2;
f2 = 0.1;
t = 128;
n = 0:1/fs:t;
x = sin(2*pi*f1*n) + sin(2*pi*f2*n);
Wn=[fp1 fp2];
Ham = fir1(30, Wn, 'stop');
[H,w] = freqz(Ham);
y = filter(Ham,1,x);
dB = mag2db(abs(H));
figure(1)
plot(w/pi, dB);
xlabel('Frequency');
ylabel('Magnitude (dB)');
figure(2)
plot(n, y);
title("Output y(n)",'fontsize',14);
xlabel('n');
ylabel('Y(n)');
However, my output signal y is almost the same as my input signal x, when I believe all that should be left is sin(2*pi*0.1*t). If I use bandstop function as "bandstop(x,Wn,fs)" the results are fine, but I cant replicate them using the upper method. If anyone knows what I am doing wrong I would appreciate any tips and advices.

採用された回答

Star Strider
Star Strider 2020 年 7 月 2 日
You need to normalise the stopband frequencies by the Nyquist frequency.
Try this:
Wn=[fp1 fp2]/(fs/2);
That worked when I tried it.
  2 件のコメント
Nediljko Ivisic
Nediljko Ivisic 2020 年 7 月 2 日
Seems to work the trick and my output is fine now, but the plot of filter magnitude looks a bit weird. I cant read out the passband frequencys from it now. Am I ploting it wrong?
Star Strider
Star Strider 2020 年 7 月 2 日
The freqz call needs to change a bit. You also need to increase the filter order, since normalising by the Nyquist frequency changes (decreases) the passband frequencies, and use filtfilt to do the actual filtering.
Change to these:
Wn=[fp1 fp2]/(fs/2);
Ham = fir1(90, Wn, 'stop');
[H,w] = freqz(Ham,1,2^16,fs);
and this:
figure(1)
plot(w, dB);
and this:
y = filtfilt(Ham,1,x);
Experiment with different filter orders to see those effects.
.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeKaiser についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by