數(shù)字信號處理實(shí)驗(yàn)快速傅里葉變換及其應(yīng)用

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

2、無限，一個(gè)時(shí)間無限的信號其頻帶寬度則為有限。因此對一個(gè)時(shí)間有限的信號，應(yīng)用DFT進(jìn)行分析，頻譜混疊難以避免。對一個(gè)時(shí)間無限的信號雖然頻帶有限，但在實(shí)際運(yùn)算中，時(shí)間總是取有限值，在將信號截?cái)嗟倪^程中，出現(xiàn)了分散的擴(kuò)展譜線的現(xiàn)象，稱之為頻譜泄露或功率泄露。（3） 柵欄效應(yīng)：DFT是對單位圓上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等份的同心圓的圓周上。兩圓上對應(yīng)的數(shù)兩量兩相乘求和，就得到全部卷積序列。這個(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序列Xa(n)=exp(-(n-p).2)/q),0=n=15 =0,其他 （b）衰減正弦序列 X(b)=exp(-an)*sin(2pi*fn),0=n=15 =0,其他 (c)三角波序列 Xb(n)=n,0=n=3 =8-n,4=n=7 =0,其他 (d)反三角波序列 Xc(n)=4-n,0=n=3 =n-4,4=n=7

4、 =0,其他上機(jī)實(shí)驗(yàn)內(nèi)容：1 觀察高斯序列的時(shí)域和幅頻特性，固定信號xa(n)中參數(shù)p=8，改變q的 值，使q分別等于2，4，8，觀察他們的時(shí)域和幅頻特性，了解當(dāng)q取不同值時(shí)，對信號序列的時(shí)域和幅頻特性影響；改變p，使p分別等于8，13，14，觀察參數(shù)p變化對信號序列的時(shí)域和幅頻特性影響，注意p等于多少時(shí)，會(huì)發(fā)生明顯的泄漏現(xiàn)象，混疊是否也隨之出現(xiàn)？記錄實(shí)驗(yàn)中觀察到的現(xiàn)象，繪出相應(yīng)的時(shí)域序列和幅頻特性曲線。代碼：%高斯序列以p，q為參數(shù)生成并輸出其時(shí)域序列波形以及其dtft；function timePlot,spectrumPlot = Gauss( p,q ) %輸入相關(guān)參數(shù)n=0:15;N

