東南大學(xué)信息學(xué)院DSP課程第一次實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告

1、數(shù)字信號(hào)處理第一次實(shí)驗(yàn)報(bào)告 -快速傅里葉變化及其應(yīng)用姓名：印友進(jìn)學(xué)號(hào)：04012540一:實(shí)驗(yàn)?zāi)康募由顚?duì)FFT的理解，熟悉matlab中的有關(guān)函數(shù)。應(yīng)用FFT對(duì)典型信號(hào)進(jìn)行頻譜分析。了解應(yīng)用FFT進(jìn)行信號(hào)頻譜分析過(guò)程中可能出現(xiàn)的問(wèn)題，以便在實(shí)際中正確應(yīng)用FFT.應(yīng)用FFT實(shí)現(xiàn)序列的線性卷積和相關(guān)。二:實(shí)驗(yàn)原理:在各種信號(hào)序列中，有限長(zhǎng)序列信號(hào)處理占有很重要地位，對(duì)有限長(zhǎng)序列，我們可以使用離散Fouier變換(DFT)。這一變換不但可以很好的反映序列的頻譜特性，而且易于用快速算法在計(jì)算機(jī)上實(shí)現(xiàn)，當(dāng)序列x(n)的長(zhǎng)度為N時(shí)，它的DFT定義為： 反變換為： 有限長(zhǎng)序列的DFT是其Z變換在單位圓上的等

2、距采樣，或者說(shuō)是序列Fourier變換的等距采樣，因此可以用于序列的譜分析。 FFT并不是與DFT不同的另一種變換，而是為了減少DFT運(yùn)算次數(shù)的一種快速算法。它是對(duì)變換式進(jìn)行一次次分解，使其成為若干小點(diǎn)數(shù)的組合，從而減少運(yùn)算量。常用的FFT是以2為基數(shù)的，其長(zhǎng)度 。它的效率高，程序簡(jiǎn)單，使用非常方便，當(dāng)要變換的序列長(zhǎng)度不等于2的整數(shù)次方時(shí)，為了使用以2為基數(shù)的FFT，可以用末位補(bǔ)零的方法，使其長(zhǎng)度延長(zhǎng)至2的整數(shù)次方。 (一)在運(yùn)用DFT進(jìn)行頻譜分析的過(guò)程中可能的產(chǎn)生三種誤差(1)混疊 序列的頻譜是被采樣信號(hào)的周期延拓，當(dāng)采樣速率不滿足Nyquist定理時(shí)，就會(huì)發(fā)生頻譜混疊，使得采樣后的信號(hào)序列

3、頻譜不能真實(shí)的反映原信號(hào)的頻譜。避免混疊現(xiàn)象的唯一方法是保證采樣速率足夠高，使頻譜混疊現(xiàn)象不致出現(xiàn)，即在確定采樣頻率之前，必須對(duì)頻譜的性質(zhì)有所了解，在一般情況下，為了保證高于折疊頻率的分量不會(huì)出現(xiàn)，在采樣前，先用低通模擬濾波器對(duì)信號(hào)進(jìn)行濾波。(2)泄漏實(shí)際中我們往往用截短的序列來(lái)近似很長(zhǎng)的甚至是無(wú)限長(zhǎng)的序列，這樣可以使用較短的DFT來(lái)對(duì)信號(hào)進(jìn)行頻譜分析，這種截短等價(jià)于給原信號(hào)序列乘以一個(gè)矩形窗函數(shù)，也相當(dāng)于在頻域?qū)⑿盘?hào)的頻譜和矩形窗函數(shù)的頻譜卷積，所得的頻譜是原序列頻譜的擴(kuò)展。泄漏不能與混疊完全分開(kāi)，因?yàn)樾孤?dǎo)致頻譜的擴(kuò)展，從而造成混疊。為了減少泄漏的影響，可以選擇適當(dāng)?shù)拇昂瘮?shù)使頻譜的擴(kuò)散減至

4、最小。(3)柵欄效應(yīng)DFT是對(duì)單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個(gè)連續(xù)函數(shù)，就一定意義上看，用DFT來(lái)觀察頻譜就好像通過(guò)一個(gè)柵欄來(lái)觀看一個(gè)圖景一樣，只能在離散點(diǎn)上看到真實(shí)的頻譜，這樣就有可能發(fā)生一些頻譜的峰點(diǎn)或谷點(diǎn)被“尖樁的柵欄”所攔住，不能別我們觀察到。 減小柵欄效應(yīng)的一個(gè)方法就是借助于在原序列的末端填補(bǔ)一些零值，從而變動(dòng)DFT的點(diǎn)數(shù)，這一方法實(shí)際上是人為地改變了對(duì)真實(shí)頻譜采樣的點(diǎn)數(shù)和位置，相當(dāng)于搬動(dòng)了每一根“尖樁柵欄”的位置，從而使得頻譜的峰點(diǎn)或谷點(diǎn)暴露出來(lái)。(二)用FFT計(jì)算線性卷積 用FFT可以實(shí)現(xiàn)兩個(gè)序列的圓周卷積。在一定的條件下，可以使圓周卷積等于線性卷積。一般情況

