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

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

j

，可見幅度頻譜擁有主極大和次極大，兩個主極大間有 |5-1|=4個極小，即有3個次級大。而對于它的相位頻譜 ，則是周期性地在π、0、π之間震蕩。對于第（2）小題，由于是離散非周期的信號。它的幅度頻譜相應地應該是周期連續信號。而它的表10n1111j11達式：X()x(n)ejn2ej12eX()2，因此主極大之間只有jnn012e12ej|0-1|=1個極小，不存在次級大。而對于它的相位頻譜，則是在一個長為2π的周期內有|11-1|=10次振蕩。而由DTFT的定義可知，頻譜都是以 2π為周期向兩邊無限延伸的 。由于DTFT是連續譜，對于計算機處理來說特別困難 ，因此我們才需要離散信號的頻譜也離散 ，由此構造出 DFT（以及為加速計算 DFT的.學習幫手..專業整理.FFT）。3.2已知有限長序列 x(n)={8,7,9,5,1,7,9,5},試分別采用 DFT和FFT求其離散傅里葉變換 X(k)的幅度、相位圖。【解答】思路：按照定義編寫 M文件即可。M文件源代碼：i)DFT函數：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;.學習幫手..專業整理.%Plotphaseangleontheright.endii)基2-FFT函數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;.學習幫手..專業整理.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.學習幫手..專業整理.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命令窗口中的運行及其結果 ：DFT：x=[8,7,9,5,1,7,9,5];DFT(8,x);.學習幫手..專業整理.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是離散信號、離散頻譜之間的映射 。在這里我們可以看到序列的頻譜也被離散化 。事實上，我們可以循著 DFT構造的方法驗證這個頻譜 ：.學習幫手..專業整理.首先，將序列做 N=8周期延拓，成為離散周期信號 。然后利用 DFS計算得到延拓后的頻譜 ：N127jknjknX()x(n)e8x(n)e4n0n0

