離散傅里葉變換和快速傅里葉變換

.專(zhuān)業(yè)整理.實(shí)驗(yàn)報(bào)告課程名稱(chēng)： 信號(hào)分析與處理 指導(dǎo)老師： 成績(jī)：__________________實(shí)驗(yàn)名稱(chēng)：離散傅里葉變換和快速傅里葉變換 實(shí)驗(yàn)類(lèi)型： 基礎(chǔ)實(shí)驗(yàn) 同組學(xué)生姓名 ：第二次實(shí)驗(yàn) 離散傅里葉變換和快速傅里葉變換一、實(shí)驗(yàn)?zāi)康?.1掌握離散傅里葉變換 （DFT）的原理和實(shí)現(xiàn) ；1.2掌握快速傅里葉變換 （FFT）的原理和實(shí)現(xiàn)，掌握用FFT對(duì)連續(xù)信號(hào)和離散信號(hào)進(jìn)行譜分析的方法 。1.3會(huì)用Matlab 軟件進(jìn)行以上練習(xí) 。裝二、實(shí)驗(yàn)原理訂2.1關(guān)于DFT的相關(guān)知識(shí)線(xiàn)序列x(n)的離散事件傅里葉變換 （DTFT）表示為X(ej)x(n)ejn，n如果()為因果有限長(zhǎng)序列，n=0,1,...,N-1，則()的DTFT表示為xnxnN1X(ej)x(n)ejn，n0x(n)的離散傅里葉變換（DFT）表達(dá)式為N1j2nkX(k)x(n)eN(k0,1,...,N1)，n0.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.序列的N點(diǎn)DFT是序列DTFT在頻率區(qū)間[0,2π]上的N點(diǎn)燈間隔采樣，采樣間隔為2π/N。通過(guò)DFT，可以完成由一組有限個(gè)信號(hào)采樣值 x(n)直接計(jì)算得到一組有限個(gè)頻譜采樣值 X(k)。X(k)的幅度譜為X(k)22(k)，其中下標(biāo)R和I分別表示取實(shí)部、虛部的運(yùn)算。Xk的相位譜為XR(k)XI()(k)arctanXI(k)。XR(k)離散傅里葉反變換 （IDFT）定義為1N1j2nkx(n)NNn0X(k)e(n0,1,...,N1)。2.2關(guān)于FFT的相關(guān)知識(shí)j2 n快速傅里葉變換 （FFT）是DFT的快速算法，并不是一個(gè)新的映射 。FFT利用了e N 函數(shù)的周期性和對(duì)稱(chēng)性以及一些特殊值來(lái)減少 DFT的運(yùn)算量，可使DFT的運(yùn)算量下降幾個(gè)數(shù)量級(jí) ，從而使數(shù)字信號(hào)處理的速度大大提高 。若信號(hào)是連續(xù)信號(hào) ，用FFT進(jìn)行譜分析時(shí) ，首先必須對(duì)信號(hào)進(jìn)行采樣 ，使之變成離散信號(hào) ，然后就可以用 FFT來(lái)對(duì)連續(xù)信號(hào)進(jìn)行譜分析 。為了滿(mǎn)足采樣定理 ，一般在采樣之前要設(shè)置一個(gè)抗混疊低通濾波器，且抗混疊濾波器的截止頻率不得高于與采樣頻率的一半 。比較DFT和IDFT的定義，兩者的區(qū)別僅在于指數(shù)因子的指數(shù)部分的符號(hào)差異和幅度尺度變換 ，因此可用FFT算法來(lái)計(jì)算 IDFT。三、實(shí)驗(yàn)內(nèi)容與相關(guān)分析 （共6道）說(shuō)明：為了便于老師查看 ，現(xiàn)將各題的內(nèi)容寫(xiě)在這里 ——題目按照 3.1、3.2、...、3.6排列。每道題包含如下內(nèi)容 ：題干、解答（思路、M文件源代碼、命令.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.窗口中的運(yùn)行及其結(jié)果 ）、分析。其中“命令窗口中的運(yùn)行及其結(jié)果 ”按照小題順序排列 ，各小題包含命令與結(jié)果（圖形或者序列）。3.1 求有限長(zhǎng)離散時(shí)間信號(hào) x(n)的離散時(shí)間傅里葉變換（DTFT）X(ejΩ)并繪圖。．．12n20n10。（1）已知x(n)；（2）已知x(n)2n0其他【解答】思路：這是DTFT的變換，按照定義編寫(xiě) DTFT的M文件即可。考慮到自變量Ω是連續(xù)的，為了方便計(jì)算機(jī)計(jì)算，計(jì)算時(shí)只取三個(gè)周期 [-2π,4π]中均勻的 1000個(gè)點(diǎn)用于繪圖 。理論計(jì)算的各序列 DTFT表達(dá)式，請(qǐng)見(jiàn)本題的分析 。M文件源代碼（我的Matlab 源文件不支持中文注釋 ，抱歉）：function DTFT(n1,n2,x)%ThisisaDTFTfunctionformyexperimentofSignalProcessing&Analysis.w=0:2*pi/1000:2*pi; %Definethebracketofomegaforplotting.X=zeros(size(w));%DefinetheinitialvaluesofX.fori=n1:n2X=X+x(i-n1+1)*exp((-1)*j*w*i); %ItisthedefinitionofDTFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);plot(w,Amp,'r');xlabel('\Omega' );ylabel('Amplification' );hold on;%Plotamplificationontheleft..學(xué)習(xí)幫手..專(zhuān)業(yè)整理.subplot(1,2,2);plot(w,Phs,'b');xlabel('\Omega' );ylabel('PhaseAngle(radian)' );holdoff;%Plotphaseangleontheright.end命令窗口中的運(yùn)行及其結(jié)果 （理論計(jì)算的各序列 DTFT表達(dá)式，請(qǐng)見(jiàn)本題的分析 ）：第（1）小題n=(-2:2);x=1.^n;DTFT(-2,2,x);544.53423.5)na1n3doat(aeci2.5l0fgnlpAmeA2sa-1hP1.5-210.5-300510-40510-5-5圖在[-2π,4π]范圍內(nèi)3個(gè)周期的幅度譜和相位譜（弧度制）第（2）小題n=(0:10);x=2.^n;DTFT(0,10,x);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.220042000318002)1600n1andoat(aegf14000nlpAmeAs1200a-1hP1000-2800-36000510-40510-5-5圖在[-2π,4π]范圍內(nèi)3個(gè)周期的幅度譜和相位譜（弧度制）【分析】對(duì)于第（1）小題，由于序列 x(n)只在有限區(qū)間（-2，-1，-，1，2）上為1，所以是離散非周期的信號(hào)。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號(hào) 。事實(shí)上，我們可計(jì)算出它的表達(dá)式 ：22j1e5j1ejnjne

5jX()x(n)eejX()nn21e1e

j

，可見(jiàn)幅度頻譜擁有主極大和次極大，兩個(gè)主極大間有 |5-1|=4個(gè)極小，即有3個(gè)次級(jí)大。而對(duì)于它的相位頻譜 ，則是周期性地在π、0、π之間震蕩。對(duì)于第（2）小題，由于是離散非周期的信號(hào)。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號(hào)。而它的表10n1111j11達(dá)式：X()x(n)ejn2ej12eX()2，因此主極大之間只有jnn012e12ej|0-1|=1個(gè)極小，不存在次級(jí)大。而對(duì)于它的相位頻譜，則是在一個(gè)長(zhǎng)為2π的周期內(nèi)有|11-1|=10次振蕩。而由DTFT的定義可知，頻譜都是以 2π為周期向兩邊無(wú)限延伸的 。由于DTFT是連續(xù)譜，對(duì)于計(jì)算機(jī)處理來(lái)說(shuō)特別困難 ，因此我們才需要離散信號(hào)的頻譜也離散 ，由此構(gòu)造出 DFT（以及為加速計(jì)算 DFT的.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.FFT）。3.2已知有限長(zhǎng)序列 x(n)={8,7,9,5,1,7,9,5},試分別采用 DFT和FFT求其離散傅里葉變換 X(k)的幅度、相位圖?！窘獯稹克悸罚喊凑斩x編寫(xiě) M文件即可。M文件源代碼：i)DFT函數(shù)：function DFT(N,x)%ThisisaDFTfunctionformyexperimentofSignalProcessing&Analysis.k=(0:N-1);%DefinevariablekforDFT.X=zeros(size(k));%DefinetheinitialvalvesofX.fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i); %ItisthedefinitionofDFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);stem(k,Amp,'.',’MarkerSize’,18);xlabel('k');ylabel('Amplification' );hold on;%Plotamplificationontheleft.subplot(1,2,2);stem(k,Phs,'*');xlabel('k');ylabel('PhaseAngle(radian)' );holdoff;.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.%Plotphaseangleontheright.endii)基2-FFT函數(shù)function myFFT(N,x)%Thisisabase-2FFTfunction.lov=(0:N-1);j1=0;fori=1:N %indexedaddressingifi<j1+1temp=x( j1+1);x(j1+1)=x(i);x(i)=temp;endk=N/2;whilek<=j1j1=j1-k;k=k/2;endj1=j1+k;enddigit=0;k=N;.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.whilek>1digit=digit+1;k=k/2;endn=N/2;%Nowwestartthe"butterfly-shaped"process.formu=1:digitdif=2^(mu-1); %Differncebetweentheindexesofthetargetvariables.idx=1;fori=1:nidx1=idx;idx2=1;forj1=1:N/(2*n)r=(idx2-1)*2^(digit-mu);wn=exp(j*(-2)*pi*r/N); %Itisthe"circulatingcoefficients".temp=x(idx);x(idx)=temp+x(idx+dif)*wn;x(idx+dif)=temp-x(idx+dif)*wn;idx=idx+1;idx2=idx2+1;endidx=idx1+2*dif;end.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.n=n/2;endAmp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.',’MarkerSize’,18);xlabel('FFTk');ylabel('FFTAmplification' );hold on;%Plottheamplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)' );holdoff;end命令窗口中的運(yùn)行及其結(jié)果 ：DFT：x=[8,7,9,5,1,7,9,5];DFT(8,x);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.60 350 240)1nandoat(ae30g0fnlpAmeAsah20P-110 -202468-3246800kk圖的幅度譜和相位譜（弧度制）FFT：x=[8,7,9,5,1,7,9,5];myFFT(8,x);603502)n40a1ndoat(aecginl300pAmeAsTahFPF20T-1FF10-20-30246802468FFTkFFTk圖3.2.2FFT算法的幅度譜和相位譜（弧度制）【分析】DFT是離散信號(hào)、離散頻譜之間的映射 。在這里我們可以看到序列的頻譜也被離散化 。事實(shí)上，我們可以循著 DFT構(gòu)造的方法驗(yàn)證這個(gè)頻譜 ：.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.首先，將序列做 N=8周期延拓，成為離散周期信號(hào) 。然后利用 DFS計(jì)算得到延拓后的頻譜 ：N127jknjknX()x(n)e8x(n)e4n0n0

51,7,9 4i,7,3,6,9 4i,7 n 0,1,...,7，從而取 DFS以上作N 8延拓的主值區(qū)間得到 DFT，與圖一致。因此計(jì)算正確。而對(duì)于FFT，我們可以看到它給出和 DFT一樣的結(jié)果，說(shuō)明了FFT算法就是 DFT的一個(gè)等價(jià)形式 。不過(guò)，由于序列不夠長(zhǎng) ，F(xiàn)FT在計(jì)算速度上的優(yōu)越性尚未凸顯 。3.3已知連續(xù)時(shí)間信號(hào)x(t)=3cos8πt,X(ω)=3[(8)(8)]，該信號(hào)從t=0開(kāi)始以采樣周期Ts=0.1s進(jìn)行采樣得到序列x(n)，試選擇合適的采樣點(diǎn)數(shù)，分別采用DFT和FFT求其離散傅里葉變換X(k)的幅度、相位圖，并將結(jié)果與X(k)的幅度、相位圖，并將結(jié)果與X(ω)相比較?！窘獯稹克悸罚捍祟}與下一題都是一樣的操作，可以在編程時(shí)統(tǒng)一用變量g（0或1）來(lái)控制是否有白噪聲。這里取g=0（無(wú)白噪聲）。另外，分別取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長(zhǎng)度的選擇與頻譜是否泄漏的關(guān)系。M文件源代碼：i)采樣函數(shù)：function xs=sampJune3(N,Ts,g)%ThisisafunctionappliedtoProblem3&4.%N:numberofsamplingpoints;Ts:samplingperiod;g=0:Nogaussiannoise;g=1:gussiannoiseexists.n=1:N;fori=1:N%Notethatimuststartat0inanalysis.ThusImadeaadaptation.sample(i)=3*cos(8*pi*Ts*(i-1))+g*randn;.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.%InMatlab,indexstartsat1,whichisnotconsistentwithourhabit.ThusImadeaadaptationincodes.%Itisasamplingprocesswith(g=1)/without(g=0)noise.endxs=sample;plot(n-1,sample, '.',’MarkerSize’,18);xlabel('n');ylabel('value');hold off;Mustuse(n-1),becauseinMatlab,indexstartsat1.endii)DFT和基2-FFT函數(shù)的代碼，請(qǐng)見(jiàn)第3.2節(jié)。不需再新編一個(gè)。命令窗口中的運(yùn)行及其結(jié)果 ：點(diǎn)采樣：>>xs=sampJune3(12,0.1,0);% 末尾的0表示無(wú)噪聲。DFT(12,xs);myFFT(12,xs);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.321eula 0v-1-2-30 2 4 6 8 10 12n圖進(jìn)行12點(diǎn)采樣得到的序列202.5182161.514)1nan12d0.5aort(aecf10g0nlpAmeA8s-0.5ah6P-14-1.52-20-2.5024681012024681012kk圖3.3.2DFT的幅度譜和相位譜（弧度制），出現(xiàn)了泄漏.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.202.5182161.514)1nanidoai12r0.5t(aecgiin100pAmeAsTa8h-0.5FFPT6F-1F4-1.52-20-2.5024681012024681012FFTkFFTk圖3.3.3FFT的幅度譜和相位譜（弧度制）。出現(xiàn)了頻譜泄漏。點(diǎn)采樣：>>xs=sampJune3(20,0.1,0);% 末尾的0表示無(wú)噪聲。DFT(20,xs);myFFT(20,xs);321eula 0v-1-2-30 2 4 6 8 10 12 14 16 18 20n圖進(jìn)行20點(diǎn)采樣得到的序列.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.354303252)nandoari201aeclgfnlpA0m15eAsahP10-15-205101520-3510152000kk圖3.3.5DFT的幅度譜和相位譜（弧度制）。頻譜無(wú)泄漏。35430325)2nandoa1irt20(aeignl0pAmeA15sTahFP-1FTF10F-25-305101520-4510152000FFTkFFTk圖3.3.6FFT的幅度譜和相位譜（弧度制）。頻譜無(wú)泄漏。28點(diǎn)采樣：>>xs=sampJune3(28,0.1,0);% 末尾的0表示無(wú)噪聲。DFT(28,xs);myFFT(28,xs);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.321eul0av-1-2-3051015202530n圖3.3.7進(jìn)行28點(diǎn)采樣得到的序列404353302)n25anidoair1t(aecl20gfnpA0meAs15ahP10-15-20102030-310203000kk圖的幅度譜和相位譜（弧度制） 。再次出現(xiàn)頻譜泄漏。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.40435330)2nanido25airt(ae1cligin20pAmeAsTa0hF15PFTFF-110-250102030-310203000FFTkFFTk圖的幅度譜和相位譜（弧度制） 。再次出現(xiàn)頻譜泄漏?！痉治觥糠謩e取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長(zhǎng)度的選擇與頻譜是否泄漏之間的關(guān)系 。現(xiàn)在與原信號(hào)頻譜X( )比較后可以得出如下結(jié)論 ：圖原信號(hào)的頻譜（由兩個(gè)沖激函數(shù)組成）原信號(hào)的頻譜是X()3[(8)(8)]，在±8π上各有一強(qiáng)度為3π的譜線(xiàn)，在其余頻率上為0?？梢?jiàn)原信號(hào)被0.1s采樣周期的采樣信號(hào)離散化之后，譜線(xiàn)以20π為周期重復(fù)，并且只在（20k±8）π(k為整數(shù)）處非0。那么，在20點(diǎn)DFT（采樣時(shí)間原信號(hào)周期的整數(shù)倍）中，只有第8根、第12根譜.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.線(xiàn)非0。而在12點(diǎn)、28點(diǎn)DFT中，由于采樣時(shí)間不是原信號(hào)周期的整數(shù)倍 ，譜線(xiàn)將向兩邊泄漏 。不過(guò)，對(duì)比12點(diǎn)采樣和 28點(diǎn)采樣，我們還可以發(fā)現(xiàn) ，28點(diǎn)采樣頻譜的主譜線(xiàn)高度是次譜線(xiàn)高度的4倍，兒12點(diǎn)采樣頻譜的主譜線(xiàn)高度是次譜線(xiàn)高度的 3倍?？梢?jiàn)，在無(wú)法保證采樣時(shí)間是信號(hào)周期整數(shù)倍的情況下，增加采樣時(shí)間有助于減輕頻譜泄漏的程度 。3.4對(duì)第3步中所述連續(xù)時(shí)間信號(hào)疊加高斯白噪聲信號(hào) ，重復(fù)第3步過(guò)程。【解答】思路：此題與上一題都是一樣的操作 ，可以在編程時(shí)統(tǒng)一用變量 g（0或1）來(lái)控制是否有白噪聲 。這里取g=1（有白噪聲）。另外，仍然分別取 12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長(zhǎng)度的選擇與頻譜是否泄漏的關(guān)系 。M文件源代碼：不需要再新編程序 。可以直接引用上面的函數(shù) ：sampJune3(N,Ts,g)，取g=1，以體現(xiàn)存在白噪聲DFT(N,x)myFFT(N,x)命令窗口中的運(yùn)行及其結(jié)果 ：點(diǎn)采樣：>>xs=sampJune3(12,0.1,1);% 末尾的1表示有噪聲。DFT(12,xs);myFFT(12,xs);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.notiacfiilpmA

