數(shù)字信號(hào)處理實(shí)驗(yàn)第一次報(bào)告實(shí)驗(yàn)三-快速傅立葉變換及其應(yīng)用

1、數(shù)字信號(hào)處理Matlab實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)三 快速傅立葉變換及其應(yīng)用姓名： 學(xué)號(hào)： 一 實(shí)驗(yàn)平臺(tái)二 實(shí)驗(yàn)?zāi)康模海?） 在理論學(xué)習(xí)的基礎(chǔ)上，通過本實(shí)驗(yàn)，加深對(duì)FFT的理解，熟悉MATLAB中的有關(guān)函數(shù)。（2） 應(yīng)用FFT對(duì)典型信號(hào)進(jìn)行頻譜分析。（3） 了解應(yīng)用FFT進(jìn)行信號(hào)頻譜分析過程中可能出現(xiàn)的問題，以便在實(shí)際中正確應(yīng)用FFT。（4） 應(yīng)用FFT實(shí)現(xiàn)序列的線性卷積和相關(guān)。三 實(shí)驗(yàn)原理：（1） 混疊：采樣序列的頻譜是被采樣信號(hào)頻譜的周期延拓，當(dāng)采樣頻率不滿足奈奎斯特采樣定理的時(shí)候，就會(huì)發(fā)生混疊，使得刺癢后的序列信號(hào)的頻譜不能真實(shí)的反映原采樣信號(hào)的頻譜。（2） 泄露：根據(jù)理論分析，一個(gè)時(shí)間的信號(hào)其頻帶寬

2、度為無限，一個(gè)時(shí)間無限的信號(hào)其頻帶寬度則為有限。因此對(duì)一個(gè)時(shí)間有限的信號(hào)，應(yīng)用DFT進(jìn)行分析，頻譜混疊難以避免。對(duì)一個(gè)時(shí)間無限的信號(hào)雖然頻帶有限，但在實(shí)際運(yùn)算中，時(shí)間總是取有限值，在將信號(hào)截?cái)嗟倪^程中，出現(xiàn)了分散的擴(kuò)展譜線的現(xiàn)象，稱之為頻譜泄露或功率泄露。（3） 柵欄效應(yīng)：DFT是對(duì)單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個(gè)連續(xù)函數(shù)，就在一定意義上看，用DFT來觀察頻譜就好象通過一個(gè)柵欄來觀看一個(gè)景象一樣，只能在離散點(diǎn)上看到真實(shí)的頻譜，這樣就有可能發(fā)生一些頻譜的峰點(diǎn)和谷點(diǎn)被“尖樁的柵欄”所擋住，不能被我們觀察到。（4） 圓周卷積：把序列X（N）分布在N等份的圓周上，而序列Y（N）經(jīng)反

3、摺后也分布在另一個(gè)具有N等份的同心圓的圓周上。兩圓上對(duì)應(yīng)的數(shù)兩量?jī)上喑饲蠛停偷玫饺烤矸e序列。這個(gè)卷積過程稱做圓周卷積。（5） 互相關(guān)函數(shù)反映了兩個(gè)序列X（N）和Y（N）的相似程度，用FFT可以很快的計(jì)算互相關(guān)函數(shù)。四 實(shí)驗(yàn)內(nèi)容：實(shí)驗(yàn)中用到的函數(shù)序列：(a) Gaussian序列exp(-(n-p).2)/q), 0=<n=<15Xa(n)=0, 其他 （b）衰減正弦序列exp(-an)*sin(2pi*fn),0=<n=<15 X(b)= 0, 其他 (c)三角波序列 n, 0=<n=<3 Xb(n)=8-n, 4=<n=<7 0, 其他 (

4、d)反三角波序列 4-n, 0=<n=<3 Xc(n)= n-4, 4=<n=<7 0, 其他五 上機(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í)域和幅頻特性影響；改變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í)域序列和幅頻特性曲線。 程序：n=0:15;xa=exp(-(n-p).2/q);subplot(2,1,

5、1);plot(n,xa);ya=fft(xa);ya=abs(ya);subplot(2,1,2);stem(n,ya); p=8,q=2（注：上面是時(shí)域，下面是頻域） P=8,q=4 P=8,q=8 P=13,q=8 P=14,q=8結(jié)論: X(n)中的參數(shù)p為高斯序列的峰值位置，q則表示高斯序列峰的尖銳度，（即峰值邊沿的陡峭度）。q值越大，時(shí)域圖中圖象越平緩，序列變化越慢；其幅頻特性圖中高頻分量越少，頻譜越窄，越不容易產(chǎn)生混疊。p值越大，序列右移，在規(guī)定的窗口內(nèi)有效值被截?cái)嗟脑蕉?。因?yàn)榇翱诮財(cái)鄷?huì)造成窗口泄露，所以我們可以在幅頻特性圖中看到，隨著p值的變大，高頻分量會(huì)增加。易出現(xiàn)泄露，當(dāng)p