51,7,9 4i,7,3,6,9 4i,7 n 0,1,...,7，從而取 DFS以上作N 8延拓的主值區間得到 DFT，與圖一致。因此計算正確。而對于FFT，我們可以看到它給出和 DFT一樣的結果，說明了FFT算法就是 DFT的一個等價形式 。不過，由于序列不夠長 ，FFT在計算速度上的優越性尚未凸顯 。3.3已知連續時間信號x(t)=3cos8πt,X(ω)=3[(8)(8)]，該信號從t=0開始以采樣周期Ts=0.1s進行采樣得到序列x(n)，試選擇合適的采樣點數，分別采用DFT和FFT求其離散傅里葉變換X(k)的幅度、相位圖，并將結果與X(k)的幅度、相位圖，并將結果與X(ω)相比較。【解答】思路：此題與下一題都是一樣的操作，可以在編程時統一用變量g（0或1）來控制是否有白噪聲。這里取g=0（無白噪聲）。另外，分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏的關系。M文件源代碼：i)采樣函數：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;.學習幫手..專業整理.%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函數的代碼，請見第3.2節。不需再新編一個。命令窗口中的運行及其結果 ：點采樣：>>xs=sampJune3(12,0.1,0);% 末尾的0表示無噪聲。DFT(12,xs);myFFT(12,xs);.學習幫手..專業整理.321eula 0v-1-2-30 2 4 6 8 10 12n圖進行12點采樣得到的序列202.5182161.514)1nan12d0.5aort(aecf10g0nlpAmeA8s-0.5ah6P-14-1.52-20-2.5024681012024681012kk圖3.3.2DFT的幅度譜和相位譜（弧度制），出現了泄漏.學習幫手..專業整理.202.5182161.514)1nanidoai12r0.5t(aecgiin100pAmeAsTa8h-0.5FFPT6F-1F4-1.52-20-2.5024681012024681012FFTkFFTk圖3.3.3FFT的幅度譜和相位譜（弧度制）。出現了頻譜泄漏。點采樣：>>xs=sampJune3(20,0.1,0);% 末尾的0表示無噪聲。DFT(20,xs);myFFT(20,xs);321eula 0v-1-2-30 2 4 6 8 10 12 14 16 18 20n圖進行20點采樣得到的序列.學習幫手..專業整理.354303252)nandoari201aeclgfnlpA0m15eAsahP10-15-205101520-3510152000kk圖3.3.5DFT的幅度譜和相位譜（弧度制）。頻譜無泄漏。35430325)2nandoa1irt20(aeignl0pAmeA15sTahFP-1FTF10F-25-305101520-4510152000FFTkFFTk圖3.3.6FFT的幅度譜和相位譜（弧度制）。頻譜無泄漏。28點采樣：>>xs=sampJune3(28,0.1,0);% 末尾的0表示無噪聲。DFT(28,xs);myFFT(28,xs);.學習幫手..專業整理.321eul0av-1-2-3051015202530n圖3.3.7進行28點采樣得到的序列404353302)n25anidoair1t(aecl20gfnpA0meAs15ahP10-15-20102030-310203000kk圖的幅度譜和相位譜（弧度制） 。再次出現頻譜泄漏。.學習幫手..專業整理.40435330)2nanido25airt(ae1cligin20pAmeAsTa0hF15PFTFF-110-250102030-310203000FFTkFFTk圖的幅度譜和相位譜（弧度制） 。再次出現頻譜泄漏。【分析】分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏之間的關系 。現在與原信號頻譜X( )比較后可以得出如下結論 ：圖原信號的頻譜（由兩個沖激函數組成）原信號的頻譜是X()3[(8)(8)]，在±8π上各有一強度為3π的譜線，在其余頻率上為0。可見原信號被0.1s采樣周期的采樣信號離散化之后，譜線以20π為周期重復，并且只在（20k±8）π(k為整數）處非0。那么，在20點DFT（采樣時間原信號周期的整數倍）中，只有第8根、第12根譜.學習幫手..專業整理.線非0。而在12點、28點DFT中，由于采樣時間不是原信號周期的整數倍 ，譜線將向兩邊泄漏 。不過，對比12點采樣和 28點采樣，我們還可以發現 ，28點采樣頻譜的主譜線高度是次譜線高度的4倍，兒12點采樣頻譜的主譜線高度是次譜線高度的 3倍。可見，在無法保證采樣時間是信號周期整數倍的情況下，增加采樣時間有助于減輕頻譜泄漏的程度 。3.4對第3步中所述連續時間信號疊加高斯白噪聲信號 ，重復第3步過程。【解答】思路：此題與上一題都是一樣的操作 ，可以在編程時統一用變量 g（0或1）來控制是否有白噪聲 。這里取g=1（有白噪聲）。另外，仍然分別取 12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏的關系 。M文件源代碼：不需要再新編程序 。可以直接引用上面的函數 ：sampJune3(N,Ts,g)，取g=1，以體現存在白噪聲DFT(N,x)myFFT(N,x)命令窗口中的運行及其結果 ：點采樣：>>xs=sampJune3(12,0.1,1);% 末尾的1表示有噪聲。DFT(12,xs);myFFT(12,xs);.學習幫手..專業整理.notiacfiilpmA

eulav1816141210864200 2

543210-1-2-3-4-50 2 4 6 8 10 12n圖進行12點采樣得到的含噪聲的序列32)1nadar(eg0nAesah-1-24681012-3246810120kk圖含噪聲序列 DFT的幅度譜和相位譜（弧度制） 。.學習幫手..專業整理.18316214)nan12i1doat(aec10giin0pAmeA8sTahFPF6T-1FF4-220-3024681012024681012FFTkFFTk圖3.4.3含噪聲FFT的幅度譜和相位譜（弧度制）。點采樣：>>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圖進行20點采樣得到的含噪聲序列.學習幫手..專業整理.30325220)1nanidoar(tae15g0fnlpAmeAsah10P-15-20-30510152005101520kk圖3.4.5含噪聲DFT的幅度譜和相位譜（弧度制）。303252)nan201doairt(aecligfn150pAmeAsTahFPF10T-1FF5-20-30510152005101520FFTkFFTk圖3.4.6含噪聲FFT的幅度譜和相位譜（弧度制）。.學習幫手..專業整理.點采樣：>>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圖進行28點采樣得到的含噪聲序列.學習幫手..專業整理.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的幅度譜和相位譜（弧度制）。.學習幫手..專業整理.依然分別取 12點、20點、28點采樣。仍然與原信號的頻譜 X( ) 3[ ( 8) ( 8)]（圖3.3.10）比較，可以得到結論：由于疊加了噪聲，所以頻譜都受到了一定的干擾。由于白噪聲在各個頻率的功率相等，因此頻譜上各處的干擾也是均勻隨機的。不過，通過對比我們可以發現，20點采樣（無噪聲時不發生泄漏的采樣方法）在存在噪聲時，仍然可以明顯區分出原信號的譜線。第二好的是28點采樣，因為采樣時間較長，即使存在頻譜泄漏也能較好地區分原信號的譜線。而最差的是12點采樣，由于噪聲的存在和嚴重的頻譜泄漏，它的次譜線與主譜線的高度相差不大，使原信號不明顯。j2k3.5已知序列x(n)4(n)3(n1)2(n2)(n3)，X(k)是x(n)的6點DFT，設WNkeN。（1）若有限長序列y(n)的6點DFT是Y(k)4kX(k)，求yn。W6()（2）若有限長序列w(n)的6點DFTW(k)是X(k)的實部，求w(n)。（3）若有限長序列q(n)的3點DFT是Q(k)X(2k)，k=0,1,2，求q(n)。【解答】思路：這是對DFT進行變換后求IDFT。考慮到IDFT和DFT定義的對稱性，可以在DFT的基礎上略加調整既可用于計算。首先，∵x(n)4(n)3(n1)2(n2)(n3)，∴它的6點采樣是序列是x(n){4,3,2,1,0,0}。值得指出的是，在Matlab中，數組的序號是從1開始的（而在信號分析中習慣從0開始），不過我在上面編程時已考慮到這一情況，具體可見實驗報告最后的“附錄”。首先生成x(n)的6點DFT，再按照各小題分別轉換 ，最后求相應的 IDFT。M文件源代碼：輸出x(n)的6點DFT的函數：.學習幫手..專業整理.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)的函數：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)的函數：function W=Convertor2(X)%Thisisamathematicalconvertorforthesubproblem(2).W=real(X);%AcquiretherealpartofX.endiv)算第（3）小題的Q(k)的函數：.學習幫手..專業整理.function Q=Convertor3(X)%Thisisamathematicalconvertorforthesubproblem(3).Q=zeros(3);% Detailedexplanationgoesherefortmp=1:3Q(tmp)=X(2*tmp);endendv)輸出IDFT的函數：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命令窗口中的運行及其結果 ：>>x=[4,3,2,1,0,0];.學習幫手..專業整理.>>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，說明是實數Columns5through62.0000+0.0000i 1.0000-0.0000i %虛部都是 0，說明是實數%事實上，在Matlab 中，由于數值計算的 截斷誤差，對原復數做乘法后 ，答案的虛部可能有一極小的量。43.532.5n2ouo1.5S10.50-0.50 1 2 3 4 5n

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

