數字信號處理課程設計實驗報告(20211126002056)

1、華 北 電 力 大 學實 驗 報 告第 1 頁實驗環(huán)境matlab 6.5 實驗名稱實驗一： fft 的應用實驗目的1、熟悉 matlab 在數字信號處理中的應用。2、掌握利用 fft 計算序列線性卷積的基本原理及編程實現。3、掌握對連續(xù)信號進行采樣的基本原理和方法，并利用fft 對信號進行頻譜分析。實驗內容1.對于兩個序列： x(n)=nr16(n),h(n)=r8(n) （1）在同一圖形窗口中繪出兩序列的時域圖形。（2）利用 fft 編程計算兩序列的線性卷積，繪出的時域圖形。2.利用 fft 對信號進行譜分析對于連續(xù)信號xa(t)=cos(2f1t) +5cos(2f2t) +cos(2

2、f3t) ，其中f1=6.5khz, f2=7khz, f3=9khz, 以采樣頻率 fs=32 khz 對其進行采樣，（1）對 xa(t) 信號采集 16 點樣本，分別作16 點和補零到 256 點的 fft，并分別繪出對應的幅頻特性曲線。（2）對 xa(t)信號采集 256 點樣本，分別作 256 點和 512 點的 fft，并分別繪出對應的幅頻特性曲線。（3）比較（ 1）和（ 2）中的結果，分析采樣點數和傅里葉變換點數對fft的影響，說明高密度頻譜和高分辨率頻譜的特點與區(qū)別。題目分析：1.利用 fft 計算線性卷積首先.周期卷積是線性卷積以l 為周期的周期延拓序列的主值序列； 兩個長度為

3、 m，n 的序列的線性卷積可用長度均為l 的循環(huán)卷積來代替，用循環(huán)卷積計算線性卷積的方法歸納如下: 將長為 n2的序列 x(n) 延長到 l, 補 l-n2 個零將長為 n1的序列 h(n) 延長到 l, 補 l-n1 個零如果 ln1+n2-1,則循環(huán)卷積與線性卷積相等, 此時, 可有 fft 計算線性卷積, 方法如下 : a.計算 x(k)=fftx(n) b.求 h(k)=ffth(n) c.求 y(k)=h(k)y(k) k=0l-1 d.求 y(n)=iffty(k) n=0l-1 可見，只要進行二次fft,一次 ifft 就可完成線性卷積計算。計算表明,l32 時, 上述計算線性卷