6、=13時(shí)，特別是p=14時(shí)，產(chǎn)生了明顯的泄露與混疊。2 觀察衰減正弦序列xb(n)的時(shí)域和幅頻特性，a=0.1,f=0.0625,檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出幅頻特性曲線，改變f，使f分別等于0.4375和0.5625，觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)的位置，有無混疊和泄漏現(xiàn)象？說明產(chǎn)生現(xiàn)象的原因。 程序： n=0:15; a=0.1;xb=exp(-a*n).*sin(2*pi*f*n);subplot(2,1,1); plot(n,xb);yb=fft(xb); yb=abs(yb);subplot(2,1,2);stem(n,yb);f=0.0625; f=0.4

7、375; f=0.5625結(jié)論：該實(shí)驗(yàn)中f=F/fs （F固有頻率 fs采樣頻率，統(tǒng)一做歸一化處理 fs=1）圖中的幅頻特性圖：當(dāng)f=0.0625時(shí)，沒有產(chǎn)生明顯的混疊和泄露；當(dāng)f=0.4375和f=0.5625時(shí)，產(chǎn)生了混疊，是因?yàn)椴粷M足奈奎斯特采樣定理的緣故;圖中后兩個(gè)序列的時(shí)域圖：因?yàn)?.4375+0.5625=1，滿足如下等式（此情況只適用于正弦序列），Xb(n)|f=0.4375=-Xb(n)|f=0.5625,即sin(2pi*fn）=-sin2pi(1-f)n,其幅頻特性是完全相同的。3. 觀察三角波和反三角波序列的時(shí)域和幅頻特性，用N=8點(diǎn)FFT分析信號(hào)序列xc(n)和xd(n

8、)的幅頻特性，觀察兩者的序列形狀和頻譜曲線有什么異同？繪出兩序列及其幅頻特性曲線。在xc(n)和xd(n)末尾補(bǔ)零，用N=32點(diǎn)FFT分析這兩個(gè)信號(hào)序列的幅頻特性，觀察頻譜特性發(fā)生了什么變化？?jī)煞N情況下的FFT頻譜還有相同之處嗎？這些變化說明了什么？ 程序：n1=0:3;xc1=n1;xd1=4-n1;n2=4:7;xc2=8-n2;xd2=n2-4;xc=xc1,xc2;xd=xd1,xd2;subplot(2,2,1);n1=0:31; n2=0:7;plot(n2,xc);yc=fft(xc,n);yc=abs(yc);subplot(2,2,2);stem(n1,yc);subplot

9、(2,2,3);plot(n2,xd);subplot(2,2,4);yd=fft(xd，n);yd=abs(yd);stem(n1,yd); n=8;(左邊是時(shí)域，右邊是頻域，下同) n=32;結(jié)論：反三角波的邊沿比較陡峭，因此它的幅頻特性曲線中高頻分量比較多。由圖知：當(dāng)N=8時(shí)，正反三角波的幅頻特性相同，因?yàn)閮烧叩臅r(shí)域只差一個(gè)相位；當(dāng)N=16時(shí)，正，反三角波的幅頻特性不同。這是因?yàn)闁艡谛?yīng)，當(dāng)N=8時(shí)，一些譜線被擋住。通過在原序列的末端補(bǔ)零，N=16，即增加采樣的點(diǎn)數(shù)和改變采樣的位置，使這些被擋住的譜線顯露出來，弱化了柵欄效應(yīng)。3 一個(gè)連續(xù)信號(hào)含兩個(gè)頻率分量，經(jīng)采樣得x(n)=sin2*0