5、=4000;w=linspace(-pi,pi,N);g=exp(-(n-p).2/q);spec=dtft(N,g);subplot(2,1,1);stem(n,g);title(strcat(當(dāng)p=,num2str(p),q=,num2str(q),時(shí)的高斯時(shí)域波形);subplot(2,1,2);plot(w,spec);title(strcat(當(dāng)p=,num2str(p),q=,num2str(q),時(shí)的高斯頻域波形);end%DTFT變換得到頻譜：function spectrum = dtft( n,x )%取4000采樣點(diǎn)；w=linspace(-pi,pi,n);spectr

6、um=zeros(1,n);for k=1:n for m=1:length(x) wn(k,m)=exp(-1j*w(k)*m); end spectrum(k)=abs(x*wn(k,:);end end實(shí)驗(yàn)結(jié)果如下：.P=8,q=2P=8,q=4;P=8,q=8;P=8,q=8;P=13,q=8P=14,q=8;分析：當(dāng)p固定，q變化時(shí)，由時(shí)域波形可知，q代表了波形以p為中心向兩側(cè)衰減的速度，發(fā)現(xiàn)當(dāng)q越大時(shí)，衰減的速度越快，反之，衰減的速度相對較慢，高頻分量越小，從頻域波形分析可知，三種情況下的形狀是相似的，但是三種情況下的頻譜在q比較的大的時(shí)候在高頻上出現(xiàn)了較少的高頻分量以及較快的高頻

7、分量的衰減。當(dāng)q固定，p變化時(shí)，由時(shí)域波形可知，p代表了波形的移序，在頻譜上，顯而 易見，在p更大的時(shí)候產(chǎn)生了更多的高頻分量。右移之后，序列相當(dāng)于被截?cái)啵?那么相對應(yīng)的其頻譜泄漏的情況更加的嚴(yán)重，從而出現(xià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)象的原因。代碼：%按照alpha和f的值生成相應(yīng)的正弦序列，畫出其時(shí)域和頻域（dtft）的波形；function

8、 output_args =sinusoidal( alpha,f)%生成序列n=0:15;N=4000;x=exp(-alpha*n).*sin(2*pi*f*n);%對離散序列進(jìn)行dtft變換spec=dtft(N,x);subplot(2,1,1);stem(n,x);title(strcat(當(dāng)=,num2str(alpha),f=,num2str(f),時(shí)的正弦時(shí)域波形);subplot(2,1,2);w=linspace(-pi,pi,N);plot(w,spec);title(strcat(當(dāng)=,num2str(alpha),f=,num2str(f),時(shí)的正弦頻域波形); en

9、da=0.1,f=0.0625;F=0.4375F=0.5625;分析：在a=0.1固定的條件下，f=0.0625時(shí)，譜峰在低頻出約等于2*pi*f，但是在f=0.4375和f=0.5626的時(shí)候，譜峰主要集中在高頻位置。后兩種情況由于不滿足nyquist采樣定律，所以i出現(xiàn)了混疊。3.代碼;function = triangle( N )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor k=0:3 x(k+1)=k;endfor k=4:7 x(k+1)=8-k;endspec=f

10、ft(x);spec=abs(spec);n=0:N-1;subplot(1,2,1);stem(n,x);title(序列波形);subplot(1,2,2);stem(n,spec);title(strcat(num2str(N),點(diǎn)FFT);End逆三角波： function = reversetriangle( N )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor k=0:3 x(k+1)=4-k;endfor k=4:7 x(k+1)=k-4;endfor k=8:N-1

11、 x(k+1)=0;endn=0:N-1;spec=fft(x);spec=abs(spec);subplot(2,1,1);stem(n,x);title(序列波形);subplot(2,1,2);stem(n,spec);title(strcat(num2str(N),點(diǎn)FFT); end正三角波逆三角波：N=32;逆三角波正三角波分析：正三角波和逆三角波的時(shí)域波形上的走勢是不一樣的。在頻域上可以發(fā)現(xiàn)，N=8時(shí)，兩者的fft互為相反數(shù)，因?yàn)閮烧叩暮蛯?shí)際上就是一個(gè)矩形窗函數(shù)，rectwin(8)，由于rectwin（8）的fft為8 0 0 0 0 0 0 0，所以除了fft（0）以外，其他

12、都為相反數(shù)。而當(dāng)N=16時(shí)，兩者之間的和為一個(gè)補(bǔ)零了的rectwin（8），對于N不同時(shí)得到的結(jié)果發(fā)現(xiàn)，實(shí)質(zhì)上同一個(gè)序列補(bǔ)零以后的頻譜并沒有發(fā)生太大的變化，補(bǔ)零后的序列相當(dāng)了有了更多的采樣點(diǎn)，在對應(yīng)位置處的fft值還是相同的，不同的是，當(dāng)N變大時(shí)，補(bǔ)零長度變大時(shí)，暴露出來的更多的采樣點(diǎn)值可以減輕fft的柵欄效應(yīng)。4、一個(gè)連續(xù)信號含兩個(gè)頻率分量，經(jīng)采樣得x(n)=sin2*0.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é)果有何不同？代碼：function output_args = test4( d

13、elta_f,n)%UNTITLED7 Summary of this function goes here% Detailed explanation goes herex=zeros(1,n);for k=0:n-1 x(k+1)=sin(2*pi*0.125*k)+cos(2*pi*(0.125+delta_f)*k);end;spec=fft(x);stem(spec);title(strcat(deltaf=,num2str(delta_f),N=,num2str(n); end N=16，deltaf=1/16;N=16，deltaf=1/64;N=128,deltaf=1/16N

14、=128,deltaf=1/64分析：由fft頻譜圖可知，在N=16時(shí)，deltaf=1/16時(shí)，fft頻譜是正確的應(yīng)該是兩個(gè)沖激，但是在deltaf=1/64時(shí)，此時(shí)的信號的帶寬減小，由于采樣點(diǎn)數(shù)是固定的，導(dǎo)致卷積的過程中出現(xiàn)了泄露的情況，因此在第二種情況下出現(xiàn)了譜線擴(kuò)散的情況。在N=128時(shí)，相當(dāng)于擴(kuò)展了加窗的寬度，因而沒有出現(xiàn)頻譜的泄漏。5、 用FFT 分別計(jì)算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16點(diǎn)循環(huán)卷積和線形卷積。代碼：n=0:15;alpha=0.1;f=0.0625;p=8;q=2;k=0:31;g=exp(-(n-p).2/q);x=e

15、xp(-alpha*n).*sin(2*pi*f*n);for i=17:32 g(i)=0; x(i)=0;endspecg=fft(g);specx=fft(x);convSpec=specg.*specx;linearConv=ifft(convSpec);stem(k,linearConv);循環(huán)卷積：代碼：n=0:15;alpha=0.1;f=0.0625;p=8;q=2;x=exp(-alpha*n).*sin(2*pi*f*n); g=exp(-(n-p).2/q);specg=fft(g);specx=fft(x);shiftSpec=specg.*specx;shiftCon

16、v=ifft(shiftSpec);stem(n,shiftConv);6、代碼：clear all;clc;N=512;kn=1:N;a=1:2*N-1;%生成512個(gè)均值為1，方差為1的隨機(jī)序列；x = 1+1.*randn(512,1);x=x;for k=1:4 t(k)=k;endfor k=5:8 t(k)=8-k;end%直接線性卷積convout=conv(x,t);specx=fft(t);subplot(3,1,1);stem(convout);title(直接線性卷積);out=zeros(1,519);out1=zeros(1,519);%重疊相加法%重疊相加法for

17、m=9:71 t(m)=0;endspec=fft(t);%重疊相加法outm=zeros(1,519);for k=1:8 over(k,:)=x(64*(k-1)+1):64*k) zeros(1,7); spectrumn(k,:)=fft(over(k,:).*spec; spectrumn(k,:)=ifft(spectrumn(k,:); changespectrum(k,:)=zeros(1,(k-1)*64) spectrumn(k,:) zeros(1,519-64*k-7); outm=outm+changespectrum(k,:);end subplot(3,1,2);

18、stem(outm);title(重疊相加法線性卷積); %重疊保留法xx=zeros(1,7),x;for k=1:8 over1(k,:)=xx(64*(k-1)+1:64*(k-1)+71);end over1(9,:)=xx(513:519) zeros(1,64);for k=1:9 speco1(k,:)=fft(over1(k,:); spectrumout1(k,:)=speco1(k,:).*spec; spectrumout1(k,:)=ifft(spectrumout1(k,:);end s=spectrumout1; out1=s(1,8:end) s(2,8:end)

19、 s(3,8:end) s(4,8:end) s(5,8:end) s(6,8:end) s(7,8:end) s(8,8:end) s(9,8:14); subplot(3,1,3); stem(out1);title(重疊保留法線性卷積); 7、線性相關(guān)：線性相關(guān)即為互相關(guān)，在此不再重述。如題8；循環(huán)相關(guān)：clear;p=8;q=2;n=0:15;xa=exp(-(n-p).2)/q);a=0.1;f=0.0625;xb=exp(-a*n).*sin(2*pi*f*n);xak=fft(xa);xbk=fft(xb);subplot(2,1,1);yn=ifft(conj(xak).*xb

20、k);stem(n,yn);xlabel(n);ylabel(幅度);title(循環(huán)相關(guān));subplot(2,1,2);yn=ifft(conj(xbk).*xak);stem(n,yn);xlabel(n);ylabel(幅度);title(循環(huán)相關(guān));8. 互相關(guān)代碼：rxyclear;clc;n=16;alpha=0.1;f=0.0625;p=8;q=2;for i =1:n g(i)=exp(-(i-p)2/q); x(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1);endfor i=n+1:2*n g(i)=0; x(i)=0;endspecx=fft(

21、x,2*n);specg=fft(g,2*n);corr=ifft(conj(specx).*specg);rm=real(corr);rm=rm(n+2:2*n) rm(1:n);m=(-n+1):(n-1);stem(m,rm);結(jié)果：Ryxclear;clc;n=16;alpha=0.1;f=0.0625;p=8;q=2;for i =1:n g(i)=exp(-(i-p)2/q); x(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1);endfor i=n+1:2*n g(i)=0; x(i)=0;endspecx=fft(x,2*n);specg=fft(g,2*n);corr=ifft(specx.*conj(specg);rm=real(corr);rm=rm(n+2:2*n) rm(1:n);m=(-n+1):(n-1);stem(m,rm);自相關(guān)、xa（n）clear;clc;n=16;alpha=0.1;f=0.0625;p=8;q=2;for i =1:n g(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1); x(i)=exp(-alpha*(i-1)*sin(2*pi*f*(i-1);endfor i=n+1:2*n g(i)=0; x(i)=0;end