數(shù)字信號處理課程設(shè)計(jì)實(shí)驗(yàn)報(bào)告(20211126002056)

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

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

3、 m，n 的序列的線性卷積可用長度均為l 的循環(huán)卷積來代替，用循環(huán)卷積計(jì)算線性卷積的方法歸納如下: 將長為 n2的序列 x(n) 延長到 l, 補(bǔ) l-n2 個零將長為 n1的序列 h(n) 延長到 l, 補(bǔ) l-n1 個零如果 ln1+n2-1,則循環(huán)卷積與線性卷積相等, 此時(shí), 可有 fft 計(jì)算線性卷積, 方法如下 : a.計(jì)算 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 可見，只要進(jìn)行二次fft,一次 ifft 就可完成線性卷積計(jì)算。計(jì)算表明,l32 時(shí), 上述計(jì)算線性卷

4、積的方法比直接計(jì)算線卷積有明顯的優(yōu)越性, 因此, 也稱上述循環(huán)卷積方法為快速卷積法。題目利用 fft 計(jì)算線性卷積以下為程序及結(jié)果：華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 2 頁實(shí)驗(yàn)結(jié)果及分析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) %補(bǔ)零操作h2=ones(1,8) zeros(1,15) x3=fft(x2) %

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

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

7、做了插值，使譜外觀平滑化。3.由于對時(shí)域數(shù)據(jù)的截短必然造成頻譜泄露，因此在頻譜中可能出現(xiàn)難以辨認(rèn)的譜峰，補(bǔ)零在一定程度上能消除這種現(xiàn)象。那么選擇 dft 時(shí) n 參數(shù)要注意 ：1.由采樣定理： fs=2fh；2.頻率分辨率 .f0=fs/n，所以一般情況給定了fh 和 f0 時(shí)也就限制了 n 范圍：n=fs/f0。當(dāng)取樣點(diǎn) n=1664 時(shí)，左上角的信號的 16 點(diǎn) fft 變換所得頻譜圖上幾乎反應(yīng)出信號的任何頻率信息，因?yàn)榉档拇笮〔⒉荒艽韺?yīng)的頻率即為信號的頻率分布。當(dāng)采樣點(diǎn) n=1664 時(shí)，可以從右上角的256 點(diǎn) fft 變換頻譜圖中清楚的看出吸納后的頻率成分（6.5khz,7kh

8、z,9khz ）,此時(shí)的頻譜圖就是高分辨率頻譜圖。當(dāng)取樣點(diǎn)足夠時(shí)，只是改變fft 變換的點(diǎn)數(shù)，在原有高分辨率頻譜圖上增加一些疊加的頻率成分，并且在信號頻率的點(diǎn)上又疊加了一些新的幅值，所以，信號的的頻率成分又在原有的高分辨率頻譜的基礎(chǔ)上變成了高密度的頻譜圖。實(shí)驗(yàn)程序：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); %橫坐標(biāo)范圍subplot(2,2,1);stem(w,abs(x1); ylabel(幅頻特性曲線 );xlabel( 采樣點(diǎn)為 16 的 16 點(diǎn) 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( 采樣點(diǎn)為 16 補(bǔ)零到 256 后的 256fft); n=0:l-1;華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 4 頁實(shí)驗(yàn)名稱設(shè)計(jì)性實(shí)驗(yàn)二： iir 數(shù)字

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

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

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

13、采用雙線性變換法來設(shè)計(jì) iir 數(shù)字濾波器。設(shè)計(jì)濾波器的步驟： 1.得到數(shù)字指標(biāo)（ wn等） 2.雙線性變換為模擬低通指標(biāo) 3.歸一化模擬擬指標(biāo) 4.利用通帶衰減與阻帶衰減的值求butterworth的階數(shù) n及歸一化的模擬系統(tǒng)函數(shù)華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 5 頁實(shí)驗(yàn)結(jié)果及分析 5.將這個再經(jīng)過去歸一化得到想要的濾波器類型 6.用雙線性變換法變?yōu)閿?shù)字濾波器程序代碼：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);%數(shù) 字 域 轉(zhuǎn) 模 擬域 % 預(yù)畸變n,wn=buttord(wp,ws,rp,rs,s);%求得濾波器的階數(shù)n與 3db點(diǎn)對應(yīng)的頻率 wn z,p,k=buttap(n);%建立傳遞函數(shù) 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; 華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 6 頁實(shí)驗(yàn)名稱設(shè)計(jì)性實(shí)驗(yàn)三：fir 數(shù)字濾波器的設(shè)計(jì)實(shí)驗(yàn)?zāi)康?、本實(shí)驗(yàn)為設(shè)計(jì)性實(shí)驗(yàn)。2、掌握用窗函數(shù)法設(shè)計(jì)fir 數(shù)字濾波器的基本原理和設(shè)計(jì)方法。3、掌握用窗函數(shù)法設(shè)計(jì)線性相位fir 數(shù)字低通濾波器的編程實(shí)現(xiàn)。實(shí)驗(yàn)內(nèi)容fir 數(shù)字濾波器的設(shè)計(jì)用窗函數(shù)法設(shè)計(jì)一個線性相位

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