10、.125n+cos2*(0.125+f)n n=0,1,N-1已知N=16，f分別為1/16和1/64，觀察其頻譜；當(dāng)N=128時(shí)，f不變，其結(jié)果有何不同？ 程序：n=0:N-1;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+f)*n);subplot(2,1,1);plot(n,x);subplot(2,1,2);y=fft(x);y=abs(y);stem(n,y);N=16，f=1/16N=16,f=1/64N=128,f=1/16N=128,f=1/64結(jié)論：當(dāng)N=16，f=1/16，N=128,f=1/16以及N=128,f=1/64時(shí)，均反應(yīng)了真實(shí)的頻譜；

11、只有當(dāng)N=16,f=1/64時(shí)，頻譜發(fā)生了嚴(yán)重的柵欄效應(yīng)。這是由于分辨率等于1/N，當(dāng)f>=1/N時(shí)，能分辨，不會(huì)發(fā)生柵欄效應(yīng)；當(dāng)f=<1/N時(shí)，不能分辨，會(huì)發(fā)生柵欄效應(yīng)。4 用FFT 分別計(jì)算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16點(diǎn)循環(huá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);ya=fft(xa);ya=abs(ya);yb=fft(xb);yb=abs(yb);y1=ya.*yb;subplot(2,1,1

12、);stem(n,y1);yaa=fft(xa，32);yaa=abs(yaa);ybb=fft(xb，32);ybb=abs(ybb);y2=yaa.*ybb;subplot(2,1,2);n=0:31;stem(n,y2);（上圖是循環(huán)卷積，下圖是線性卷積）結(jié)論：比較圖中線性卷積與圓周卷積序列： Xa(n)（序列長(zhǎng)度為N1）與Xb(n)(序列長(zhǎng)度為N2)的N點(diǎn)圓周卷積序列（當(dāng)N<N1+N2-1）,即為將Xa(n)與Xb(n)線性卷積序列中序號(hào)從N到N1+N2-1的序列疊加到原序列序號(hào)從0到N-1的地方。5 產(chǎn)生一512點(diǎn)的隨即序列xe(n)并用xc(n)和 xe(n)做線形卷積，觀察

13、卷積前后xe(n) 頻譜的變化。要求將xe(n)分成8段，分別采用重疊相加法和重疊保留法。用重疊保留法和重疊相加法實(shí)現(xiàn)線形卷積的過程為： Xc(n)序列長(zhǎng)度為8，Xe(n)序列長(zhǎng)度為512，分Xe(n)序列為8段，每段長(zhǎng)度為64，則每段序列與Xc(n)序列卷積后的長(zhǎng)度為72，總長(zhǎng)度為520。（湊成2的整數(shù)倍） 程序：（重疊相加法）e=rand(1,512);n1=0:3;xc1=n1;n2=4:7;xc2=8-n2;xc=xc1,xc2;yc=fft(xc,72);/將短序列補(bǔ)零后做72點(diǎn)的FFTxe1=xe(1:64);ye1=fft(xe1,72);/對(duì)長(zhǎng)序列第一段做72點(diǎn)的FFTy1=y

14、e1.*yc;/將上述兩個(gè)FFT相乘y1=y1,zeros(1,448);/補(bǔ)上448個(gè)零，以便相加，以下7段重復(fù)上述過程xe2=xe(65:128);ye2=fft(xe2,72);y2=ye2.*yc;y2=zeros(1,64),y2,zeros(1,384);xe3=xe(129:192);ye3=fft(xe3,72);y3=ye3.*yc;y3=zeros(1,128),y3,zeros(1,320);xe4=xe(193:256);ye4=fft(xe4,72);y4=ye4.*yc;y4=zeros(1,192),y4,zeros(1,256);xe5=xe(257:320);

