加窗FIR濾波器設計實驗【內含實際使用filter的例子，討論了filter、fftfilt、filtfilt的差別】
轉自：http://blog.sina.com.cn/s/blog_6e8d34350100ng3f.html
題一：
利用加窗傅立葉級數法，設計一個具有如下指標的線性相位FIR低通濾波器：通帶截止頻率在4rad/s處，阻帶截止頻率在6rad/s處，最大通帶衰減為0.2dB，最小阻帶衰減為42dB，抽樣率為18rad/s。利用下面的各個窗函式進行設計：海明窗、漢寧窗和布萊克曼窗。對於每種情況，給出衝擊響應的係數並畫出設計的濾波器的增益響應。分析你的結果。不要使用M檔案fir1。
分析：由於題中給定的頻率為實際頻率，而數字濾波器一般使用歸一化頻率，故應該先對給定頻率資訊進行歸一化。
%程式設計：
close all
Omegap=4;
Omegas=6;
Omegat=18;
wp=2*pi*Omegap/Omegat;%0.44pi
ws=2*pi*Omegas/Omegat;%0.667pi %歸一化角頻率
alphap=0.2;
alphas=42;
NUM_Hamming=1;
NUM_Hann=2;
NUM_Blackman=3;
c=[3.32,3.11,5.56]*pi;
wc=(ws+wp)/2;
delt_w=ws-wp;
freq_labels={'用海明窗設計的FIR頻率響應','用漢寧窗設計的FIR頻率響應',...
'用布萊克曼窗設計的FIR頻率響應'};
ht_labels={'海明窗FIR衝激響應','漢寧窗FIR衝激響應','布萊克曼窗FIR響應'};
for filter_kind=NUM_Hamming:NUM_Blackman,
M=ceil(c(filter_kind)/delt_w);%向上取整
N=2*M+1;
switch filter_kind,
case NUM_Hamming,
win=hamming(N);
display(['海明窗生成的衝激響應係數：',' (階數N=',num2str(N),')']);
figure
case NUM_Hann,
win=hann(N);
display(['漢寧窗生成的衝激響應係數：',' (階數N=',num2str(N),')']);
figure
case NUM_Blackman,
win=blackman(N);
display(['布萊克曼窗生成的衝激響應係數：','(階數N=',num2str(N),')']);
figure
otherwise disp('error');
end
n=-M:M;
hd=sin(wc*n)./(pi*n);
hd(find(n==0))=wc*cos(wc*0)/pi;
%因為n=0時公式“sin(wc*n)./(pi*n)”中分母為0，計算結果為NaN,所以應該另用
%一種方法單獨計算n=0時的hd值。這裡採用洛必達法則對原hd公式進行分子分母微分
%得hd=wc*cos(w%c*n)./pi，n->0;這樣可以求得當n=0時原函式值。(也可以將n加上eps
%—-系統精度，使n最大程度接近0而不等0）
ht=hd.*win';
display(['',num2str(ht)]); %列印顯示衝激響應係數
%wvtool(ht);
subplot(1,2,1);
plot(n,ht,'.-') %畫衝激響應係數圖
title(ht_labels(filter_kind));
xlabel('n','FontSize',12);
ylabel('ht','fontsize',12);
grid on
[h,w]=freqz(ht,1,512);
W=w/pi;
H=20*log10(abs(h));
subplot(1,2,2); %畫頻響圖
hold on
title(freq_labels(filter_kind));
plot(W,H);
xlabel(['\pi (\omega/\omega','s)']);
ylabel('增益(dB)');
grid on
%為了更直觀的看到所設計FIR是否達到指標，可分別在圖中標出相關點
%從而可以對不同窗函式進行對比。標示過程如下：
dotp=round(mean(find(w>wp-0.1&w<wp+0.1)));%找到點wp的座標
dots=round(mean(find(w>ws-0.1&w<ws+0.1)));%找到點ws的座標
plot(W(dotp),H(dotp),'.r','MarkerSize',25);%對應處打上顯目的紅點
plot(W(dots),H(dots),'.r','MarkerSize',25);
text(W(dotp),H(dotp),['[','',num2str(round(W(dotp)*1000)/1000),...
',',num2str(round(H(dotp)*1000)/1000),'dB]'],'FontSize',15);%標示wp在頻響圖
%中的位置，從中可以%看來相應增益
text(W(dots),H(dots),['[','',num2str(round(W(dots)*100)/100),...
',',num2str(round(H(dots)*100)/100),'dB]'],'fontsize',15); %標示ws在頻響圖
%的位置
hold off
end
%產生的各種窗函式的係數：
0.0014702 -0.0013161 -0.001885 0.0038559 0.0022981 -0.0097177 3.3611e-017 0.019275 -0.0091462 -0.031341 0.031509 0.043366 -0.083816 -0.052269 0.31032 0.55556 0.31032 -0.052269 -0.083816 0.043366 0.031509 -0.031341 -0.0091462 0.019275 3.3611e-017 -0.0097177 0.0022981 0.0038559 -0.001885 -0.0013161 0.0014702
漢寧窗生成的衝激響應係數： (階數N=29)
-0 -0.0001973 0.0011375 0.0010796 -0.0059013 2.3913e-017 0.015232 -0.0077763 -0.028084 0.029338 0.041522 -0.081865 -0.051739 0.30954 0.55556 0.30954 -0.051739 -0.081865 0.041522 0.029338 -0.028084 -0.0077763 0.015232 2.3913e-017 -0.0059013 0.0010796 0.0011375 -0.0001973 -0
布萊克曼窗生成的衝激響應係數：(階數N=53)
-1.6732e-019 -5.7522e-006 -6.1629e-005 0.00011007 0.00021128 -0.00048407 -0.0003015 0.0013 -9.0902e-018 -0.002622 0.0012439 0.0042178 -0.0041222 -0.0053512 0.0092484 0.00464 -0.016847 5.1274e-017 0.026475 -0.011539 -0.036934 0.035186 0.046454 -0.087054 -0.053145 0.3116 0.55556 0.3116 -0.053145 -0.087054 0.046454 0.035186 -0.036934 -0.011539 0.026475 5.1274e-017 -0.016847 0.00464 0.0092484 -0.0053512 -0.0041222 0.0042178 0.0012439 -0.002622 -9.0902e-018 0.0013 -0.0003015 -0.00048407 0.00021128 0.00011007 -6.1629e-005 -5.7522e-006 -1.6732e-019
海明窗所得FIR的衝激響應圖及對應頻率響應圖：
其中，紅點標示處為wp：[0.445,-0.027dB] ;ws :[0.67,-50.79dB]
漢寧窗所得FIR的衝激響應及相應頻率響應：
Wp：[0.445,-0.067dB] ;ws :[0.64,-45.07dB]
布萊克曼窗所得FIR的衝激響應及相應頻率響應：
Wp :[0.445,-0.001dB] ;ws :[0.67,-78.18dB]
從以上三圖中可以看出，對於給定的指標[-0.2dB,42dB]，各窗函式所得FIR在[wp,ws]處的增益分別為：海明窗[-0.027dB,-50.79dB]
漢寧窗[-0.067dB,-45.07dB]
布萊克曼窗[-0.001dB,-78.18dB]
對比原指標，可以看出，漢寧窗比較接近指標，同時其階數（N=29）也是最少的。但就衰減效能而言，顯然海明窗的-50.79dB及布萊克曼窗的-78.18dB比漢寧窗的-45.07dB要好。從以上三圖還可以看出，海明窗阻帶波紋較大，而布萊克曼窗阻帶衰減較好但波紋較密，漢寧窗處於海明窗與布萊克曼窗的中間。
綜合上述特點，可以得出如下結論：對於加窗法FIR濾波器的設計，如果是要設計要求不高的濾波器，可以採用海明窗，因為海明窗阻帶衰減較好但阻帶波紋相對較大，而且濾波器階數相對不大。如果是要設計要求較高同時又不計成本及複雜性的FIR濾波器可以採用布萊克曼加窗，因為布萊克曼加窗法可以產生較好的阻帶衰減效能。雖然布萊克曼阻帶波紋較密，但隨頻率的增加波紋衰減越厲害，這對阻帶衰減有一定好處。不過，在階數方面，布萊克曼窗的階數又是最大的，這就會造成系統複雜性及成本的增加。綜合考慮阻帶效果及成本、複雜性等因素，加窗FIR設計最好採用漢寧窗，這樣可以在階數及阻帶衰減上有一個較好的取捨。
題二：
使用M檔案fir1，設計一個48階FIR帶通濾波器，通帶邊界的歸一化頻率為0.35和0.65。假設一個訊號，其中含f1=1Hz，f2=10Hz，f3=20Hz，三種頻率成分訊號的取樣頻率為50Hz。使用fftfilt函式濾波，將原訊號與通過濾波器的訊號進行比較。
%程式設計：
close all
w1=0.35;
w2=0.65;
wn=[w1,w2];
N=48;
b=fir1(N,wn,'bandpass');%得FIR係數
fvtool(b,1); %視覺化數字濾波器工具，可以同時得到頻響、相位、群延時、相位函式等圖形
f1=1;
f2=10;
f3=20;
fs=50;
% Omega1=2*pi*f1/fs;
% Omega2=2*pi*f2/fs;
% Omega3=2*pi*f3/fs;%歸一化角頻率
% t=-150:150;%頻率已歸一化，故間隔應為1.
% y1=cos(Omega1*t)+cos(Omega2*t)+cos(Omega3*t);
t=-1:1/fs:1;
y1=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); %這裡生成y1的方法也可以用上面註釋掉的部分的方法，實踐證明這兩種方法所得波形是一樣的（但前者時間軸必須進行一定的調整）
y2=fftfilt(b,y1); %濾波
subplot(2,1,1);
plot(t,y1,'.-');
title('原訊號');
xlabel('t');
ylabel('y1');
grid on
subplot(2,1,2);
plot(t,y2,'.-');
title('濾波後');
xlabel('t');
ylabel('y2');
grid on
%為更直觀地對濾波前後訊號進行分析，下面進行頻譜分析
figure
Len=length(y1);
NFFT=2^(nextpow2(Len));
Y1=fft(y1,NFFT); %快速傅立葉變換
Y2=fft(y2,NFFT);
f=((-NFFT/2:NFFT/2-1)/NFFT)*fs;%頻率軸調整使之與實際頻率對應
subplot(2,1,1)
semilogy(f,fftshift(abs(Y1))/NFFT,'.-r');
title('原訊號y1的頻譜','color','red','fontsize',15);
ylabel('頻譜幅度（log）');
xlabel('頻率（Hz）','color','red','fontsize',15);
set(gca,'xlim',[0,f(end)],'ylim',[0,max(abs(Y1)/NFFT)]);
grid on
subplot(2,1,2);
semilogy(f,fftshift(abs(Y2))/NFFT,'.-r');
title('濾波後的y2的頻譜','color','red','fontsize',15);
ylabel('頻譜幅度（log）');
xlabel('頻率（Hz）','color','red','fontsize',15);
set(gca,'xlim',[0,f(end)],'ylim',[0,max(abs(Y1)/NFFT)]);
grid on
所設計的FIR濾波器的頻率響應：
相位響應：（從下圖可以看出FIR濾波器優越的線性相位效能）
原訊號及濾波後的波形示於如下：
從上圖可以得到如下引申及結論：
（1）、按照題目要求，Wn分別為歸一化的：0.35(pi rad/s)、0.65(pi rad/s)，按照歸一化公式：
可知， 、
ω1=0.35pi，ω2=0.65pi，代入得：F1=8.75Hz、F2=16.25Hz。綜上，對於通帶為[ω1, ω2]的帶通濾波器，題中原訊號經過濾波後理論上被保留的應為10Hz的訊號。即上圖中“濾波後”理論上應該為10Hz。從圖中可以看出，事實也是如此的。
（2）、從上圖可以看出：“濾波後”的[-1,0.6]區間，波形出現了嚴重失真。原本應該為近似三角波的波形卻變為平坦地近乎直線。對於此現象剛開始很費解，故我們對已經接觸過的三個濾波函式（filter、filtfilt、fftfilt）進行了探索，以試圖找出三個濾波函式的本質區別及產生圖中失真的原因。首先，我們將本次實驗程式中的fftfilt分別用其他兩個函式代替，發現：對於filter函式，其結果與fftfilt的結果（無論是時域還是頻域上）是一樣的，而當使用filtfilt函式代替時，出現如下錯誤資訊：“Data must have length more than 3 times filter order.”意為使用filtfilt時，資料長度必須大於濾波器階數的3倍，改變長度又可以運行了.這說明程式中資料長度是有條件限制的.為更深層地找出它們的區別，我們通過找help文件及網上搜集的資料得到如下結論：
這三個函式的不同之處基本在於其實現演算法的差異。
(a)filter，名曰：“Generic filtering function”又稱一維數字濾波器（相對於filter2而言）、直接II型濾波器，也即為通用的數字濾波器。其演算法是通過如下公式“a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na)”進行濾波處理的。使用的是遞推演算法，引用網上的說法如下：
“ 用filter([1,1],[1,-0.5,0.06],x);來計算系統響應時，是假設系統初始狀態為全0
首先計算x(0)輸入系統的響應：y(0)=x(0)
然後 y(1)=0.5y(0)+x(1)+x(0)
接著 y(2)=0.5y(1)-0.6y(0)+x(2)+x(1)
以此類推，遞推的算出所有x（n）對應的系統輸出y(n)
這樣算出的y和x是等長度的”
這也就不難理解圖中的濾波波形了。本來對加窗FIR濾波器的係數ht就近似於sinc函式，其前面幾個係數（如海明窗生成的衝激響應係數： (階數N=31)
0.0014702 -0.0013161 -0.001885 0.0038559 0.0022981 -0.0097177 3.3611e-017 0.019275 -0.0091462 -0.031341…… ）都很小，甚至趨近於0，這樣在進行濾波時，由於x(n)要與係數相乘，自然所得結果剛開始會很小，由於使用遞推演算法，當濾波器處理的資料增加時，相加和的波動就越來越大，這才會漸漸出現正常的濾波波形。
(b)filtfilt，名曰：零相位數字濾波。其濾波演算法是基於filter而來的。只是filtfilt實現了零相位。其基本實現過程為先讓訊號用filter濾波，再將訊號時域反轉再次通過filter濾波，由《數字訊號處理》的知識可知，這樣兩次濾波後幅度響應為單個filter幅度的平方而相位實現了零相位。因為同樣使用filter的演算法，故在使用filtfilt濾波後，波形同樣會出現filter濾波的失真現象，即剛開始時輸出訊號近似為直線。
(c)fftfilt,名曰：基於FFT和重疊相加法的FIR濾波。英文描述為：“Filters with an FIR filter using the FFT” matlab描述為：“Overlap-add method of FIR filtering using FFT.”從網上資料來看，fftfilt似乎為專門的FIR濾波函式並一般應用於音訊訊號處理中。按照網上的說法，fftfilt只對FIR濾波器有效（從其輸入係數也可以看出這一點）。但對於具體的fftfilt實現公式無論是matlab自帶的help還是網上的介紹資料都比較少。故在此未能對fftfilt函式做更多討論。
綜合上述資料，對三個函式的總體感覺是：filter既能進行IIR濾波又能進行FIR濾波，而fftfilt雖然名曰“基於FFT和重疊相加法的FIR濾波”但總體感覺與filter相差無幾，只不過fftfilt只能對FIR濾波器有效。而對於filtfilt來說，它是一個使用了filter、比filter略高一個層次的濾波器，它可以實現零相位，這是它的優點，但對於普通濾波器來說，我們要知道的正是濾波器的幅度和相位特點，實現零相位會掩蓋濾波器的一些原有內涵，故filtfilt沒有filter更廣的適用性，它應該說是為特意達到零相位濾波而設計的。通過以上的學習，我們知道，已經接觸到的這三個濾波函式在功能各有千秋，雖沒有完全找到相關資料，但基本上都從原理上認識了三個濾波函式，這為以後的使用提供了很重要的參考。
做完了相關濾波函式的探討這個插曲後，現在讓我們把內容繼續迴歸到本次實驗的後續分析部分上來：
下圖是為了更準確地分析濾波前後訊號特點而專門增加的頻譜分析部分。圖中示出了濾波前後的頻譜分佈，兩個圖均以實際頻率為橫座標。
從第一個圖可以看出原訊號由於取樣頻率的影響存在許多噪聲紋波，但由於抽樣率仍符合奈奎斯特准則，故原訊號中的幾個頻率不會受太多影響，也就是說，雖然圖中有這麼多的噪聲紋波但訊號仍是能有效表達原訊號的。
在此我們重點分析第二個圖形（濾波後頻譜）。從第二個圖可以看出在0到10Hz，頻譜比較平緩而在10Hz到20Hz卻出現大量波紋。這點是否符合理論呢？首先，讓我們從頻響圖（本題第一個圖）上分析。從頻響圖可以看出，該帶濾波器大致以0.5pi rad/s（歸一化）為中心，左右成對稱特點。而下面第一個圖（原訊號）也大致以10Hz呈一定的對稱性，而10Hz恰好位於帶通濾波器的通帶內。由公式：Y(jw)=X(jw)H(jw)可知，濾波後頻譜實際上是頻率響應函式與原訊號頻譜函式的乘積，那麼這樣一來，從理論上來說，本題中訊號濾波後頻譜圖也應該有一定的對稱性。也就是說，在濾波後頻譜圖中，0到10Hz範圍內也應該像10到20Hz波形相似有大量的波紋，而實際卻是平緩的，這是為什麼呢？
這裡首先引用上面剛才對filter、fftfilt等三個函式進行分析的相關內容。從上面的分析可知使用這些函式進行濾波時，由於演算法上的原因，濾波後訊號在最前面的部分訊號是平坦的、趨於直線的失真訊號，那麼當對該段訊號進行頻譜分析時，自然得到的頻率就是比較穩定的“低頻”成份了，頻譜圖中的低頻部分就主要是這些因素造成的。
另外這裡還有一種理解：由於題中濾波設計的阻帶位於[0.35,0.65]之間，轉換為以50Hz為取樣率的實際頻率就是：[8.75,16.25]Hz之間，可見，對於原訊號中包含的10Hz頻率，其分佈靠近8.75Hz截止頻率，也就是說，10Hz以下的頻率衰減很大，而10到16.25Hz的噪聲基本沒被衰減，這也就造成濾波後10Hz以下頻譜平滑而10Hz以上波紋較多這一“不對稱”現象。個人覺得這種理解更有道理些。
當然，從本圖也可以更清楚地看到原訊號1Hz與20Hz頻率成份基本被濾除而10Hz頻率成分被保留的這一結果。
作者：華南理工大學08電信微電2班---飛鳥