數(shù)字濾波器的設(shè)計(jì)實(shí)驗(yàn)

1、實(shí)驗(yàn)二 IIR數(shù)字濾波器的設(shè)計(jì)實(shí)驗(yàn)內(nèi)容及步驟:數(shù)字濾波器的性能指標(biāo)：通帶臨界頻率fp、阻帶臨界頻率fr；通帶內(nèi)的最大衰減Ap；阻帶內(nèi)的最小衰減Ar；采樣周期T；(1)、fp=0.3KHz,Ap=0.8dB, fr=0.2KHz,Ar=20dB,T=1ms；設(shè)計(jì)一Chebyshev高通濾波器；觀察其通帶損耗和阻帶衰減是否滿足要求。程序如下：fp=300; fr=200; Ap=0.8; Ar=20; T=0.001;fs=1/T; wp=2*pi*fp*T;wr=2*pi*fr*T; Wp=2/T*tan(wp/2);Wr=2/T*tan(wr/2);N,Wn=cheb1ord(Wp,Wr,Ap

2、,Ar,'s');B,A = cheby1(N,Ap,Wn,'high','s'); num,den=bilinear(B,A,1/T); h,w=freqz(num,den);plot(w*fs/(2*pi),20*log10(abs(h); %衰減及頻率都用歸一化的1為單位顯示axis(0,500,-30,0);title('Chebyshev高通濾波器');xlabel('頻率');ylabel('衰減');grid on;根據(jù)下圖知道通帶損耗與阻帶衰減滿足要求(2)、fp=0.2KHz,Ap

3、=1dB, fr=0.3KHz,Ar=25dB,T=1ms；分別用脈沖響應(yīng)不變法及雙線性變換法設(shè)計(jì)一Butterworth數(shù)字低通濾波器，觀察所設(shè)計(jì)數(shù)字濾波器的幅頻特性曲線，記錄帶寬和衰減量，檢查是否滿足要求。比較這兩種方法的優(yōu)缺點(diǎn)。程序如下：fp=200; fr=300; Ap=1;Ar=25; T=0.001;fs=1/T; wp=2*pi*fp*T;wr=2*pi*fr*T; Wp=2/T*tan(wp/2);Wr=2/T*tan(wr/2);N,Wn=buttord(Wp,Wr,Ap,Ar,'s');B,A = butter(N,Wn,'s'); num

4、1,den1=impinvar(B,A,1/T); %脈沖響應(yīng)不變法得出設(shè)計(jì)的傳遞函數(shù) num2,den2=bilinear(B,A,1/T); %雙線性變換法得出設(shè)計(jì)的傳遞函數(shù)h1,w=freqz(num1,den1); plot(w*fs/(2*pi),20*log10(abs(h2),w*fs/(2*pi),20*log10(abs(h1),'r.');grid on; %衰減及頻率都用歸一化的1為單位顯示axis(0,500,-30,0);title('Butterworth低通濾波器（紅線脈沖響應(yīng)不變法 藍(lán)線雙線性變換法）');xlabel('

5、ƵÂÊ');ylabel('Ë¥¼õ');grid on; 優(yōu)缺點(diǎn)：采用脈沖響應(yīng)不變法優(yōu)點(diǎn)： 1.h(n)完全模仿模擬濾波器的單位抽樣響應(yīng)時(shí)域逼近良好 2線性相位模擬濾波器轉(zhuǎn)變?yōu)榫€性相位數(shù)字濾波器 缺點(diǎn)： 1.對(duì)時(shí)域的采樣會(huì)造成頻域的“混疊效應(yīng)”，故有可能使所設(shè)計(jì)數(shù)字濾波器的頻率響應(yīng)與原來模擬濾波器的頻率響應(yīng)相差很大2不能用來設(shè)計(jì)高通和帶阻濾波器。只適用于限帶的低通、帶通濾波器采用雙線性變換法優(yōu)點(diǎn)：1避免了頻率響應(yīng)的混迭失真現(xiàn)象 2在特定數(shù)字濾波器和特定模擬濾波器處，頻率響應(yīng)是嚴(yán)格相等的，

6、它可以較準(zhǔn)確地控制截止頻率的位置。 3它是一種簡(jiǎn)單的代數(shù)關(guān)系，設(shè)計(jì)十分方便。 缺點(diǎn)： 1除了零頻率附近，w與W之間嚴(yán)重非線性，即線性相位模擬濾波器變?yōu)榉蔷€性相位數(shù)字濾波器 2要求模擬濾波器的幅頻響應(yīng)為分段常數(shù)型，不然會(huì)產(chǎn)生畸變 3對(duì)于分段常數(shù)型AF濾波器，經(jīng)雙線性變換后，仍得到幅頻特性為分段常數(shù)的 DF.但在各個(gè)分段邊緣的臨界頻率點(diǎn)產(chǎn)生畸變，這種頻率的畸變，可通過頻率預(yù)畸變加以校正(3)、利用雙線性變換法分別設(shè)計(jì)滿足下列指標(biāo)的Butterworth型和Chebyshev型數(shù)字低通濾波器，并作圖驗(yàn)證設(shè)計(jì)結(jié)果。fp=1.2kHz, Ap0.5dB, fr=2KHz, Ar40dB, fs=8KHz

7、程序如下：fp=1200;fr=2000; Ap=0.5; Ar=40; fs=8000; T=1/fs; wp=2*pi*fp*T; wr=2*pi*fr*T; Wp=2/T*tan(abs(wp/2); Wr=2/T*tan(abs(wr/2);N,Wn=buttord(Wp,Wr,Ap,Ar,'s');B,A = butter(N,Wn,'s'); num2,den2=bilinear(B,A,1/T); %雙線性變換法得到設(shè)計(jì)的傳遞函數(shù)h2,w=freqz(num2,den2); figure(1);subplot(2,1,1)plot(w/pi*fs/

8、2,10*log10(abs(h2).2); axis(0 fs/2 -50 0); title('Butterworth低通濾波器');xlabel('頻率 Hz'); ylabel('衰減 dB'); grid on; subplot(2,1,2)N,Wn=cheb1ord(Wp,Wr,Ap,Ar,'s');B,A = cheby1(N,Ap,Wn,'s'); num,den=bilinear(B,A,1/T); h1,w=freqz(num,den);plot(w/pi*fs/2,10*log10(abs(h

9、1).2); axis(0 fs/2 -50 0); title('chebyshev低通濾波器');xlabel('頻率 Hz'); ylabel('衰減 dB');grid on;    (4)、利用雙線性變換法設(shè)計(jì)一Butterworth型數(shù)字帶通濾波器，已知fs=30KHz，其等效的模擬濾波器指標(biāo)為Ap3dB, 2KHzf3KHz, Ar5dB, f6KHz, Ar20dB, f1.5KHz  程序如下：f1=2000; f3=3000;fsl=1500; fsh=6000;rp=3; rs=20;

10、Fs=30000;wp1=2*pi*f1/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=wp1 wp3;ws=wsl wsh;wap=2*Fs*tan(wp./2);was=2*Fs*tan(ws./2);n,wn=buttord(wap,was,rp,rs,'s');z,p,k=buttap(n);bp,ap=zp2tf(z,p,k);bw=wap(2)-wap(1);w0=sqrt(wap(2)*was(1);bs,as=lp2bp(bp,ap,w0,bw);bz1,az1=bilinear(bs,as,Fs);h

11、,w=freqz(bz1,az1,256,Fs);plot(w,20*log10(abs(h);axis(0 8000 -30 0); title('Butterworth帶通濾波器');xlabel('頻率 Hz');ylabel('衰減 dB');grid on; 1 雙線性變換法中和之間的關(guān)系是非線性的，在實(shí)驗(yàn)中你注意到這種非線性關(guān)系了嗎？從哪幾種數(shù)字濾波器的幅頻特性曲線中可以觀察到這種非線性關(guān)系？   答：在雙線性變換法中，模擬頻率與數(shù)字頻率不再是線性關(guān)系，所以一個(gè)線性相位模擬器經(jīng)過雙線性變換后得到的數(shù)字濾波

12、器不再保持原有的線性相位了。如以上實(shí)驗(yàn)過程中，采用雙線性變化法設(shè)計(jì)的butter和cheby1數(shù)字濾波器，從圖中可以看到這種非線性關(guān)系。  2能否利用公式完成脈沖響應(yīng)不變法的數(shù)字濾波器設(shè)計(jì)？為什么？  答：不能，這樣會(huì)使得H(z)中的z以對(duì)數(shù)形式出現(xiàn)，使H(z)的分子分母不再是z的有理多項(xiàng)式，這樣會(huì)給系統(tǒng)的分析和實(shí)現(xiàn)帶來很大的困難。實(shí)驗(yàn)三、FIR數(shù)字濾波器的設(shè)計(jì)實(shí)驗(yàn)內(nèi)容及步驟(1) N=15, 。用Hanning窗設(shè)計(jì)一線性相位帶通濾波器，觀察它的實(shí)際3dB和20dB帶寬。N=45，重復(fù)這一設(shè)計(jì)，觀察幅頻和相位特性的變化，注意長(zhǎng)度N變化的影響；程序如下：

13、N1=15;w1=0.3;w2=0.5;wn=w1,w2; b1=fir1(N1-1,wn,hanning(N1); %用Hanning窗作為沖擊響應(yīng)的窗函數(shù)figure(1);freqz(b1,1);title('N=15 hanning窗幅頻相頻特性'); N2=45;b2=fir1(N2-1,wn,hanning(N2);figure(2);freqz(b2,1);title('N=45 hanning窗幅頻相頻特性');由圖可知：隨著N值的增大，取樣值增大，主瓣的寬度減小，但是幅頻特性與相頻特性曲線的波動(dòng)頻率會(huì)增加。(2) 分別改用矩形窗和Blackma

14、n窗，設(shè)計(jì)(1)中的帶通濾波器，觀察并記錄窗函數(shù)對(duì)濾波器幅頻特性的影響，比較三種窗的特點(diǎn)；程序如下：N1=15;w1=0.3; w2=0.5; wn=w1,w2; b1=fir1(N1-1,wn,boxcar(N1); %用boxcar作為沖擊響應(yīng)的窗函數(shù)figure(1);freqz(b1,1);title('N=15 boxcar窗幅頻相頻特性'); N2=45;b2=fir1(N2-1,wn,boxcar(N2);figure(2);freqz(b2,1);title('N=45 boxca窗幅頻相頻特性'); b3=fir1(N1-1,wn,blackm

15、an(N1); %用blackman作為沖擊響應(yīng)的窗函數(shù)figure(3);freqz(b3,1);title('N=15 blackman窗幅頻相頻特性'); N2=45;b4=fir1(N2-1,wn,blackman(N2);figure(4);freqz(b4,1);title('N=45 blackman窗幅頻相頻特性');分析：由圖可知，在三個(gè)窗函數(shù)中，Blackman窗的衰減性最好，矩形窗的衰減最差。，而矩形窗的過渡帶帶寬最小，Blackman窗的過渡帶帶寬最大。Hanning窗的性能比較平均。其實(shí)，改善阻帶衰減的一種方法是加寬過渡帶寬，以犧牲過渡

16、帶寬換取阻帶衰減的增大。也就是以增加主瓣寬度為代價(jià)來降低旁瓣。    (3)用Kaiser窗設(shè)計(jì)一專用線性相位濾波器，N=40， 如圖，當(dāng)0=4，6，10時(shí)，分別設(shè)計(jì)，比較它們的幅頻和相頻特性，注意0取不同值時(shí)的影響；程序如下：N=40;beta=4,6,10; wn=0.2,0.4,0.6,0.8; s='rgb' for i=1:3 window=kaiser(N,beta(i); b1=fir1(N-1,wn,window); H,w=freqz(b1,1); magH=abs(H);angH=angle(H); subplot(2,1,1)

17、 plot(w/pi, magH,s(i);hold on; grid on xlim(0,1); ylim(0,1); xlabel('w/（2*pi）');ylabel('|H( ejw)|') title('N=40 幅度部分'); legend('=4','=6','=10'); subplot(2,1,2); plot(w/pi,angH,s(i);hold on; grid on; xlim(0,1); xlabel('w/(2*pi)');ylabel('(w)

18、'); title('N=40 相頻部分' );end由下圖可知，隨著的增大，幅頻曲線變得越平緩，但其邊峰瓣值的衰減越大，相頻特性變得越好。(4)、用頻率采樣法設(shè)計(jì)(3)中的濾波器，過渡帶分別設(shè)一個(gè)過渡點(diǎn)，令H(k)=0.5。比較兩種不同方法的結(jié)果；分析：，。而在通帶兩側(cè)的過渡帶需要插入0.5。采用第一種相位特性，因此，可以得到以下式子程序如下：Hr=zeros(1,3),0.5,ones(1,5),0.5,0,0.5,ones(1,5),0.5,zeros(1,5),-0.5,-ones(1,5),-0.5,0,-0.5,-ones(1,5),-0.5,0,0;N =

19、 40;k = 0:N-1; subplot(2,2,1);stem(k,Hr,'.');grid on;xlabel('n'),ylabel('h(n)'),title('插值圖'); wm = 2*pi*k./N;Hd = Hr.*exp(-j*(N-1)*wm*0.5);hh = real(ifft(Hd);w1 = linspace(0,pi-0.1,1000);H1 = freqz(hh,1,w1);amgH = abs(H1);dbH = 20*log10(amgH);angH = angle(H1); subplot

20、(2,2,2);plot(w1/pi,dbH),grid on;title('幅度響應(yīng)');xlabel('w/pi'),ylabel('|H( ejw)|/dB')axis(0 1 -100 20);subplot(2,2,3);plot(w1/pi,angH),grid on;title('相頻特性');xlabel('w/pi'),ylabel('')subplot(2,2,4);plot(w1/pi,amgH),grid on;title('幅頻特性');xlabel(

21、9;w/pi'),ylabel('|H( ejw)|')（5）、用雷米茲(Remez)交替算法設(shè)計(jì)(3)中的濾波器，并比較(3)、(4)、(5)三種方法的結(jié)果。 程序如下：f=0 0.38 0.42 0.78 0.82 1;A=0 0 1 1 0 0;weigh=1 8 1; b=remez(40,f,A,weigh);%h,w=freqz(b,1,512,'whole');magh=abs(h);dbh=20*log10(magh);phah=angle(h); subplot(2,1,1)plot(w/(2*pi),dbh);grid on

22、;xlabel('w/(2*pi)'),ylabel('|H( ejw)|/dB'),title('幅頻特性');axis(0 1 -70 5) subplot(2,1,2)plot(w/(2*pi),phah);grid on;xlabel('w/（2*pi）'),ylabel('(w)'),title('相頻特性'); 四、實(shí)驗(yàn)思考1. 定性地說明用本實(shí)驗(yàn)程序設(shè)計(jì)的FIR濾波器的3dB截止頻率在什么位置？它等于理想頻率響應(yīng)Hd(ej)的截止頻率嗎？ 答：在之間，大致等于理想頻率響應(yīng)Hd

23、(ej)的截止頻率。2. 如果沒有給定h(n)的長(zhǎng)度N，而是給定了通帶邊緣截止頻率c和阻帶臨界頻率p，以及相應(yīng)的衰減，你能根據(jù)這些條件用窗函數(shù)法設(shè)計(jì)線性相位FIR低通濾波器嗎？答：能，可以根據(jù)通帶邊緣截止頻率c和阻帶臨界頻率p，以及相應(yīng)的衰減與所選用的窗函數(shù)來大致判斷其長(zhǎng)度N，從而設(shè)計(jì)線性相位FIR低通濾波器。實(shí)驗(yàn)四、用MATLAB信號(hào)處理工具箱進(jìn)行濾波器設(shè)計(jì)訓(xùn)練針對(duì)一個(gè)含有5Hz、50Hz和100Hz的混和正弦波信號(hào)，設(shè)計(jì)一個(gè)FIR帶通濾波器并對(duì)其進(jìn)行濾波。利用MATLAB實(shí)現(xiàn)的三種方法：程序設(shè)計(jì)法、 FDATool設(shè)計(jì)法和SPTool設(shè)計(jì)法。參數(shù)要求：采樣頻率fs=100Hz，通帶下限截止

24、頻率fc1=10 Hz，通帶上限截止頻率fc2=20 Hz，過渡帶寬6 Hz，通阻帶波動(dòng)0.01，采用凱塞窗設(shè)計(jì)。方法1 ：程序設(shè)計(jì)法fsamp = 100;fcuts = 4 10 20 26;mags = 0 1 0;devs = 0.01 0.01 0.01;n,wn,beta,ftype = kaiserord(fcuts,mags,devs,fsamp);n = n + rem(n,2);hh = fir1(n,wn,ftype,kaiser(n+1,beta),'noscale');H,f = freqz(hh,1,1024,fsamp); figure(1)plo

25、t(f,abs(H),grid on;xlabel('f/Hz'),ylabel('|H( ejw)|'),title('kaiser數(shù)字帶通濾波器');用一個(gè)序列測(cè)試濾波器性能figure(2)N = 256;n = 0:1:N-1;t =n/fsamp;xn = sin(2*pi*5*t)+sin(2*pi*15*t)+sin(2*pi*35*t); subplot(2,2,1)stem(n,xn,'.');grid on;xlabel('n');ylabel('x(n)');title(

26、9;序列x(n)');axis(0 256 -2 2); X = fft(xn,N);magX = abs(X);f = n*fsamp/N; subplot(2,2,2);plot(f(1:N/2),magX(1:N/2);grid on;%axis(0 fs 0 1);xlabel('f/Hz');ylabel('|X(ejw)|');title('序列x(n)幅頻特性'); x = xn;y = filter(hh,1,x);Y = fft(y,N);magY = abs(Y); subplot(2,2,3);stem(n,y,&#

27、39;.');grid on;xlabel('n');ylabel('y(n)');title('濾波后的序列y(n)');axis(0 256 -2 2); subplot(2,2,4);plot(f(1:N/2),magY(1:N/2);grid on;xlabel('f/Hz');ylabel('|Y(ejw)|');title('序列y(n)幅頻特性');方法2：利用FDATool設(shè)計(jì)法在MATLAB命令窗口輸入FDATool后回車就會(huì)彈出FDATool界面。  帶通濾波器設(shè)計(jì):   已知濾波器的階數(shù)n=38，beta=3.4。(1) 首先在Filter Type中選擇Bandpass；(2) 在Design Method選項(xiàng)中選擇FIR Window，(3) 接著在Window選項(xiàng)中選取Kaiser，Beta值為3.4；(4) 指定Filter Order項(xiàng)中的Specify order為38；采樣頻率Fs=100Hz,截止頻率Fc1=10Hz,Fc2=20Hz。(5) 設(shè)置完以后點(diǎn)擊窗口下方的Design Filter，在窗口上方就會(huì)看