5、，設(shè)兩個(gè)序列的長(zhǎng)度分別為N1和N2，要使圓周卷積等于線性卷積的充要條件是FFT的長(zhǎng)度NN1N2對(duì)于長(zhǎng)度不足N的兩個(gè)序列，分別將他們補(bǔ)零延長(zhǎng)到N。 當(dāng)兩個(gè)序列中有一個(gè)序列比較長(zhǎng)的時(shí)候，我們可以采用分段卷積的方法。有兩種方法：（1）重疊相加法。將長(zhǎng)序列分成與短序列相仿的片段，分別用FFT對(duì)它們作線性卷積，再將分段卷積各段重疊的部分相加構(gòu)成總的卷積輸出。（2）重疊保留法。這種方法在長(zhǎng)序列分段時(shí)，段與段之間保留有互相重疊的部分，在構(gòu)成總的卷積輸出時(shí)只需將各段線性卷積部分直接連接起來(lái)，省掉了輸出段的直接相加。（三）用FFT計(jì)算相關(guān)函數(shù) 兩個(gè)長(zhǎng)為的實(shí)離散時(shí)間序列與的互相關(guān)函數(shù)定義為 ：的離散傅里葉變換為：

6、 當(dāng)時(shí)，得到的自相關(guān)函數(shù)為： 利用FFT求兩個(gè)有限長(zhǎng)序列線性相關(guān)的步驟（設(shè)長(zhǎng)，長(zhǎng)）： （1）為了使兩個(gè)有限長(zhǎng)序列的線性相關(guān)可用其圓周相關(guān)代替而不產(chǎn)生混淆，選擇周期，以便使用FFT，將，補(bǔ)零至長(zhǎng)為。（2）用FFT計(jì)算（3）（4）對(duì)作IFFT；取后項(xiàng)，得；取前項(xiàng)，得。 三、實(shí)驗(yàn)內(nèi)容及步驟實(shí)驗(yàn)中用到的信號(hào)序列：a)高斯（Gaussian）序列 b)衰減正弦序列 c)三角波序列 d)反三角波序列 上機(jī)實(shí)驗(yàn)內(nèi)容：（1）觀察高斯序列的時(shí)域和幅頻特性，固定信號(hào)xa(n)中參數(shù)p=8，改變q的值，使q分別等于2，4，8，觀察它們的時(shí)域和幅頻特性，了解當(dāng)q取不同值時(shí)，對(duì)信號(hào)序列的時(shí)域幅頻特性的影響；固定q=8，

7、改變p，使p分別等于8，13，14，觀察參數(shù)p變化對(duì)信號(hào)序列的時(shí)域及幅頻特性的影響，觀察p等于多少時(shí)，會(huì)發(fā)生明顯的泄漏現(xiàn)象，混疊是否也隨之出現(xiàn)？記錄實(shí)驗(yàn)中觀察到的現(xiàn)象，繪出相應(yīng)的時(shí)域序列和幅頻特性曲線。參數(shù)p=8不變分析：保持p=8不變，改變q，時(shí)域上的波形展開(kāi)，幅度最大值不變（p=8時(shí)幅度為1）。頻域上頻譜壓縮，幅度增大。分析：保持q=8不變，改變p：時(shí)域上，波形發(fā)生偏移（幅度并未變化）。 當(dāng)p=13、14時(shí)出現(xiàn)頻譜泄漏和混疊，p=14時(shí)泄漏現(xiàn)象明顯。這是因?yàn)闀r(shí)域信號(hào)被截?cái)啵囟痰男盘?hào)并不能如實(shí)反映原信號(hào)。信號(hào)截短的過(guò)程中矩形函數(shù)頻譜的作用下，使得原有的信號(hào)頻譜出現(xiàn)了較大的波動(dòng)和分散的頻譜