4、積的方法比直接計算線卷積有明顯的優(yōu)越性, 因此, 也稱上述循環(huán)卷積方法為快速卷積法。題目利用 fft 計算線性卷積以下為程序及結果：華 北 電 力 大 學實 驗 報 告第 2 頁實驗結果及分析n=1:15; x=n*ones(1,n); m=8; h=ones(1,8); subplot(3,1,1) stem(n,x);title(x(n)=n*r16(n); subplot(3,1,2) stem(0:7,h);title(h(n)=r8(n); %stem表示離散序列x2=0:15 zeros(1,7) %補零操作h2=ones(1,8) zeros(1,15) x3=fft(x2) %

5、計算 x2 的快速傅氏變換x3 h3=fft(h2) y2=x3.*h3 y=ifft(y2) subplot(3,1,3) stem(0:22,y);title( 卷積結果 ); ylabel(對數幅度、 db);xlabel( 以pi 為單位的頻率 ); 題目分析：利用 fft 進行頻譜觀測分析頻率分辨率可以理解為在使用dft 時，在頻率軸上的所能得到的最小頻率間隔 f0=fs/n=1/nts=1/t, 其中 n 為采樣點數， fs 為采樣頻率， ts 為采樣間隔。因為最小的頻率差值為7khz6.5khz=0.5khz， 采樣頻率為 32k,所以最小的樣本數目應當是32/0.5=64 個，

6、當采樣點不足時，必然發(fā)生混疊失真，即不能觀測出原信號的頻率分布。不是采樣點數越多，頻率分辨力就越高，因為一段數據拿來就確定了時間t，注意： f0=1/t，而 t=nts，增加 n 必然減小 ts ，因此，增加 n 時 f0 是不變的。只有增加點數的同時導致增加了數據長度t 才能使分辨率越好。 還有容易搞混的一點，我們在做dft 時，常常 在有效數據后面補零達到對頻譜做某種改善的目的，我們常常認為這是增加了n，從而使頻率分辨率變好了，其實不是這樣的，補零并沒有增加有效數據的長度，仍然為 t。但是補零其實有其他好處：1.使數據 n 為 2 的整次冪，便于使用fft。2.補零后，其實是對dft 結果

7、做了插值，使譜外觀平滑化。3.由于對時域數據的截短必然造成頻譜泄露，因此在頻譜中可能出現難以辨認的譜峰，補零在一定程度上能消除這種現象。那么選擇 dft 時 n 參數要注意 ：1.由采樣定理： fs=2fh；2.頻率分辨率 .f0=fs/n，所以一般情況給定了fh 和 f0 時也就限制了 n 范圍：n=fs/f0。當取樣點 n=1664 時，左上角的信號的 16 點 fft 變換所得頻譜圖上幾乎反應出信號的任何頻率信息，因為幅值的大小并不能代表對應的頻率即為信號的頻率分布。當采樣點 n=1664 時，可以從右上角的256 點 fft 變換頻譜圖中清楚的看出吸納后的頻率成分（6.5khz,7kh

8、z,9khz ）,此時的頻譜圖就是高分辨率頻譜圖。當取樣點足夠時，只是改變fft 變換的點數，在原有高分辨率頻譜圖上增加一些疊加的頻率成分，并且在信號頻率的點上又疊加了一些新的幅值，所以，信號的的頻率成分又在原有的高分辨率頻譜的基礎上變成了高密度的頻譜圖。實驗程序：n=16;l=256;fs=32000;f1=6500;f2=7000;f3=9000;n2=512; t=1/fs;ws=2*pi*fs; n=0:n-1; x=cos(2*pi*f1*n*t)+5*cos(2*pi*f2*n*t)+cos(2*pi*f3*n*t); x1=x(1:n); x1=fft(x1,n); w=(0:n

9、-1)*ws/n)/(2*pi); %橫坐標范圍subplot(2,2,1);stem(w,abs(x1); ylabel(幅頻特性曲線 );xlabel( 采樣點為 16 的 16 點 fft 頻率 hz); x2=x(1:n) zeros(1,l-n); x2=fft(x2,l); w=(0:l-1)*ws/l)/(2*pi); subplot(2,2,3);stem(w,abs(x2); ylabel(幅頻特性曲線 );xlabel( 采樣點為 16 補零到 256 后的 256fft); n=0:l-1;華 北 電 力 大 學實 驗 報 告第 4 頁實驗名稱設計性實驗二： iir 數字

10、濾波器的設計實驗目的1、本實驗為設計性實驗。2、掌握用雙線性變換法設計iir 數字濾波器的基本原理和設計方法。3、掌握用雙線性變換法設計iir 數字 butterworth 濾波器的原理和設計方法。實驗內容iir 數字濾波器的設計用雙線性變換法設計一個iir 數字 butterworth低通濾波器。技術指標為：通帶截止頻率 fp=1khz ，阻帶截止頻率fs=1.5khz ，通帶衰減 rp1db，阻帶衰減 rs 40db ，采樣頻率 fs=10khz 。繪出濾波器的幅頻特性曲線和相頻特性曲線，判斷設計是否符合要求。題目分析：雙線性變換法 ，其主要的 優(yōu)點：避免了頻率響應的混疊現象。這是因為s

11、平面與 z 平面是單值的一一對應關系。 s平面整個 j軸單值地對應于z 平面單位圓一周，即頻率軸是單值變換關系。在零頻率附近，模擬角頻率與數字頻率之間的變換關系接近于線性關系；但當進一步增加時，增長得越來越慢，最后當時，終止在折疊頻率 =處，因而雙線性變換就不會出現由于高頻部分超過折疊頻率而混淆到低頻部分去的現象，從而消除了頻率混疊現象。但是雙線性變換的這個特點是靠頻率的嚴重非線性關系而得到的。由于這種頻率之間的非線性變換關系，就產生了新的問題。首先，一個線性相位的模擬濾波器經雙線性變換后得到非線性相位的數字濾波器，不再保持原有的線性相位了；其次，這種非線性關系要求模擬濾波器的幅頻響應必須是分

12、段常數型的，即某一頻率段的幅頻響應近似等于某一常數（這正是一般典型的低通、高通、帶通、帶阻型濾波器的響應特性），不然變換所產生的數字濾波器幅頻響應相對于原模擬濾波器的幅頻響應會有畸變，雙線性變換法幅度和相位特性的非線性映射對于分段常數的濾波器，雙線性變換后 ，仍得到幅頻特性為分段常數的濾波器，但是各個分段邊緣的臨界頻率點產生了畸變，這種 頻率的畸變 ，可以通過頻率的 預畸來加以校正 。預畸變就是： 將臨界模擬頻率事先加以畸變，然后經變換后正好映射到所需要的數字頻率上。由于用脈沖響應不變法設計iir 數字濾波器，存在從s 平面到 z 平面的多值映射的缺點，而且還有頻譜混疊現象，所以本次課程設計就

13、采用雙線性變換法來設計 iir 數字濾波器。設計濾波器的步驟： 1.得到數字指標（ wn等） 2.雙線性變換為模擬低通指標 3.歸一化模擬擬指標 4.利用通帶衰減與阻帶衰減的值求butterworth的階數 n及歸一化的模擬系統(tǒng)函數華 北 電 力 大 學實 驗 報 告第 5 頁實驗結果及分析 5.將這個再經過去歸一化得到想要的濾波器類型 6.用雙線性變換法變?yōu)閿底譃V波器程序代碼：fp=1000;%通帶截止頻率fs=1500;%阻帶截止頻率rp=1;rs=40;fs=10000; wp=2*pi*fp/fs;ws=2*pi*fs/fs; t=1;fs=1/t;wp=(2/t)*tan(wp/2)

14、;ws=(2/t)*tan(ws/2);%數 字 域 轉 模 擬域 % 預畸變n,wn=buttord(wp,ws,rp,rs,s);%求得濾波器的階數n與 3db點對應的頻率 wn z,p,k=buttap(n);%建立傳遞函數 ha(s) 。bz,az=zp2tf(z,p,k); b,a=lp2lp(bz,az,wn); b,a=bilinear(b,a,fs); h,w=freqz(b,a,512);%從幫助里面找到返回的單位是hz，而不是 rad/s h=(h(1:501);w=(w(1:501); mag=abs(h);pha=angle(h); subplot(2,1,1);plo

15、t(w/pi,mag);ylabel(幅度);xlabel(以pi為單位的頻率); title(幅頻特性曲線 );grid; subplot(2,1,2);plot(w/pi,pha);title(相頻特性曲線 ); ylabel(相位);xlabel(以pi為單位的頻率 ); grid; 華 北 電 力 大 學實 驗 報 告第 6 頁實驗名稱設計性實驗三：fir 數字濾波器的設計實驗目的1、本實驗為設計性實驗。2、掌握用窗函數法設計fir 數字濾波器的基本原理和設計方法。3、掌握用窗函數法設計線性相位fir 數字低通濾波器的編程實現。實驗內容fir 數字濾波器的設計用窗函數法設計一個線性相位

16、fir 數字低通濾波器。技術指標為：通帶截止角頻率 p=0.2，阻帶截止角頻率 s=0.3 ，通帶衰減 rp1db ，阻帶衰減 rs 40db 。繪出濾波器的幅頻特性曲線和相頻特性曲線，判斷設計是否符合要求。根據相同的濾波器要求， 選用不同的窗函數進行設計， 比較各種窗函數對 fir數字濾波器頻率特性的影響。實驗分析：我們考慮窗函數主要是以下幾點：1.主瓣 寬度 b最小 （相當于矩形窗時的 4/n， 頻域兩個過零點間的寬度） 。2.最大旁瓣 峰值 a 最?。ㄟ@樣旁瓣泄露小，一些高頻分量損失少了）。3.旁瓣 譜峰漸近 衰減速度 d 最大（同樣是減少旁瓣泄露） 。在此，總結幾種很常用的窗函數的優(yōu)缺

17、點：矩形窗： b=4/n a=-13db d=-6db/oct 三角窗： b=8/n a=-27db d=-12db/oct 漢寧窗： b=8/n a=-32db d=-18db/oct 海明窗： b=8/n a=-43db d=-6db/oct 布萊克曼窗： b=12/n a=-58db d=-18db/oct 可以看出，矩形窗有最窄的主瓣，但是旁瓣泄露嚴重。漢寧窗和海明窗雖主瓣較寬，但是旁瓣泄露少，是常選用的窗函數。實驗中將用數據驗證幾種窗函數的特性，選擇了：矩形窗、布萊克曼窗、凱澤窗設計核心 是從給定的頻率特性通過加窗確定有限長單位脈沖響應序列h(n) 利用窗函數設計濾波器，是讓待設計的

18、濾波器去逼近理想特性，理想低通濾波器的頻率特性應該是振幅特性在通帶內為1，阻帶內為 0，在通帶內的相位特性與 w 成線性關系。設計步驟 ：1.確定數字濾波器的性能要求，確定濾波器的臨界頻率，濾波器的長度。2.根據性能要求， 合理選擇單位脈沖響應， 從而確定理想頻率響應的幅頻特性和相頻特性3.選擇適合的窗函數起初所需設計的fir 濾波器的單位脈沖響應。若不滿足，可以改變傳函數的形式或者是長度n 華 北 電 力 大 學實 驗 報 告第 7 頁實驗結果及分析布萊克曼窗實驗程序：wp=0.2*pi;ws=0.3*pi;w=ws-wp;n=ceil(11*pi/w)+1;n=0:1:n-1; wc=(w

19、s+wp)/2;n1=(n-1)/2;n=0:1:n-1;m=n-n1+eps; hd=sin(wc*m)./(pi*m);stem(n,hd);w_bla=(blackman(n); stem(n,w_bla);h=hd.*w_bla;stem(n,h);h,w=freqz(h); mag=abs(h);db=20*log10(mag); subplot(1,2,1);plot(w/pi,db);title( 幅度響應 db); xlabel(以pi 為單位的頻率 );ylabel(對數幅度 ); grid on;pha=angle(h); subplot(1,2,2);plot(w/pi,

20、pha);title(相位響應 ); xlabel(以pi 為單位的頻率 );ylabel(相位);grid on; delta_w=2*pi/1000;rp=-(min(db(1:1:wp/delta_w+1) as=-round(max(db(ws/delta_w+1:1:501)legend(布萊克曼窗 ) rp = 0.0011 as = 53 從結果可以看出， rp、 as 都滿足實驗要求，旁瓣幅度較小，且衰減較快華 北 電 力 大 學實 驗 報 告第 8 頁實驗結果及分析矩形窗實驗程序：wp=0.2*pi;ws=0.3*pi;w=ws-wp;n=ceil(1.8*pi/w)+1;

21、n=0:1:n-1;wc=(ws+wp)/2;n1=(n-1)/2;n=0:1:n-1; m=n-n1+eps;hd=sin(wc*m)./(pi*m);stem(n,hd);w_rec=(rectwin(n); stem(n,w_rec);h=hd.*w_rec;stem(n,h);h,w=freqz(h);mag=abs(h); db=20*log10(mag); subplot(1,2,1);plot(w/pi,db);title( 幅度響應 db); xlabel(以pi 為單位的頻率 );ylabel(對數幅度 ); grid on;pha=angle(h); subplot(1,2

22、,2);plot(w/pi,pha);title(相位響應 ); xlabel(以pi 為單位的頻率 );ylabel(相位);grid on; delta_w=2*pi/1000;rp=-(min(db(1:1:wp/delta_w+1) as=-round(max(db(ws/delta_w+1:1:501)legend(矩形窗 ) rp = 0.6459 as = 17 從結果可以看出， rp 滿足實驗要求，但 as 并不滿足實驗要求。旁瓣幅度較大，且衰減較慢華 北 電 力 大 學實 驗 報 告第 9 頁實驗結果及分析實驗程序：%用 kaiser 窗函數設計 fir 數字濾波器wp=0.2*pi;ws=0.3*pi;as=40;tr_width=ws-wp; n=ceil(as-7.95)/(14.36*tr_width/(2*pi)+1)+1;n=0:1:n-1; beta=0.1102*(as-8.7)wc=(ws+wp)/2;hd=ideal_1p(wc,n); w_kai=(kaiser(n,beta);h=hd.*w_kaidb,mag,pha,grd,w=freqz_m(h,1); delta_w=2*pi/1000;rp=-(min(db(1:1:wp/delta_w+1) as=-round(ma