1、實驗4 基于MATLAB的FIR數字濾波器設計實驗目的：加深對數字濾波器的常用指標和設計過程的理解。實驗原理：低通濾波器的常用指標：通帶邊緣頻率，阻帶邊緣頻率 ，通帶起伏，通帶峰值起伏，阻帶起伏，最小阻帶衰減。數字濾波器有IIR和FIR兩種類型，它們的特點和設計方法不同。 在MATLAB中，可以用b=fir1(N,Wn,ftype,taper) 等函數輔助設計FIR數字濾波器。N代表濾波器階數；Wn代表濾波器的截止頻率(歸一化頻率)，當設計帶通和帶阻濾波器時，Wn為雙元素相量；ftype代表濾波器類型，如high高通，stop帶阻等；taper為窗函數類型，默認為海明窗，窗系數需要實現用窗函數

2、blackman, hamming,hanning chebwin, kaiser產生。例1 用凱塞窗設計一FIR低通濾波器，通帶邊界頻率 ，阻帶邊界頻率 ，阻帶衰減 不小于50dB。 解 首先由過渡帶寬和阻帶衰減 來決定凱塞窗的N和 ， ，  上圖給出了以上設計的頻率特性，(a) 為N=30直接截取的頻率特性(b)為凱塞窗設計的頻率特性。凱塞窗設計對應的MATLAB程序為： wn=kaiser(30,4.55); nn=0:1:29; alfa=(30-1)/2; hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa); h=hd.*wn; h1,w1=fre

3、qz(h,1); 或者：b = fir1(29,0.4,kaiser(30,4.55); h1,w1=freqz(b,1); plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10); grid;xlabel('歸一化頻率/p') ;ylabel('幅度/dB') ;還可以使用n,Wn,beta,ftype = kaiserord(f,a,dev)函數來估計濾波器階數等，得到凱塞窗濾波器：fcuts = 0.3 0.5; %歸一化頻率omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,bet

4、a,ftype = kaiserord(fcuts,mags,devs); %計算出凱塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); freqz(hh);   實際中，一般調用MATLAB信號處理工具箱函數remezord來計算等波紋濾波器階數N和加權函數W()，調用函數remez可進行等波紋濾波器的設計，直接求出濾波器系數。函數remezord中的數組fedge為通帶和阻帶邊界頻率，數組mval是兩個邊界處的幅值，而數組dev是通帶和阻帶的波動，fs是采樣頻率單位為Hz。例2

5、 利用雷米茲交替算法設計等波紋濾波器，設計一個線性相位低通FIR數字濾波器，其指標為：通帶邊界頻率fc=800Hz，阻帶邊界fr=1000Hz，通帶波動 阻帶最小衰減At=40dB，采樣頻率fs=4000Hz。 解 在MATLAB中可以用remezord 和remez兩個函數設計，其結果如圖2，MATLAB程序如下： fedge=800 1000; mval=1 0; dev=0.0559 0.01;                 fs=400