8、擴(kuò)展，即產(chǎn)生了頻譜泄漏。頻譜泄漏也會(huì)引起混疊，由于泄漏使信號(hào)的頻譜展寬，如果它的高頻成分超過(guò)了折疊頻率（fs/2），就造成了混疊，這種混疊的可能性在矩形窗截短時(shí)尤為明顯。源程序：p=8;q=2;n=0:15;x=exp(-(n-p).2/q);w=0:0.01:2*pi;F=zeros(size(w);for k=1:16 F=F+x(k)*exp(-i*w*k);endf=fft(x);subplot(3,3,1);stem(n,x);xlabel(n);ylabel(時(shí)域特性);title(p=8 q=2);subplot(3,3,2);plot(w,abs(F);xlabel(w);yl

9、abel(DTFT頻域特性);title(p=8 q=2);subplot(3,3,3);stem(n,abs(f);xlabel(n);ylabel(FFT頻域特性);title(p=8 q=2);p=8;q=4;n=0:15;x=exp(-(n-p).2/q);w=0:0.01:2*pi;F=zeros(size(w);for k=1:16 F=F+x(k)*exp(-i*w*k);endf=fft(x);subplot(3,3,4);stem(n,x);xlabel(n);ylabel(時(shí)域特性);title(p=8 q=4);subplot(3,3,5);plot(w,abs(F);x

10、label(w);ylabel(DTFT頻域特性);title(p=8 q=4);subplot(3,3,6);stem(n,abs(f);xlabel(k);ylabel(FFT頻域特性);title(p=8 q=4);p=8;q=8;n=0:15;x=exp(-(n-p).2/q);w=0:0.01:2*pi;F=zeros(size(w);for k=1:16 F=F+x(k)*exp(-i*w*k);endf=fft(x);subplot(3,3,7);stem(n,x);xlabel(n);ylabel(時(shí)域特性);title(p=8 q=8);subplot(3,3,8);plot

11、(w,abs(F);xlabel(w);ylabel(DTFT頻域特性);title(p=8 q=8);subplot(3,3,9);stem(n,abs(f);xlabel(k);ylabel(FFT頻域特性);title(p=8 q=8);p=13;q=8;n=0:15;x=exp(-(n-p).2/q);w=-pi:0.01:pi;F=zeros(size(w);for k=1:16 F=F+x(k)*exp(-i*w*k);endf=fft(x);subplot(3,3,4);stem(n,x);xlabel(n);ylabel(時(shí)域特性);title(p=13 q=8);subplo

12、t(3,3,5);plot(w,abs(F);xlabel(w);ylabel(DTFT頻域特性);title(p=13 q=8);subplot(3,3,6);stem(n,abs(f);xlabel(k);ylabel(FFT頻域特性);title(p=13 q=8);p=14;q=8;n=0:15;x=exp(-(n-p).2/q);w=-pi:0.01:pi;F=zeros(size(w);for k=1:16 F=F+x(k)*exp(-i*w*k);endf=fft(x);subplot(3,3,7);stem(n,x);xlabel(n);ylabel(時(shí)域特性);title(p

13、=14 q=8);subplot(3,3,8);plot(w,abs(F);xlabel(w);ylabel(DTFT頻域特性);title(p=14 q=8);subplot(3,3,9);stem(n,abs(f);xlabel(k);ylabel(FFT頻域特性);title(p=14 q=8);(2)觀察衰減正弦序列xb(n)的時(shí)域和幅頻特性，a=0.1，f=0.0625，檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出幅頻特性曲線，改變f，使f分別等于0.4375和0.5625，觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)位置，有無(wú)混疊和泄漏現(xiàn)象？說(shuō)明產(chǎn)生現(xiàn)象的原因。分析：保持a不變，當(dāng)f很

14、小時(shí)，頻譜會(huì)在正負(fù)頻率上產(chǎn)生相對(duì)稱(chēng)的峰值（圖中取的是0到2之間的頻譜圖）。當(dāng)f=0.4375,頻譜發(fā)生了混疊現(xiàn)象，峰值有所偏移。 當(dāng)f=0.5625,不滿足取樣定理（1=N1+N2-1）。源程序：n=0:15;m=0:30; %L=16+16-1=31線性長(zhǎng)度；p=8;q=2;xa=exp(-(n-p).2/q);a=0.1;f=0.0625;xb=exp(-a*n).*sin(2*pi*f*n);fa=fft(xa,16);fa_2=fft(xa,31);fb=fft(xb,16);fb_2=fft(xb,31);F_line=ifft(fa.*fb);F_cirle=ifft(fa_2.*

15、fb_2);subplot(2,1,1)stem(n,F_line);title(循環(huán)卷積);subplot(2,1,2)stem(m,F_cirle);title(線性卷積);（6）產(chǎn)生一512點(diǎn)的隨機(jī)序列xe(n),并用xc(n)和xe(n)作線性卷積，觀察卷積前后xe(n)頻譜的變化。要求將xe(n)分成8段，分別采用重疊相加法和重疊保留法。源程序：clearxe=rand(1,512)-0.5;for n=0:3 xc(n+1)=n;endfor n=4:7 xc(n+1)=8-n;endG=fft(xc,71);%512/8+8-1=71;w=-pi:0.01:pi;F=zeros(

16、size(w);for r=1:512 F=F+xe(r)*exp(-i*w*r);endsubplot(3,1,1);plot(w,abs(F);title(xe頻域特性);%重疊相加法for m=1:8 f_64=xe(64*(m-1)+1):64*m);%每段寫(xiě)入64個(gè)隨機(jī)值 k(m,:)=ifft(fft(f_64,71).*G);%每個(gè)xe（m）與x(c)(補(bǔ)零后)的卷積endy(1:64)=k(1,1:64);%第一段前64個(gè)for n=1:7 y(64*n+1:64*n+7)=k(n,65:71)+k(n+1,1:7);%從第2段起每段重疊部分的值為前后段相加； y(64*n+8

17、:64*n+71)=k(n+1,8:71);%其余正常賦值；endY=fft(y);i=0:(512+8-1)-1);subplot(3,1,2);stem(i,y);title(重疊相加法的時(shí)域特性);subplot(3,1,3);stem(i,abs(Y);title(重疊相加法的頻域特性);源程序：clearxe=rand(1,512)-0.5;for n=0:3 xc(n+1)=n;endfor n=4:7 xc(n+1)=8-n;endG=fft(xc,71);w=-pi:0.01:pi;F=zeros(size(w);for r=1:512 F=F+xe(r)*exp(-i*w*r

18、);endsubplot(3,1,1);plot(w,abs(F);title(xe頻域特性);%重疊保留法for m=1:8 if(m=1) f(1:7)=0;f(8:71)=xe(1:64); k(m,:)=ifft(fft(f,71).*G); else f=xe(64*(m-1)-6:64*(m-1)+64); k(m,:)=ifft(fft(f,71).*G); endend y(1:71)=k(1,1:71);for d=2:8 y(64*(d-1)+8:64*(d-1)+71)=k(d,8:71);endY=fft(y);i=0:(512+8-1)-1);subplot(3,1,

19、2);stem(i,y);title(重疊保留法的時(shí)域特性);subplot(3,1,3);stem(i,abs(Y);title(重疊保留法的頻域特性);（7）用FFT分別計(jì)算xa(n)(p=8,q=2)和xb(n)（a0.1，f0.0625）的16點(diǎn)循環(huán)相關(guān)和線性相關(guān)，問(wèn)一共有多少種結(jié)果，它們之間有何共同點(diǎn)？源程序：n=0:15;p=8;q=2;xa=exp(-(n-p).2/q);a=0.1;f=0.0625;xb=exp(-a*n).*sin(2*pi*f*n);fa=fft(xa);fb=fft(xb);rm=real(ifft(conj(fa).*fb);subplot(2,2,1

20、);stem(n,rm);xlabel(m);ylabel(幅度);title(循環(huán)相關(guān)1);fa_2=fft(xa,2*16);fb_2=fft(xb,2*16);rm_2=real(ifft(conj(fa_2).*fb_2);rm_2=rm_2(16+2):(16*2) rm_2(1:16);m=(-16+1):(16-1);subplot(2,2,2);stem(m,rm_2);xlabel(m);ylabel(幅度);title(線性相關(guān)1);RM=real(ifft(conj(fb).*fa);%相關(guān)函數(shù)公式中fa、fb調(diào)換；subplot(2,2,3);stem(n,RM);xl

21、abel(m);ylabel(幅度);title(循環(huán)相關(guān)2);RM_2=real(ifft(conj(fb_2).*fa_2);RM_2=RM_2(16+2):(16*2) RM_2(1:16);subplot(2,2,4);stem(m,RM_2);xlabel(m);ylabel(幅度);title(線性相關(guān)2);（8）用FFT分別計(jì)算xa(n)(p=8,q=2)和xb(n)（a0.1，f0.0625）的自相關(guān)函數(shù)。源程序：n=0:15;p=8;q=2;xa=exp(-(n-p).2/q);a=0.1;f=0.0625;xb=exp(-a*n).*sin(2*pi*f*n);fa=fft(xa);fb=fft(xb);rm=real(ifft(conj(fa).*fa);subplot(2,2,1);stem(n,rm);xla