1.5000+0.0000i 1.0000+0.0000i 1.0000+0.0000i%虛部都是 0，說明是實數1.5000+0.0000i% 虛部都是0，說明是實數；中，由于數值計算的截斷誤差 ，對原復數做乘法后 ，答案的虛部可能有一極小的答案：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的復數，而教材（實驗指導第 109頁）并未要求作圖，這里略去。.學習幫手..專業整理.答案：q(n)={1.5，-0.1667-0.2887i，0.7500-1.2990i，0.8333，-0.5000-0.8660i，1.0833-1.8764i}【分析】對原序列進行DFT運算后，可以得到X(k)={10，3.5-4.33i，2.5-0.866i，2，2.5+0.866i，3.5+4.330i}。第（1）小題，根據DFT的性質可以判斷，對原頻譜乘上旋轉因子W64k之后進行IDFT得到的y(n)，就是對原序列做圓周移位：y(n)15X(k)ejk0n0,0,4,3,2,1。Nk0第（2）小題，由于對原頻譜取了實部，那么根據DFT的奇偶虛實性知，得到的w(n)也是實數的。第（3）小題，對原信號進行了尺度變換（“抽取”），導致丟失了一些譜線，使得無法通過IDFT得到原來的序列()。說明頻譜記錄了原有信號的信息，若頻譜發生變化，則對應的時域信號也隨之改變。xn3.6已知信號x(t) sin(2f1t) sin(2f2t) sin(2f3t)，其中f1=4Hz、f2=4.02Hz、f3=5Hz，采用采樣頻率為20Hz進行采樣，求1）當采樣長度N分別為512和2048情況下x(t)的幅度頻譜；．．．．2）當采樣長度N為32，且增補N個零點、4N個零點、8N個零點、16N個零點情況下x(t)的幅度頻．．．譜。．【解答】思路：采樣是有限且離散的 ，用DFT（FFT算法）計算頻譜，以便得到離散的頻譜 ，并且具有較高速度 。20Hz對應的采樣周期 Ts=0.05s。M文件源代碼：i)采樣函數（其中Plus表示采樣后補零的個數 ）function xs=sampJune6(N,Plus).學習幫手..專業整理.%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函數中繪制相位頻譜的命令 ，使它的最后部分修改如下 ：原來的：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..學習幫手..專業整理.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命令窗口中的運行及其結果 ：第（1）小題x512=sampJune6(512,0);x2048=sampJune6(2048,0);myFFT(512,x512);myFFT(2048,x2048);.學習幫手..專業整理.512Points300250200ntacfiilmAFF1005000 50 100 150 200 250 300 350 400 450 500FFTk圖進行512點采樣得到的頻譜2048Points12001000800noc600pATF40020000200400600800100012001400160018002000FFTk圖進行2048點采樣得到的頻譜第（2）小題>>x32p1N=sampJune6(32,32*1);%32點采樣，補零N=32個，共64個數據點>>x32p4N=sampJune6(32,32*4);%32點采樣，補零4N=128個，共160個數據點>>x32p8N=sampJune6(32,32*8);%32點采樣，補零8N=256個，共288個數據點>>x32p16N=sampJune6(32,32*16);%32點采樣，補零16N=128個，共544個數據點>>myFFT(32+32*1,x32p1N);>>myFFT(32+32*4,x32p4N);.學習幫手..專業整理.myFFT(32+32*8,x32p8N);myFFT(32+32*16,x32p16N);353025n tafilpmAF151050010203040506070FFTk圖采樣N=32點，補零N=32點，共64點的頻譜35302520afilpmATFF

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