6、0; N,fpts,mag,wt=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); h,w=freqz(b,1,256); plot(w*2000/pi,20*log10(abs(h); grid;                    xlabel('頻率/Hz') ; ylabel('幅度/dB'); &#

7、160;                   一、實驗內容： 利用MATLAB編程設計一個數字帶通濾波器，指標要求如下：通帶邊緣頻率：，通帶峰值起伏：。阻帶邊緣頻率：，最小阻帶衰減： 。 分別用窗函數法和等波紋濾波器法設計兩種FIR數字濾波器。實驗要求：給出FIR數字濾波器的沖激響應，繪出它們的幅度和相位頻響曲線，討論它們各自的實現形式和特點。1-1）用窗函數法實現：調用函數n,wn,bta,ftype=kaiserord（f

8、，a，dev，fs）參數：f=0.3 0.45 0.65 0.8為對應數字頻率， a=0 1 0為由f指定的各個頻帶上的幅值向量，一般只有0和1表示；和f長度關系為（2*a的長度）2=（f的長度）devs=0.01 0.1087 0.01用于指定各個頻帶輸出濾波器的頻率響應與其期望幅值之間的最大輸出誤差或偏差，長度與a相等，計算公式：阻帶衰減誤差=，通帶波動衰減誤差=fs缺省值為2HZ>> n,wn,bta,ftype=kaiserord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用kaiserord函數估計出濾波器階數n和beta參數&

9、gt;> h1=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale');>> hh1,w1=freqz(h1,1,256);>> figure(1)>> subplot(2,1,1)>> plot(w1/pi,20*log10(abs(hh1)>> grid>> xlabel('歸一化頻率w');ylabel('幅度/db');>> subplot(2,1,2)>> plot(w1/pi,angle(hh1)>

10、> grid>> xlabel('歸一化頻率w');ylabel('相位/rad');h1 = Columns 1 through 8 0.0041 0.0055 -0.0091 -0.0018 -0.0056 -0.0000 0.0391 -0.0152 Columns 9 through 16 -0.0381 0.0077 -0.0293 0.0940 0.0907 -0.2630 -0.0517 0.3500 Columns 17 through 24 -0.0517 -0.2630 0.0907 0.0940 -0.0293 0.007

11、7 -0.0381 -0.0152 Columns 25 through 31 0.0391 -0.0000 -0.0056 -0.0018 -0.0091 0.0055 0.0041 圖4-1如果直接用freqz(h1,1,256);畫圖得： 1-2）用等波紋法設計：調用函數n,fpts,mag,wt=remezord(f，a，dev)f=0.3 0.45 0.65 0.8a=0 1 0dev=0.01 0.1087 0.01其含義同函數n,wn,bta,ftype=kaiserord（f，a，dev，fs）中的參數相同。 >> n,fpts,mag,wt=remezord(0.

12、3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用remezord函數估算出remez函數要用到的階n、歸一化頻帶邊緣矢量fpts、頻帶內幅值響應矢量mag及加權矢量w，使remez函數設計出的濾波器滿足f、a及dev指定的性能要求。>> h2=remez(n,fpts,mag,wt);%設計出等波紋濾波器>> hh2,w2=freqz(h2,1,256);>> figure(2)>> subplot(2,1,1)>> plot(w2/pi,20*log10(abs(hh2)>> grid

13、>> xlabel('歸一化頻率w');ylabel('幅度/db');>> subplot(2,1,2)>> plot(w2/pi,angle(hh2)>> grid>> xlabel('歸一化頻率w');ylabel('相位/rad');>> h2h2 = Columns 1 through 9 -0.0013 0.0092 -0.0255 -0.0642 0.1177 0.0922 -0.2466 -0.0466 0.3116 Columns 10 th

14、rough 17 -0.0466 -0.2466 0.0922 0.1177 -0.0642 -0.0255 0.0092 -0.0013 圖4-2用freqz(h2,1,256);直接得圖：二、對課本作業9.23畫圖。2-1）用漢寧窗實現：1）手動計算>> n=0.001:58.001;>> hd=sin(0.18125*pi*(n-29)./(pi*(n-29);>> win=0.5+0.5*cos(2*pi*(n-29)/58);>> h1=2*cos(pi*(n-29)/2).*hd.*win;>> hh1,w1=freqz(

15、h1,1,256);>> figure(1)>> subplot(2,1,1)>> plot(w1,20*log10(abs(hh1)>> grid>> xlabel('數字頻率w/rad');ylabel('幅度/db');>> subplot(2,1,2)>> plot(w1,angle(hh1)>> grid>> xlabel('數字頻率w/rad');ylabel('相位/rad'); 圖4-32-2）用自帶的fir

16、1函數：>> n=59;>> wn=3/8 5/8;>> h2=fir1(n,wn,'bandpass',hann(n+1);>> hh2,w2=freqz(h2,1,256);>> figure(2)>> subplot(2,1,1)>> plot(w2,20*log10(abs(hh2)>> xlabel('數字頻率w/rad');ylabel('幅度/db');>> grid>> subplot(2,1,2)>>

17、; plot(w2,angle(hh2)>> xlabel('數字頻率w/rad');ylabel('相角/rad');>> grid 圖4-42-3）用等波紋法設計：調用函數n,fpts,mag,wt=remezord(f，a，dev)f=0.2625 0.375 0.625 0.7375a=0 1 0dev=0.01 0.1087 0.01 >> n,fpts,mag,wt=remezord(0.2625 0.375 0.625 0.7375,0 1 0,0.01 0.1087 0.01); >> h2=rem

18、ez(n,fpts,mag,wt);>> hh2,w2=freqz(h2,1,256);>> figure(2)>> subplot(2,1,1)>> plot(w2/pi,20*log10(abs(hh2)>> grid>> xlabel('歸一化頻率w');ylabel('幅度/db');>> subplot(2,1,2)>> plot(w2/pi,angle(hh2)>> grid>> xlabel('歸一化頻率w');y

19、label('相位/rad');>> h2畫圖為：另外帶通濾波器還可以用低通和高通級聯的方法實現，因為步驟過多，這里不作討論。三、實驗分析：1）總結：FIR濾波器實現一般采用窗函數法和等紋波設計法。窗函數法還包含兩個分支，一種是用公式先手動算出N值和其他對應得窗函數參數值，再代入窗函數和fir1實現，一種是用函數*rord估算出N和相應參數再用fir1實現。不過要注意*rord會低估或高估階次n，可能會使濾波器達不到指定的性能，這時應稍微增加或降低階次。如果截止頻率在0或Nyquist頻率附近，或者設定的dev值較大，則得不到正確結果。2）濾波器實現形式及特點： 由

20、于一般的濾波器在利用窗函數是其通帶波紋和阻帶波紋不同（一般為第一個阻帶波紋最大）因此，在滿足第一個阻帶衰減旁瓣時，比其頻率高的旁瓣，它們的衰減都大大超出要求。而根據阻帶衰減與項數的近似關系可得當阻帶衰減越大，所需項數越多。2）等波紋設計和窗函數設計不同之處在于實現形式及特點： 窗函數設計是通過最小平方積分辦法來設計的即該濾波器的誤差為： 即要求最小方法來設計濾波器，這樣的濾波器更忠實于理想濾波器（即濾波系數更接近于理想濾波器。證明如下： = = = = 因此，幅度頻譜差值越小，實際濾波器就越接近理想濾波器。 而等波紋濾波器是通過最大加權誤差最小化來實現，其誤差為： 要求該誤差最小來實現濾波器，

21、得出來的濾波系數較窗函數設計相差較遠。以下通過對題目中的h1及h2作比較。% sigsum是用來對數組各元素進行求和function y=sigsum(n1,n2,n,x);y=0; for i=n1+1-min(n):n2+1-min(n) y=y+x(i);end>> n=0.001:30.001;>> h=2*cos(0.55*pi*(n-15).*sin(0.175*pi*(n-15)./(pi*(n-15);>> delta1=h-h1;>> n=0.001:16.001;>> h=2*cos(0.55*pi*(n-15).*sin(0.175*pi*(n-15)./(pi*(n-15);>> delta2=h-h2;>> y1=sigsum(0,30,0:30,(abs(delta1).2)/31;>> y2=sigsum(0,16,0:16,(abs(delta2).2)/17;>> y1y1 = 1.9099e-004>> y2y2 = 0.0278由此得到用窗函數實現的濾波系數比用等波紋濾波器系數的每一項更接近于理想濾波器（y1為用窗函數