eulav1816141210864200 2

543210-1-2-3-4-50 2 4 6 8 10 12n圖進(jìn)行12點(diǎn)采樣得到的含噪聲的序列32)1nadar(eg0nAesah-1-24681012-3246810120kk圖含噪聲序列 DFT的幅度譜和相位譜（弧度制） 。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.18316214)nan12i1doat(aec10giin0pAmeA8sTahFPF6T-1FF4-220-3024681012024681012FFTkFFTk圖3.4.3含噪聲FFT的幅度譜和相位譜（弧度制）。點(diǎn)采樣：>>xs=sampJune3(20,0.1,1);% 末尾的1表示有噪聲。DFT(20,xs);myFFT(20,xs);4321eula 0v-1-2-3-40 2 4 6 8 10 12 14 16 18 20n圖進(jìn)行20點(diǎn)采樣得到的含噪聲序列.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.30325220)1nanidoar(tae15g0fnlpAmeAsah10P-15-20-30510152005101520kk圖3.4.5含噪聲DFT的幅度譜和相位譜（弧度制）。303252)nan201doairt(aecligfn150pAmeAsTahFPF10T-1FF5-20-30510152005101520FFTkFFTk圖3.4.6含噪聲FFT的幅度譜和相位譜（弧度制）。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.點(diǎn)采樣：>>xs=sampJune3(28,0.1,0);% 末尾的1表示有噪聲。DFT(28,xs);myFFT(28,xs);4321eula 0v-1-2-3-40 5 10 15 20 25 30n圖進(jìn)行28點(diǎn)采樣得到的含噪聲序列.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.403530n 25oitacifil 20pmA1510500403530noi 25tacfiilp 20mATFF 1510500【分析】

321)naida 0(elgnA-1sahP-2-3102030-41020300kk圖含噪聲DFT的幅度譜和相位譜（弧度制）。43)2naidar(e1lgnAesa0hPTFF-1-2102030-31020300FFTkFFTk圖3.4.9含噪聲FFT的幅度譜和相位譜（弧度制）。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.依然分別取 12點(diǎn)、20點(diǎn)、28點(diǎn)采樣。仍然與原信號(hào)的頻譜 X( ) 3[ ( 8) ( 8)]（圖3.3.10）比較，可以得到結(jié)論：由于疊加了噪聲，所以頻譜都受到了一定的干擾。由于白噪聲在各個(gè)頻率的功率相等，因此頻譜上各處的干擾也是均勻隨機(jī)的。不過(guò)，通過(guò)對(duì)比我們可以發(fā)現(xiàn)，20點(diǎn)采樣（無(wú)噪聲時(shí)不發(fā)生泄漏的采樣方法）在存在噪聲時(shí)，仍然可以明顯區(qū)分出原信號(hào)的譜線(xiàn)。第二好的是28點(diǎn)采樣，因?yàn)椴蓸訒r(shí)間較長(zhǎng)，即使存在頻譜泄漏也能較好地區(qū)分原信號(hào)的譜線(xiàn)。而最差的是12點(diǎn)采樣，由于噪聲的存在和嚴(yán)重的頻譜泄漏，它的次譜線(xiàn)與主譜線(xiàn)的高度相差不大，使原信號(hào)不明顯。j2k3.5已知序列x(n)4(n)3(n1)2(n2)(n3)，X(k)是x(n)的6點(diǎn)DFT，設(shè)WNkeN。（1）若有限長(zhǎng)序列y(n)的6點(diǎn)DFT是Y(k)4kX(k)，求yn。W6()（2）若有限長(zhǎng)序列w(n)的6點(diǎn)DFTW(k)是X(k)的實(shí)部，求w(n)。（3）若有限長(zhǎng)序列q(n)的3點(diǎn)DFT是Q(k)X(2k)，k=0,1,2，求q(n)?！窘獯稹克悸罚哼@是對(duì)DFT進(jìn)行變換后求IDFT?？紤]到IDFT和DFT定義的對(duì)稱(chēng)性，可以在DFT的基礎(chǔ)上略加調(diào)整既可用于計(jì)算。首先，∵x(n)4(n)3(n1)2(n2)(n3)，∴它的6點(diǎn)采樣是序列是x(n){4,3,2,1,0,0}。值得指出的是，在Matlab中，數(shù)組的序號(hào)是從1開(kāi)始的（而在信號(hào)分析中習(xí)慣從0開(kāi)始），不過(guò)我在上面編程時(shí)已考慮到這一情況，具體可見(jiàn)實(shí)驗(yàn)報(bào)告最后的“附錄”。首先生成x(n)的6點(diǎn)DFT，再按照各小題分別轉(zhuǎn)換 ，最后求相應(yīng)的 IDFT。M文件源代碼：輸出x(n)的6點(diǎn)DFT的函數(shù)：.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.function X=ExportDFT(N,x)%ThisisaDFTfunctionthatexportsthesolutiontovectorX.k=(0:N-1); %DefinevariablekforDFT.X=zeros(size(k)); %DefinetheinitialvalvesofX.fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i); %ItisthedefinitionofDFT.endendii)算第（1）小題的Y(k)的函數(shù)：function Y=Convertor1(X)%Thisisamathematicalconvertorforthesubproblem(1).fork=1:6Y(k)=exp((-j)*2*pi*(-4*(k-1))/6)*X(k);%Herewemustuse(k-1),becauseinMatlabindexstartsat1.endend算第（2）小題的W(k)的函數(shù)：function W=Convertor2(X)%Thisisamathematicalconvertorforthesubproblem(2).W=real(X);%AcquiretherealpartofX.endiv)算第（3）小題的Q(k)的函數(shù)：.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.function Q=Convertor3(X)%Thisisamathematicalconvertorforthesubproblem(3).Q=zeros(3);% Detailedexplanationgoesherefortmp=1:3Q(tmp)=X(2*tmp);endendv)輸出IDFT的函數(shù)：function x=ExportIDFT(N,X)%ThisistheIDFTfunctionformyexperiment.n=(0:N-1);%DefinevariablenforIDFT.x=zeros(size(n));%Definetheinitialvalvesofx.fork=0:N-1x=x+X(k+1)*exp(j*2*k*pi/N*n);endx=x/N;a=real(x);%WeMUSTusereal(x),thoughwemayALREADYknowxisreal.%It'scausedbynumericcalculation(notanalyticcalculation)inMatlab.stem(n,a,'.','MarkerSize',18);xlabel('n');ylabel('Solution');end命令窗口中的運(yùn)行及其結(jié)果 ：>>x=[4,3,2,1,0,0];.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.>>X=ExportDFT(6,x);第（1）小題Y=Convertor1(X);y=ExportIDFT(6,Y)y=Columns1through40.0000+0.0000i 0.0000+0.0000i 4.0000+0.0000i 3.0000+0.0000i%虛部都是 0，說(shuō)明是實(shí)數(shù)Columns5through62.0000+0.0000i 1.0000-0.0000i %虛部都是 0，說(shuō)明是實(shí)數(shù)%事實(shí)上，在Matlab 中，由于數(shù)值計(jì)算的 截?cái)嗾`差，對(duì)原復(fù)數(shù)做乘法后 ，答案的虛部可能有一極小的量。43.532.5n2ouo1.5S10.50-0.50 1 2 3 4 5n