17、點(diǎn)：矩形窗： 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 可以看出，矩形窗有最窄的主瓣，但是旁瓣泄露嚴(yán)重。漢寧窗和海明窗雖主瓣較寬，但是旁瓣泄露少，是常選用的窗函數(shù)。實(shí)驗(yàn)中將用數(shù)據(jù)驗(yàn)證幾種窗函數(shù)的特性，選擇了：矩形窗、布萊克曼窗、凱澤窗設(shè)計(jì)核心 是從給定的頻率特性通過加窗確定有限長單位脈沖響應(yīng)序列h(n) 利用窗函數(shù)設(shè)計(jì)濾波器，是讓待設(shè)計(jì)的

18、濾波器去逼近理想特性，理想低通濾波器的頻率特性應(yīng)該是振幅特性在通帶內(nèi)為1，阻帶內(nèi)為 0，在通帶內(nèi)的相位特性與 w 成線性關(guān)系。設(shè)計(jì)步驟 ：1.確定數(shù)字濾波器的性能要求，確定濾波器的臨界頻率，濾波器的長度。2.根據(jù)性能要求， 合理選擇單位脈沖響應(yīng)， 從而確定理想頻率響應(yīng)的幅頻特性和相頻特性3.選擇適合的窗函數(shù)起初所需設(shè)計(jì)的fir 濾波器的單位脈沖響應(yīng)。若不滿足，可以改變傳函數(shù)的形式或者是長度n 華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 7 頁實(shí)驗(yàn)結(jié)果及分析布萊克曼窗實(shí)驗(yàn)程序：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( 幅度響應(yīng) db); xlabel(以pi 為單位的頻率 );ylabel(對數(shù)幅度 ); grid on;pha=angle(h); subplot(1,2,2);plot(w/pi,

20、pha);title(相位響應(yīng) ); 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 從結(jié)果可以看出， rp、 as 都滿足實(shí)驗(yàn)要求，旁瓣幅度較小，且衰減較快華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 8 頁實(shí)驗(yàn)結(jié)果及分析矩形窗實(shí)驗(yàn)程序：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( 幅度響應(yīng) db); xlabel(以pi 為單位的頻率 );ylabel(對數(shù)幅度 ); grid on;pha=angle(h); subplot(1,2

22、,2);plot(w/pi,pha);title(相位響應(yīng) ); 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 從結(jié)果可以看出， rp 滿足實(shí)驗(yàn)要求，但 as 并不滿足實(shí)驗(yàn)要求。旁瓣幅度較大，且衰減較慢華 北 電 力 大 學(xué)實(shí) 驗(yàn) 報(bào) 告第 9 頁實(shí)驗(yàn)結(jié)果及分析實(shí)驗(yàn)程序：%用 kaiser 窗函數(shù)設(shè)計(jì) fir 數(shù)字濾波器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