15、ye5=fft(xe5,72);y5=ye5.*yc;y5=zeros(1,256),y5,zeros(1,192);xe6=xe(321:384);ye6=fft(xe6,72);y6=ye6.*yc;y6=zeros(1,320),y6,zeros(1,128);xe7=xe(385:448);ye7=fft(xe7,72);y7=ye7.*yc;y7=zeros(1,384),y7,zeros(1,64);xe8=xe(449:512);ye8=fft(xe8,72);y8=ye8.*yc;y8=zeros(1,448),y8;y=y1+y2+y3+y4+y5+y6+y7+y8;/將這8

16、個(gè)序列相加，便可得到最終的結(jié)果。y=abs(y);n=1:520;plot(n,y)（重疊保留法）xe=rand(1,512);xe=zeros(1,8),xe,zeros(1,56)/長(zhǎng)序列前添8個(gè)零，后添56個(gè)零，構(gòu)成576點(diǎn)的序列n1=0:3; xc1=n1;n2=4:7;xc2=8-n2;xc=xc1,xc2;yc=fft(xc,72);對(duì)短序列做72點(diǎn)的FFTxe1=xe(1:72);/將所得序列分成8段，每段序列長(zhǎng)度為72ye1=fft(xe1,72);對(duì)長(zhǎng)序列的第一段做72點(diǎn)的FFTy1=ye1.*yc;將上述兩段序列相乘y1=y1(9:72);取第一段所得結(jié)果的后64點(diǎn)，以下七

17、段同上述布驟。xe2=xe(73:144);ye2=fft(xe2,72);y2=ye2.*yc;y2=y2(9:72);xe3=xe(145:216);ye3=fft(xe3,72);y3=ye3.*yc;y3=y3(9:72);xe4=xe(216:288);ye4=fft(xe4,72);y4=ye4.*yc;y4=y4(9:72);xe5=xe(289:360);ye5=fft(xe5,72);y5=ye5.*yc;y5=y5(9:72);xe6=xe(361:432);ye6=fft(xe6,72);y6=ye6.*yc;y6=y6(9:72);xe7=xe(433:504);ye7

18、=fft(xe7,72);y7=ye7.*yc;y7=y7(9:72);xe8=xe(505:576);ye8=fft(xe8,72);y8=ye8.*yc;y=y1,y2,y3,y4,y5,y6,y7,y8/將上述8段合并，變可得到最終結(jié)果y=abs(y);n=1:520;plot(n,y)結(jié)論：比較圖中序列的線形卷積頻譜：原序列的頻譜曲線較線性卷積序列的頻譜曲線陡峭，即一個(gè)長(zhǎng)序列與一個(gè)短序列作線性卷積，短序列就相當(dāng)于一個(gè)低通濾波器，濾除長(zhǎng)序列的一部分高頻分量；6 用FFT分別計(jì)算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16點(diǎn)循環(huán)相關(guān)和線形相關(guān)，問一共有多少

19、種結(jié)果，他們之間有何異同點(diǎn)。程序：function y=t27N=16;n=0:N-1;m=(-N+1):(N-1);xa=exp(-(n-8).2/2);xb=exp(-0.1*n).*sin(2*pi*0.0625*n);Xa1=abs(fft(xa,N);Xb1=abs(fft(xb,N);Xa2=fft(xa,2*N);Xb2=fft(xb,2*N);rm1=real(ifft(conj(Xa1).*Xb1);rm2=real(ifft(conj(Xa2).*Xb2);rm2=rm2(N+2:2*N) rm2(1:N);subplot(2,1,1)stem(n,rm1)subplot(2,1,2)stem(m,rm2)(上面是循環(huán)相關(guān)，下面是線性相關(guān))由上圖可以看到，16點(diǎn)的循環(huán)相關(guān)由于高斯噪聲的干擾，衰減發(fā)生了微小的變化，時(shí)間位置不對(duì)了。并且線形相關(guān)32點(diǎn)，循環(huán)相關(guān)只有16點(diǎn)。7 用FFT分別計(jì)算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.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);k1=length(xa);k2=length(x