完整版教程下載地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547
第36章 FIR濾波器的Matlab設計(含低通，高通，帶通和帶阻)
本章節講解FIR濾波器的Matlab設計。主要是函數fir1和fir2的使用。
36.1 窗函數
36.2 fir1函數
36.2 fir2函數
36.4 總結
36.1 窗函數
在數字信號處理中不可避免地要用到數據截取的問題。例如，在應用DFT的時候，數據x(n)總是有限長的，在濾波器設計中遇到了對理想濾波器抽樣響應h(n)的截取問題，在功率譜估計中也要遇到對自相關函數的截取問題。總之，我們在實際工作中所能處理的離散序列總是有限長，把一個長序列變換成有限長的序列不可避免的要用到窗函數。因此，窗函數本身的研究及其應用是信號處理中的一個基本問題。
不同的窗函數對信號頻譜的影響是不一樣的，這主要是因為不同的窗函數，產生泄漏的大小不一樣，頻率分辨能力也不一樣。信號的截斷產生了能量泄漏，而用FFT算法計算頻譜又產生了柵欄效應，從原理上講這兩種誤差都是不能消除的，但是我們可以通過選擇不同的窗函數對它們的影響進行抑制。(矩形窗主瓣窄，旁瓣大，頻率識別精度最高，幅值識別精度最低；布萊克曼窗主瓣寬，旁瓣小，頻率識別精度最低，但幅值識別精度最高)。
對於窗函數的選擇，應考慮被分析信號的性質與處理要求。如果僅要求精確讀出主瓣頻率，而不考慮幅值精度，則可選用主瓣寬度比較窄而便於分辨的矩形窗，例如測量物體的自振頻率等；如果分析窄帶信號，且有較強的干擾噪聲，則應選用旁瓣幅度小的窗函數，如漢寧窗、三角窗等；對於隨時間按指數衰減的函數，可采用指數窗來提高信噪比。
- 矩形窗：
矩形窗屬於時間變量的零次冪窗。矩形窗使用最多，習慣上不加窗就是使信號通過了矩形窗。這種窗的優點是主瓣比較集中，缺點是旁瓣較高，並有負旁瓣，導致變換中帶進了高頻干擾和泄漏，甚至出現負譜現象。
- 三角窗：
三角窗亦稱費傑（Fejer）窗，是冪窗的一次方形式。與矩形窗比較，主瓣寬約等於矩形窗的兩倍，但旁瓣小，而且無負旁瓣。
- 漢寧窗：
漢寧窗又稱升余弦窗，漢寧窗可以看作是3個矩形時間窗的頻譜之和，或者說是 3個 sinc（t）型函數之和，而括號中的兩項相對於第一個譜窗向左、右各移動了 π/T，從而使旁瓣互相抵消，消去高頻干擾和漏能。可以看出，漢寧窗主瓣加寬並降低，旁瓣則顯著減小，從減小泄漏觀點出發，漢寧窗優於矩形窗．但漢寧窗主瓣加寬，相當於分析帶寬加寬，頻率分辨力下降。
- 海明窗：
海明窗也是余弦窗的一種，又稱改進的升余弦窗。海明窗與漢寧窗都是余弦窗，只是加權系數不同。海明窗加權的系數能使旁瓣達到更小。分析表明，海明窗的第一旁瓣衰減為一42dB．海明窗的頻譜也是由3個矩形時窗的頻譜合成，但其旁瓣衰減速度為20dB/（10oct），這比漢寧窗衰減速度慢。海明窗與漢寧窗都是很有用的窗函數。
- 高斯窗：
三角窗亦稱費傑（Fejer）窗，是冪窗的一次方形式。與矩形窗比較，主瓣寬約等於矩形窗的兩倍，但旁瓣小，而且無負旁瓣。
還有很多其它的窗口這里就不做介紹了，需更詳細的了解的話，可以看matlab中help文檔中的如下部分：
或者直接在命令窗口輸入windowDesigner可以打開窗口工具：
打開后界面如下：
36.2 fir1函數
36.2.1 fir1函數介紹
函數fir1用來設計標准頻率響應的基於窗函數的FIR濾波器，可實現加窗線性相位FIR濾波器設計。
語法：
b = fir1(n,Wn)
b = fir1(n,Wn,'ftype')
b = fir1(n,Wn,window)
b = fir1(n,Wn,'ftype',window)
b = fir1(...,'normalization')
其中，n：為了濾波器的階數；
Wn：為濾波器的截止頻率；
ftype：參數用來決定濾波器的類型，當ftype=high時，可設計高通濾波器，當ftype=stop時，可設計帶阻濾波器。Window參數用來指導濾波器采用的窗函數類型。其默認值為漢明（Hamming）窗。
使用fir1函數可設計標准的低通，高通，帶通和帶阻濾波器。濾波器的系數包含在返回值b中，可表示為：
b(z) = b(1) + b(2)z-1 + …… +b(n+1)z-n
(1) 采用漢明窗設計低通FIR濾波器
使用b=fir1(n, Wn)可得到低通濾波器。其中， 0<=Wn<=1, Wn=1相當於0.5。其語法格式為
b=fir1(n, Wn)
(2) 采用漢明窗設計高通FIR濾波器
在b=fir1(n, Wn, 'ftype')中，當ftype=high時，可設計高通濾波器。其語法格式為
b=fir1(n, Wn, 'high')
(3) 采用漢明窗設計帶通FIR濾波器
在b=fir1(n, Wn)中，當Wn=[W1 W2]時，fir1函數可得到帶通濾波器，其通帶為W1 < W < W2
W1 和 W2分別為通帶的下限頻率和上限頻率。其語法格式為
b=fir1(n, [W1 W2])
(4) 采用漢明窗設計帶阻FIR濾波器
在b = fir1(n,Wn,'ftype')中，當ftype=stop，Wn=[W1 W2]時，fir1函數可得到帶阻濾波器，其語法格式為
b=fir1(n, [W1 W2], 'stop')
(5) 采用其他窗口函數設計FIR濾波器
使用Window參數，可以用其他窗口函數設計出各種加窗濾波器，Window參數可采用的窗口函數有Boxcar，Hamming，Bartlett，Blackman，Kasier和Chebwin等。其默認時為Hamming窗。例如，采用Bartlett窗設計帶阻濾波器，其語法結構為
b=fir1(n, [W1 W2], 'stop', Bartlett[n+1])
注意：用fir1函數設計高通和帶阻濾波器時，所使用的階數n應為偶數，當輸入的階數n為奇數時，fir1函數會自動將階數增加1形成偶數。
36.2.2 fir1設計低通濾波器實例
下面我們通過一個實例來講解fir1的用法。原始信號是由50Hz正弦波和200Hz的正弦波組成，將200Hz的正弦波當做噪聲濾掉，下面通過函數fir1設計一組低通濾波器系數，其階數是30，截止頻率為0.25（也就是125Hz）。Matlab運行代碼如下：
%**************************************************************************************** % FIR低通濾波器設計 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*50*t); %信號50Hz正弦波 Signal_Noise=sin(2*pi*200*t); %噪聲200Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; b = fir1(30, 0.25); %30階FIR低通濾波器，截止頻率125Hz %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab的運行結果如下：
從運行結果的FFT和信噪比來看，濾波效果比較明顯。
36.2.3 fir1設計高通濾波器實例
下面我們通過一個實例來講解fir1的高通濾波器的用法。原始信號是由50Hz正弦波和200Hz的正弦波組成，將50Hz的正弦波當做噪聲濾掉，下面通過函數fir1設計一組高通濾波器系數，其階數是30，截止頻率為0.25（也就是125Hz）。Matlab運行代碼如下：
%**************************************************************************************** % FIR高通濾波器設計 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*200*t); %信號200Hz正弦波 Signal_Noise=sin(2*pi*50*t); %噪聲50Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; b = fir1(30, 0.25, 'high'); %30階FIR低通濾波器，截止頻率125Hz %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab的運行結果如下：
從運行結果的FFT和信噪比來看，濾波效果比較明顯。
36.2.4 fir1設計帶通濾波器實例
下面我們通過一個實例來講解fir1的帶通濾波器的用法。原始信號是由50Hz正弦波和200Hz的正弦波組成，設計通帶為125Hz到300Hz，下面通過函數fir1設計一組帶通濾波器系數，其階數是30，通帶為0.25 < W <0.6。Matlab運行代碼如下：
%**************************************************************************************** % FIR帶通濾波器設計 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*200*t); %信號200Hz正弦波 Signal_Noise=sin(2*pi*50*t); %噪聲50Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; b = fir1(30, [0.25 0.6]); %30階FIR低通濾波器，截止頻率125Hz %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab運行結果如下：
從運行結果的FFT和信噪比來看，濾波效果比較明顯。
36.2.5 fir1設計帶阻濾波器實例
下面我們通過一個實例來講解fir1的帶阻濾波器的用法。原始信號是由50Hz正弦波和200Hz的正弦波組成，設計阻帶為125Hz到300Hz，下面通過函數fir1設計一組帶阻濾波器系數，其階數是30，阻帶為0.25 < W <0.6。Matlab運行代碼如下：
%**************************************************************************************** % FIR帶阻濾波器設計 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*50*t); %信號50Hz正弦波 Signal_Noise=sin(2*pi*200*t); %噪聲200Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; b = fir1(30, [0.25 0.6], 'stop'); %30階FIR低通濾波器，截止頻率125Hz %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab運行結果如下：
從運行結果的FFT和信噪比來看，濾波效果比較明顯。
36.2.6 切比雪夫窗口函數設計帶通濾波器實例
下面我們通過一個實例來講解fir1設計切比雪夫窗口的的帶通濾波器。原始信號是由50Hz正弦波和200Hz的正弦波組成，設計通帶為125Hz到300Hz，下面通過函數fir1設計一組帶通濾波器系數，其階數是30，通帶為0.25 < W <0.6，並且具有25db波紋的切比雪夫窗。Matlab運行代碼如下：
%**************************************************************************************** % 切比雪夫窗口函數設計帶通濾波器 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*200*t); %信號200Hz正弦波 Signal_Noise=sin(2*pi*50*t); %噪聲50Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; Window = chebwin(31, 25); %25db的切比雪夫窗 b = fir1(30, [0.25 0.6], Window); %30階FIR低通濾波器，截止頻率125Hz %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab運行結果如下：
通過歸一化頻率可以看出切比雪夫窗口是有一定紋波的。不過從FFT結果和信噪比來看，通過切比雪夫窗口做的濾波效果也是比較明顯的。
36.3 fir2函數
36.3.1 fir2函數介紹
函數fir2用來設計有任意頻率響應的各種加窗FIR濾波器。
語法：
b = fir2(n,f,m)
b = fir2(n,f,m,window)
b = fir2(n,f,m,npt)
b = fir2(n,f,m,npt,window)
b = fir2(n,f,m,npt,lap)
b = fir2(n,f,m,npt,lap,window)
參數n為濾波器的階數。
參數f為頻率點矢量，且f[0, 1], f=1對應於0.5fs。矢量f按升序排列，且第一個元素必須為0，最后一個必須為1，並可以包含重復的頻率點。
參數m為幅度點矢量，在矢量m中包含了與f相對應的期望得到的濾波器幅度。
參數Window用來指導所使用的窗函數類型，其默認值為漢明窗。
參數npt用來指定fir2函數對頻率響應進行內插的點數。
參數lap用來指定fir2函數在重復頻率點附近插入的區域大小。
36.3.2 fir2設計低通濾波器
fir2函數是用來設計任意頻率響應的各種加窗FIR濾波器，此函數使用也比較簡單，但是要采樣的頻率點和幅值不好把握，關於這個函數我們僅提供一個低通濾波器的設計。
原始信號是由50Hz正弦波和200Hz的正弦波組成，將200Hz的正弦波當做噪聲濾掉，下面通過函數fir2進行設計。其中頻率點矢量和幅度點矢量配置如下：
F = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1];
A = [1 1 1 1 0 0 0 0 0 0];
Matlab運行的代碼如下：
%**************************************************************************************** % fir2設計低通濾波器 %*************************************************************************************** fs=1000; %設置采樣頻率 1k N=1024; %采樣點數 n=0:N-1; t=0:1/fs:1-1/fs; %時間序列 f=n*fs/N; %頻率序列 Signal_Original=sin(2*pi*50*t); %信號50Hz正弦波 Signal_Noise=sin(2*pi*200*t); %噪聲200Hz正弦波 Mix_Signal=Signal_Original+Signal_Noise; %將信號Signal_Original和Signal_Original合成一個信號進行采樣 subplot(221); plot(t, Mix_Signal); %繪制信號Mix_Signal的波形 xlabel('時間'); ylabel('幅值'); title('原始信號'); grid on; subplot(222); y=fft(Mix_Signal, N); %對信號 Mix_Signal做FFT plot(f,abs(y)); xlabel('頻率/Hz'); ylabel('振幅'); title('原始信號FFT'); grid on; F = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1]; %表示要采樣的點 A = [1 1 1 1 0 0 0 0 0 0]; %表示采樣點的幅值 b = fir2(30, F, A); %30階FIR低通濾波器 %y2= filter(b, 1, x); y2=filtfilt(b,1,x); %經過FIR濾波器后得到的信號 Ps=sum(Signal_Original.^2); %信號的總功率 Pu=sum((y2-Signal_Original).^2); %剩余噪聲的功率 SNR=10*log10(Ps/Pu); %信噪比 y3=fft(y2, N); %經過FIR濾波器后得到的信號做FFT subplot(223); plot(f,abs(y3)); xlabel('頻率/Hz'); ylabel('振幅'); title('濾波后信號FFT'); grid on; [H,F]=freqz(b,1,512); %通過fir1設計的FIR系統的頻率響應 subplot(224); plot(F/pi,abs(H)); %繪制幅頻響應 xlabel('歸一化頻率'); title(['Order=',int2str(30),' SNR=',num2str(SNR)]); grid on;
Matlab運行結果如下：
從FFT結果和信噪比來看，fir2任意濾波器設計效果也是比較明顯的。
36.4 總結
本章節主要講解了函數fir1和函數fir2的使用，想深入的掌握這兩個函數，還需要大家多多練習。