答案：y(n)={0,0,4,3,2,1}圖輸出的y(n)，這是對(duì)x(n)的圓周移位。第（2）小題W=Convertor2(X);w=ExportIDFT(6,W)w=.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.Columns1through44.0000+0.0000iColumns5through61.0000+0.0000i%事實(shí)上，在Matlab量。43.532.5notiul 2oS1.510.500 1 2

1.5000+0.0000i 1.0000+0.0000i 1.0000+0.0000i%虛部都是 0，說(shuō)明是實(shí)數(shù)1.5000+0.0000i% 虛部都是0，說(shuō)明是實(shí)數(shù)；中，由于數(shù)值計(jì)算的截?cái)嗾`差 ，對(duì)原復(fù)數(shù)做乘法后 ，答案的虛部可能有一極小的答案：w(n)={0,0,4,3,2,1}圖輸出的w(n)。3 4 5n第（3）小題Q=Convertor3(X);q=ExportIDFT(6,Q)q=Columns1through41.5000-0.0000i-0.1667-0.2887i0.7500-1.2990i0.8333-0.0000iColumns5through6-0.5000-0.8660i 1.0833-1.8764i這里的答案都是幅值 、相位均非 0的復(fù)數(shù)，而教材（實(shí)驗(yàn)指導(dǎo)第 109頁(yè)）并未要求作圖，這里略去。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.答案：q(n)={1.5，-0.1667-0.2887i，0.7500-1.2990i，0.8333，-0.5000-0.8660i，1.0833-1.8764i}【分析】對(duì)原序列進(jìn)行DFT運(yùn)算后，可以得到X(k)={10，3.5-4.33i，2.5-0.866i，2，2.5+0.866i，3.5+4.330i}。第（1）小題，根據(jù)DFT的性質(zhì)可以判斷，對(duì)原頻譜乘上旋轉(zhuǎn)因子W64k之后進(jìn)行IDFT得到的y(n)，就是對(duì)原序列做圓周移位：y(n)15X(k)ejk0n0,0,4,3,2,1。Nk0第（2）小題，由于對(duì)原頻譜取了實(shí)部，那么根據(jù)DFT的奇偶虛實(shí)性知，得到的w(n)也是實(shí)數(shù)的。第（3）小題，對(duì)原信號(hào)進(jìn)行了尺度變換（“抽取”），導(dǎo)致丟失了一些譜線(xiàn)，使得無(wú)法通過(guò)IDFT得到原來(lái)的序列()。說(shuō)明頻譜記錄了原有信號(hào)的信息，若頻譜發(fā)生變化，則對(duì)應(yīng)的時(shí)域信號(hào)也隨之改變。xn3.6已知信號(hào)x(t) sin(2f1t) sin(2f2t) sin(2f3t)，其中f1=4Hz、f2=4.02Hz、f3=5Hz，采用采樣頻率為20Hz進(jìn)行采樣，求1）當(dāng)采樣長(zhǎng)度N分別為512和2048情況下x(t)的幅度頻譜；．．．．2）當(dāng)采樣長(zhǎng)度N為32，且增補(bǔ)N個(gè)零點(diǎn)、4N個(gè)零點(diǎn)、8N個(gè)零點(diǎn)、16N個(gè)零點(diǎn)情況下x(t)的幅度頻．．．譜。．【解答】思路：采樣是有限且離散的 ，用DFT（FFT算法）計(jì)算頻譜，以便得到離散的頻譜 ，并且具有較高速度 。20Hz對(duì)應(yīng)的采樣周期 Ts=0.05s。M文件源代碼：i)采樣函數(shù)（其中Plus表示采樣后補(bǔ)零的個(gè)數(shù) ）function xs=sampJune6(N,Plus).學(xué)習(xí)幫手..專(zhuān)業(yè)整理.%ThisisafunctionappliedtoProblem6%N:samlingpoints;Plus:quantityofadditinalzeropoints.Ts=1/20;w1=2*pi*4;w2=2*pi*4.02;w3=2*pi*5;sample=zeros(1,N+Plus);n=(1:N);sample=sin(w1*Ts*(n-1))+sin(w2*Ts*(n-1))+sin(w3*Ts*(n-1));forp=(N+1):(N+Plus)sample(p)=0;%Addzeropoints.endxs=sample;%Returnendii)由于只要求顯示幅度頻譜 ，所以刪去 myFFT函數(shù)中繪制相位頻譜的命令 ，使它的最后部分修改如下 ：原來(lái)的：function myFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification' );hold on;%Plottheamplification..學(xué)習(xí)幫手..專(zhuān)業(yè)整理.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)' );holdoff;end修改后的：function myFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification' );%Plottheamplification.end命令窗口中的運(yùn)行及其結(jié)果 ：第（1）小題x512=sampJune6(512,0);x2048=sampJune6(2048,0);myFFT(512,x512);myFFT(2048,x2048);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.512Points300250200ntacfiilmAFF1005000 50 100 150 200 250 300 350 400 450 500FFTk圖進(jìn)行512點(diǎn)采樣得到的頻譜2048Points12001000800noc600pATF40020000200400600800100012001400160018002000FFTk圖進(jìn)行2048點(diǎn)采樣得到的頻譜第（2）小題>>x32p1N=sampJune6(32,32*1);%32點(diǎn)采樣，補(bǔ)零N=32個(gè)，共64個(gè)數(shù)據(jù)點(diǎn)>>x32p4N=sampJune6(32,32*4);%32點(diǎn)采樣，補(bǔ)零4N=128個(gè)，共160個(gè)數(shù)據(jù)點(diǎn)>>x32p8N=sampJune6(32,32*8);%32點(diǎn)采樣，補(bǔ)零8N=256個(gè)，共288個(gè)數(shù)據(jù)點(diǎn)>>x32p16N=sampJune6(32,32*16);%32點(diǎn)采樣，補(bǔ)零16N=128個(gè)，共544個(gè)數(shù)據(jù)點(diǎn)>>myFFT(32+32*1,x32p1N);>>myFFT(32+32*4,x32p4N);.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.myFFT(32+32*8,x32p8N);myFFT(32+32*16,x32p16N);353025n tafilpmAF151050010203040506070FFTk圖采樣N=32點(diǎn)，補(bǔ)零N=32點(diǎn)，共64點(diǎn)的頻譜35302520afilpmATFF

151050020406080100120140160FFTk圖采樣N=32點(diǎn)，補(bǔ)零 4N=128點(diǎn)，共160點(diǎn)的頻譜.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.353025n taficilpmATFF 151050050100150200250300FFTk圖采樣N=32點(diǎn)，補(bǔ)零 8N=32點(diǎn)，共288點(diǎn)的頻譜353025ni 20fiilpmATFF 10500 100 200 300 400 500 600FFTk圖采樣N=32點(diǎn)，補(bǔ)零16N=32點(diǎn)，共544點(diǎn)的頻譜【分析】請(qǐng)注意，題目只要求繪制幅度頻譜 。第（1）小題：首先，由于采樣時(shí)間都不是原有信號(hào)周期的整數(shù)倍 ，兩個(gè)采樣方式對(duì)應(yīng)的頻譜均發(fā)生了泄漏 。不過(guò)由于2048點(diǎn)采樣對(duì)應(yīng)的采樣時(shí)間較長(zhǎng) ，它頻譜泄漏的程度比 512點(diǎn)采樣輕。其次，由于20Hz的2048點(diǎn)采樣的頻率分辨率為 20/2048=0.0098Hz<0.2Hz ，因此放大頻譜圖之.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.后我們可以看到 4Hz、4.02Hz和5Hz對(duì)應(yīng)的譜線(xiàn)，而512點(diǎn)采樣的頻率分辨率為 20/512=0.039Hz>0.2Hz，因此4Hz和4.02Hz對(duì)應(yīng)的譜線(xiàn)無(wú)法區(qū)分 。第（2）小題：首先，由于采樣時(shí)間都不是原有信號(hào)周期的整數(shù)倍 ，頻譜均發(fā)生了泄漏 。而且由于采樣時(shí)間較短 ，頻譜泄漏比第（1）小題的兩個(gè)頻譜更加嚴(yán)重 。其次，由于都是 32點(diǎn)采樣，因此實(shí)際的頻率分辨率較低 ，無(wú)法區(qū)分 4Hz和4.02Hz對(duì)應(yīng)的譜線(xiàn)。最后，雖然都是 32點(diǎn)采樣，但由于補(bǔ) 0個(gè)數(shù)的不同，各頻譜譜線(xiàn)間距各不相同 。例如，補(bǔ)零最多的序列一共有 544個(gè)數(shù)據(jù)點(diǎn)，因此譜線(xiàn)間距小 。由此還可以得出結(jié)論 ：數(shù)據(jù)點(diǎn)個(gè)數(shù)越多 ，則頻譜越傾向于連續(xù)。可見(jiàn)，當(dāng)采樣時(shí)間不是原信號(hào)周期整數(shù)倍而且采樣時(shí)間較短時(shí) ，頻譜泄漏相當(dāng)嚴(yán)重的 ，所有的頻率上都有了幅值即能量 ，可見(jiàn)當(dāng)取樣信號(hào)的樣點(diǎn)數(shù)取得不夠時(shí) ，原信號(hào)所攜帶的信息就不能被完全取得 。而若將取樣信號(hào)補(bǔ)零 ，由圖可見(jiàn)信號(hào)的能量相應(yīng)的泄漏到了幾乎所有頻率上了 ，這樣所得的信號(hào)仍然嚴(yán)重失真，因此不能靠將信號(hào)補(bǔ)零這樣的方法來(lái)取得更精確的采樣信號(hào) 。要想獲得不泄漏的頻譜 ，在采樣頻率不變的前提下 ，必須使采樣時(shí)間等于原信號(hào)周期的整數(shù)倍 ，或者盡量延長(zhǎng)采樣時(shí)間以減少泄漏 。四、實(shí)驗(yàn)體會(huì)4.1關(guān)于各個(gè)實(shí)驗(yàn)的分析，請(qǐng)見(jiàn)第3部分每道題的末尾。4.2在Matlab中，數(shù)組的序號(hào)是從1開(kāi)始的，這與信號(hào)處理時(shí)通常的序號(hào)起點(diǎn)（0）不一致。我在編程充分注意到了這個(gè)問(wèn)題。4.3由于Matlab進(jìn)行數(shù)值計(jì)算的過(guò)程中存在截?cái)嗾`差，所以最后算得的值并不是準(zhǔn)確值。例如，對(duì)一個(gè)復(fù)數(shù)z，即使f(z)的虛部為零，但由于截?cái)嗾`差的存在（特別是z的虛部為無(wú)窮小數(shù)時(shí)），最終f(z)值的虛部可能是一個(gè)極小的非零值，從而在顯示時(shí)出現(xiàn)“零虛部”（例如，2+0.0000i)。.學(xué)習(xí)幫手..專(zhuān)業(yè)整理.4.4通過(guò)利用軟件對(duì)離散信號(hào)的各種變換 DTFT、DFT以及其快速算法 FFT進(jìn)行計(jì)算，使得在實(shí)驗(yàn)中比較難以實(shí)現(xiàn)的信號(hào)分析過(guò)程 （離散信號(hào)的采集和顯示都是比較困難的 ）在計(jì)算機(jī)計(jì)算中實(shí)現(xiàn) ，證明了理論的正確性，說(shuō)明仿真計(jì)算是一種十分有效的輔助手段 。4.5通過(guò)這次實(shí)驗(yàn)和上次實(shí)驗(yàn) 《信號(hào)的采集與恢復(fù) 》我知道了，要想盡量不失真地取得一個(gè)信號(hào)的頻譜 （低混疊、低泄漏），應(yīng)該盡量滿(mǎn)足以下條件 ：1）使用的開(kāi)關(guān)函數(shù)要盡量接近理想沖激串；2）采樣頻率要高于原始信號(hào)的奈奎斯特頻率。對(duì)于頻譜不受限的信號(hào)，為了避免頻譜混疊，應(yīng)該使用低通濾波器進(jìn)行濾波；3）對(duì)于頻帶不受限的信號(hào)，抗混疊濾波器要盡量接近理想濾波器。4）采樣的持續(xù)時(shí)間最好能夠是原信號(hào)周期的整數(shù)倍，一避免頻譜泄漏。而當(dāng)不知道原信號(hào)的周期（或者周期不穩(wěn)定）時(shí)，就要通過(guò)延長(zhǎng)采樣時(shí)間來(lái)盡量減少泄漏，從而突出原信號(hào)的譜線(xiàn)。5）當(dāng)信號(hào)混有白噪聲時(shí)，就更應(yīng)注意減少頻譜的泄漏和混疊，否則信號(hào)分析更加困難，甚至可能會(huì)使原信號(hào)被誤差“淹沒(méi)”。6）若原信號(hào)有多個(gè)頻率成分，應(yīng)該盡量提高采樣的頻率分辨率，以區(qū)分出更細(xì)微的頻率差異。4.6在實(shí)驗(yàn)中，在計(jì)算2048點(diǎn)采樣時(shí)，初步體會(huì)到了 FFT算法的優(yōu)越性 。在我的計(jì)算機(jī)上 ，F(xiàn)FT算法的確比原始的DFT更快。不過(guò)由于采樣點(diǎn)數(shù)較少 ，這一差別僅限于幾秒鐘 。在采樣點(diǎn)更多時(shí) ，F(xiàn)FT在速度上的優(yōu)越性應(yīng)該能進(jìn)一步突出 。4.7實(shí)驗(yàn)中遇到的問(wèn)題及其解決 ：實(shí)驗(yàn)中有些 M 文件代碼總是出錯(cuò) 。解決方