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

1、實(shí)驗(yàn)報(bào)告課程名稱： 信號分析與處理 指導(dǎo)老師： 成績：_實(shí)驗(yàn)名稱：離散傅里葉變換和快速傅里葉變換 實(shí)驗(yàn)類型： 基礎(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對連續(xù)信號和離散信號進(jìn)行譜分析的方法。1.3 會(huì)用Matlab軟件進(jìn)行以上練習(xí)。二、實(shí)驗(yàn)原理2.1關(guān)于DFT的相關(guān)知識序列x(n)的離散事件傅里葉變換（DTFT）表示為，如果x(n)為因果有限長序列，n=0,1,.,N-1，則x(n)的DTFT表示為，x(n)的離散傅里葉變換（DFT）表達(dá)式

2、為，序列的N點(diǎn)DFT是序列DTFT在頻率區(qū)間0,2上的N點(diǎn)燈間隔采樣，采樣間隔為2/N。通過DFT，可以完成由一組有限個(gè)信號采樣值x(n)直接計(jì)算得到一組有限個(gè)頻譜采樣值X(k)。X(k)的幅度譜為，其中下標(biāo)R和I分別表示取實(shí)部、虛部的運(yùn)算。X(k)的相位譜為。離散傅里葉反變換（IDFT）定義為。 2.2關(guān)于FFT的相關(guān)知識快速傅里葉變換（FFT）是DFT的快速算法，并不是一個(gè)新的映射。FFT利用了函數(shù)的周期性和對稱性以及一些特殊值來減少DFT的運(yùn)算量，可使DFT的運(yùn)算量下降幾個(gè)數(shù)量級，從而使數(shù)字信號處理的速度大大提高。若信號是連續(xù)信號，用FFT進(jìn)行譜分析時(shí)，首先必須對信號進(jìn)行采樣，使之變成離

3、散信號，然后就可以用FFT來對連續(xù)信號進(jìn)行譜分析。為了滿足采樣定理，一般在采樣之前要設(shè)置一個(gè)抗混疊低通濾波器，且抗混疊濾波器的截止頻率不得高于與采樣頻率的一半。 比較DFT和IDFT的定義，兩者的區(qū)別僅在于指數(shù)因子的指數(shù)部分的符號差異和幅度尺度變換，因此可用FFT算法來計(jì)算IDFT。3、 實(shí)驗(yàn)內(nèi)容與相關(guān)分析（共6道）說明：為了便于老師查看，現(xiàn)將各題的內(nèi)容寫在這里題目按照3.1、3.2、.、3.6排列。每道題包含如下內(nèi)容：題干、解答（思路、M文件源代碼、命令窗口中的運(yùn)行及其結(jié)果）、分析。其中“命令窗口中的運(yùn)行及其結(jié)果”按照小題順序排列，各小題包含命令與結(jié)果（圖形或者序列）。3.1 求有限長離散時(shí)

4、間信號x(n)的離散時(shí)間傅里葉變換（DTFT）X(ej)并繪圖。（1） 已知；（2）已知?！窘獯稹克悸罚哼@是DTFT的變換，按照定義編寫DTFT的M文件即可?？紤]到自變量是連續(xù)的，為了方便計(jì)算機(jī)計(jì)算，計(jì)算時(shí)只取三個(gè)周期-2,4中均勻的1000個(gè)點(diǎn)用于繪圖。理論計(jì)算的各序列DTFT表達(dá)式，請見本題的分析。M文件源代碼（我的Matlab源文件不支持中文注釋，抱歉）：function DTFT(n1,n2,x)%This is a DTFT function for my experiment of Signal Processing & Analysis.w=0:2*pi/1000:2*p

5、i;%Define the bracket of omega for plotting.X=zeros(size(w);%Define the initial values of X.for i=n1:n2 X=X+x(i-n1+1)*exp(-1)*j*w*i);%It is the definition of DTFT.endAmp=abs(X);%Acquire the amplification.Phs=angle(X);%Acquire the phase angle (radian).subplot(1,2,1);plot(w,Amp,'r'); xlabel(&#

6、39;Omega');ylabel('Amplification');hold on;%Plot amplification on the left.subplot(1,2,2);plot(w,Phs,'b');xlabel('Omega');ylabel('Phase Angle (radian)');hold off;%Plot phase angle on the right.end命令窗口中的運(yùn)行及其結(jié)果（理論計(jì)算的各序列DTFT表達(dá)式，請見本題的分析）：第（1）小題>> n=(-2:2);>&g

7、t; x=1.n;>>DTFT(-2,2,x);圖3.1.1在-2,4范圍內(nèi)3個(gè)周期的幅度譜和相位譜（弧度制）第（2）小題>> n=(0:10);>> x=2.n;>> DTFT(0,10,x);圖3.1.2在-2,4范圍內(nèi)3個(gè)周期的幅度譜和相位譜（弧度制）【分析】對于第（1）小題，由于序列x(n)只在有限區(qū)間（-2，-1，-，1，2）上為1，所以是離散非周期的信號。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號。事實(shí)上，我們可計(jì)算出它的表達(dá)式：，可見幅度頻譜擁有主極大和次極大，兩個(gè)主極大間有|5-1|=4個(gè)極小，即有3個(gè)次級大。而對于它的相位頻譜，則是周

8、期性地在-、0、之間震蕩。對于第（2）小題，由于是離散非周期的信號。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號。而它的表達(dá)式：，因此主極大之間只有|0-1|=1個(gè)極小，不存在次級大。而對于它的相位頻譜，則是在一個(gè)長為2的周期內(nèi)有|11-1|=10次振蕩。而由DTFT的定義可知，頻譜都是以2為周期向兩邊無限延伸的。由于DTFT是連續(xù)譜，對于計(jì)算機(jī)處理來說特別困難，因此我們才需要離散信號的頻譜也離散，由此構(gòu)造出DFT（以及為加速計(jì)算DFT的FFT）。3.2已知有限長序列x(n)=8,7,9,5,1,7,9,5,試分別采用DFT和FFT求其離散傅里葉變換X(k)的幅度、相位圖。【解答】思路：按照定義編寫M

9、文件即可。M文件源代碼：i) DFT函數(shù)：function DFT(N,x)%This is a DFT function for my experiment of Signal Processing & Analysis.k=(0:N-1);%Define variable k for DFT.X=zeros(size(k);%Define the initial valves of X.for i=0:N-1 X=X+x(i+1)*exp(-1)*j*2*k*pi/N*i);%It is the definition of DFT.endAmp=abs(X);%Acquire th

10、e amplification.Phs=angle(X);%Acquire the phase angle (radian).subplot(1,2,1);stem(k,Amp,'.',MarkerSize,18); xlabel('k');ylabel('Amplification');hold on;%Plot amplification on the left.subplot(1,2,2);stem(k,Phs,'*');xlabel('k');ylabel('Phase Angle (radian)

11、');hold off;%Plot phase angle on the right.endii) 基2-FFT函數(shù)function myFFT(N,x)%This is a base-2 FFT function.lov=(0:N-1);j1=0;for i=1:N %indexed addressing if i<j1+1 temp=x(j1+1); x(j1+1)=x(i); x(i)=temp; end k=N/2; while k<=j1 j1=j1-k; k=k/2; end j1=j1+k;enddigit=0;k=N;while k>1 digit=d

12、igit+1; k=k/2;endn=N/2;% Now we start the "butterfly-shaped" process.for mu=1:digit dif=2(mu-1);%Differnce between the indexes of the target variables. idx=1; for i=1:n idx1=idx; idx2=1; for j1=1:N/(2*n) r=(idx2-1)*2(digit-mu); wn=exp(j*(-2)*pi*r/N);%It is the "circulating coefficient

13、s". temp=x(idx); x(idx)=temp+x(idx+dif)*wn; x(idx+dif)=temp-x(idx+dif)*wn; idx=idx+1; idx2=idx2+1; end idx=idx1+2*dif; end n=n/2;endAmp=abs(x);%Acquire the amplification.Phs=angle(x);%Acquire the phase angle (radian).subplot(1,2,1);stem(lov,Amp,'.',MarkerSize,18);xlabel('FFT k')

14、;ylabel('FFT Amplification');hold on;%Plot the amplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFT k');ylabel('FFT Phase Angle (radian)');hold off;end命令窗口中的運(yùn)行及其結(jié)果：DFT：>> x=8,7,9,5,1,7,9,5;>> DFT(8,x);圖3.2.1 DFT的幅度譜和相位譜（弧度制）FFT：>> x=8,7,9,5,1,7

15、,9,5;>> myFFT(8,x);圖3.2.2 FFT算法的幅度譜和相位譜（弧度制）圖3.2.1 DFT的幅度譜和相位譜（相位是弧度制的）【分析】DFT是離散信號、離散頻譜之間的映射。在這里我們可以看到序列的頻譜也被離散化。事實(shí)上，我們可以循著DFT構(gòu)造的方法驗(yàn)證這個(gè)頻譜：首先，將序列做N=8周期延拓，成為離散周期信號。然后利用DFS計(jì)算得到延拓后的頻譜：，從而取DFS的主值區(qū)間得到DFT，與圖一致。因此計(jì)算正確。而對于FFT，我們可以看到它給出和DFT一樣的結(jié)果，說明了FFT算法就是DFT的一個(gè)等價(jià)形式。不過，由于序列不夠長，F(xiàn)FT在計(jì)算速度上的優(yōu)越性尚未凸顯。3.3已知連續(xù)

16、時(shí)間信號x(t)=3cos8t, X()=，該信號從t=0開始以采樣周期Ts=0.1 s進(jìn)行采樣得到序列x(n)，試選擇合適的采樣點(diǎn)數(shù)，分別采用DFT和 FFT求其離散傅里葉變換X(k)的幅度、相位圖，并將結(jié)果與X(k)的幅度、相位圖，并將結(jié)果與X()相比較?！窘獯稹克悸罚捍祟}與下一題都是一樣的操作，可以在編程時(shí)統(tǒng)一用變量g（0或1）來控制是否有白噪聲。這里取g=0（無白噪聲）。 另外，分別取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長度的選擇與頻譜是否泄漏的關(guān)系。M文件源代碼：i)采樣函數(shù)：function xs=sampJune3(N,Ts,g)%This is a function appl

17、ied to Problem 3 & 4.%N: number of sampling points; Ts: sampling period; g=0: No gaussian noise; g=1: gussian noise exists.n=1:N;for i=1:N%Note that i must start at 0 in analysis. Thus I made a adaptation.sample(i)=3*cos(8*pi*Ts*(i-1)+g*randn;%In Matlab, index starts at 1, which is not consisten

18、t with our habit. Thus I made a adaptation in codes. %It is a sampling process with(g=1)/without(g=0) noise.endxs=sample;plot(n-1,sample,'.',MarkerSize,18);xlabel('n');ylabel('value');hold off;% Must use (n-1), because in Matlab, index starts at 1.endii)DFT和基2-FFT函數(shù)的代碼，請見第3.2

19、節(jié)。不需再新編一個(gè)。命令窗口中的運(yùn)行及其結(jié)果：12點(diǎn)采樣：>> xs=sampJune3(12,0.1,0);%末尾的0表示無噪聲。>> DFT(12,xs);>> myFFT(12,xs);圖3.3.1 進(jìn)行12點(diǎn)采樣得到的序列圖3.3.2 DFT的幅度譜和相位譜（弧度制），出現(xiàn)了泄漏圖3.3.3 FFT的幅度譜和相位譜（弧度制）。出現(xiàn)了頻譜泄漏。20點(diǎn)采樣：>> xs=sampJune3(20,0.1,0);%末尾的0表示無噪聲。>> DFT(20,xs);>> myFFT(20,xs);圖3.3.4 進(jìn)行20點(diǎn)采樣得

20、到的序列圖3.3.5 DFT的幅度譜和相位譜（弧度制）。頻譜無泄漏。圖3.3.6 FFT的幅度譜和相位譜（弧度制）。頻譜無泄漏。 28點(diǎn)采樣：>> xs=sampJune3(28,0.1,0);%末尾的0表示無噪聲。>> DFT(28,xs);>> myFFT(28,xs);圖3.3.7 進(jìn)行28點(diǎn)采樣得到的序列圖3.3.8 DFT的幅度譜和相位譜（弧度制）。再次出現(xiàn)頻譜泄漏。圖3.3.9 FFT的幅度譜和相位譜（弧度制）。再次出現(xiàn)頻譜泄漏?！痉治觥?分別取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長度的選擇與頻譜是否泄漏之間的關(guān)系?，F(xiàn)在與原信號頻譜比較后可以得

21、出如下結(jié)論：圖3.3.10 原信號的頻譜（由兩個(gè)沖激函數(shù)組成）原信號的頻譜是，在±8上各有一強(qiáng)度為3的譜線，在其余頻率上為0?？梢娫盘柋?.1 s采樣周期的采樣信號離散化之后，譜線以20為周期重復(fù)，并且只在（20k±8） (k為整數(shù)）處非0。那么，在20點(diǎn)DFT（采樣時(shí)間原信號周期的整數(shù)倍）中，只有第8根、第12根譜線非0。而在12點(diǎn)、28點(diǎn)DFT中，由于采樣時(shí)間不是原信號周期的整數(shù)倍，譜線將向兩邊泄漏。 不過，對比12點(diǎn)采樣和28點(diǎn)采樣，我們還可以發(fā)現(xiàn)，28點(diǎn)采樣頻譜的主譜線高度是次譜線高度的4倍，兒12點(diǎn)采樣頻譜的主譜線高度是次譜線高度的3倍?？梢姡跓o法保證采樣時(shí)間

22、是信號周期整數(shù)倍的情況下，增加采樣時(shí)間有助于減輕頻譜泄漏的程度。3.4對第3步中所述連續(xù)時(shí)間信號疊加高斯白噪聲信號，重復(fù)第3步過程。【解答】思路：此題與上一題都是一樣的操作，可以在編程時(shí)統(tǒng)一用變量g（0或1）來控制是否有白噪聲。這里取g=1（有白噪聲）。 另外，仍然分別取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣，以考察采樣長度的選擇與頻譜是否泄漏的關(guān)系。M文件源代碼：不需要再新編程序。可以直接引用上面的函數(shù)：sampJune3(N,Ts,g)，取g=1，以體現(xiàn)存在白噪聲DFT(N,x)myFFT(N,x)命令窗口中的運(yùn)行及其結(jié)果：12點(diǎn)采樣：>> xs=sampJune3(12,0.1,1);%

23、末尾的1表示有噪聲。>> DFT(12,xs);>> myFFT(12,xs);圖3.4.1 進(jìn)行12點(diǎn)采樣得到的含噪聲的序列圖3.4.2 含噪聲序列DFT的幅度譜和相位譜（弧度制）。圖3.4.3 含噪聲FFT的幅度譜和相位譜（弧度制）。20點(diǎn)采樣：>> xs=sampJune3(20,0.1,1);%末尾的1表示有噪聲。>> DFT(20,xs);>> myFFT(20,xs);圖3.4.4 進(jìn)行20點(diǎn)采樣得到的含噪聲序列圖3.4.5 含噪聲DFT的幅度譜和相位譜（弧度制）。 圖3.4.6 含噪聲FFT的幅度譜和相位譜（弧度制）。2

24、8點(diǎn)采樣：>> xs=sampJune3(28,0.1,0);%末尾的1表示有噪聲。>> DFT(28,xs);>> myFFT(28,xs);圖3.4.7 進(jìn)行28點(diǎn)采樣得到的含噪聲序列圖3.4.8 含噪聲DFT的幅度譜和相位譜（弧度制）。圖3.4.9 含噪聲FFT的幅度譜和相位譜（弧度制）?！痉治觥恳廊环謩e取12點(diǎn)、20點(diǎn)、28點(diǎn)采樣。仍然與原信號的頻譜（圖3.3.10）比較，可以得到結(jié)論：由于疊加了噪聲，所以頻譜都受到了一定的干擾。由于白噪聲在各個(gè)頻率的功率相等，因此頻譜上各處的干擾也是均勻隨機(jī)的。不過，通過對比我們可以發(fā)現(xiàn)，20點(diǎn)采樣（無噪聲時(shí)不發(fā)生

25、泄漏的采樣方法）在存在噪聲時(shí)，仍然可以明顯區(qū)分出原信號的譜線。第二好的是28點(diǎn)采樣，因?yàn)椴蓸訒r(shí)間較長，即使存在頻譜泄漏也能較好地區(qū)分原信號的譜線。而最差的是12點(diǎn)采樣，由于噪聲的存在和嚴(yán)重的頻譜泄漏，它的次譜線與主譜線的高度相差不大，使原信號不明顯。3.5已知序列，X(k)是x(n)的6點(diǎn)DFT，設(shè)。（1） 若有限長序列y(n)的6點(diǎn)DFT是，求y(n)。（2） 若有限長序列w(n)的6點(diǎn)DFT W(k)是的實(shí)部，求w(n)。（3） 若有限長序列q(n)的3點(diǎn)DFT是，k=0,1,2，求q(n)?！窘獯稹克悸罚哼@是對DFT進(jìn)行變換后求IDFT?？紤]到IDFT和DFT定義的對稱性，可以在DFT的

26、基礎(chǔ)上略加調(diào)整既可用于計(jì)算。首先，它的6點(diǎn)采樣是序列是。值得指出的是，在Matlab中，數(shù)組的序號是從1開始的（而在信號分析中習(xí)慣從0開始），不過我在上面編程時(shí)已考慮到這一情況，具體可見實(shí)驗(yàn)報(bào)告最后的“附錄”。 首先生成x(n)的6點(diǎn)DFT，再按照各小題分別轉(zhuǎn)換，最后求相應(yīng)的IDFT。M文件源代碼：i) 輸出x(n)的6點(diǎn)DFT的函數(shù)：function X = ExportDFT(N,x)%This is a DFT function that exports the solution to vector X.k=(0:N-1); %Define variable k for DFT.X=ze

27、ros(size(k); %Define the initial valves of X.for i=0:N-1 X=X+x(i+1)*exp(-1)*j*2*k*pi/N*i); %It is the definition of DFT.endendii)算第（1）小題的Y(k)的函數(shù)：function Y = Convertor1(X)%This is a mathematical convertor for the subproblem (1).for k=1:6 Y(k)=exp(-j)*2*pi*(-4*(k-1)/6)*X(k);%Here we must use (k-1),be

28、cause in Matlab index starts at 1.endendiii)算第（2）小題的W(k)的函數(shù)：function W = Convertor2(X)%This is a mathematical convertor for the subproblem (2).W=real(X);%Acquire the real part of X.endiv)算第（3）小題的Q(k)的函數(shù)：function Q = Convertor3(X)%This is a mathematical convertor for the subproblem (3).Q=zeros(3);% D

29、etailed explanation goes herefor tmp=1:3 Q(tmp)=X(2*tmp);endendv)輸出IDFT的函數(shù)：function x = ExportIDFT(N,X)%This is the IDFT function for my experiment.n=(0:N-1);%Define variable n for IDFT.x=zeros(size(n);%Define the initial valves of x.for k=0:N-1 x=x+X(k+1)*exp(j*2*k*pi/N*n);endx=x/N;a=real(x);%We MU

30、ST use real(x), though we may ALREADY know x is real. %It's caused by numeric calculation (not analytic calculation) in Matlab.stem(n,a,'.','MarkerSize',18);xlabel('n');ylabel('Solution');end命令窗口中的運(yùn)行及其結(jié)果：>> x=4,3,2,1,0,0;>> X=ExportDFT(6,x);第（1）小題>&

31、gt; Y=Convertor1(X);>> y=ExportIDFT(6,Y)y = Columns 1 through 4 0.0000 + 0.0000i 0.0000 + 0.0000i 4.0000 + 0.0000i 3.0000 + 0.0000i%虛部都是0，說明是實(shí)數(shù) Columns 5 through 6 2.0000 + 0.0000i 1.0000 - 0.0000i %虛部都是0，說明是實(shí)數(shù)%事實(shí)上，在Matlab中，由于數(shù)值計(jì)算的截?cái)嗾`差，對原復(fù)數(shù)做乘法后，答案的虛部可能有一極小的量。答案：y(n)=0,0,4,3,2,1圖3.5.1 輸出的y(n)，這

32、是對x(n)的圓周移位。第（2）小題>> W=Convertor2(X);>> w=ExportIDFT(6,W)w = Columns 1 through 4 4.0000 + 0.0000i 1.5000 + 0.0000i 1.0000 + 0.0000i 1.0000 + 0.0000i%虛部都是0，說明是實(shí)數(shù) Columns 5 through 6 1.0000 + 0.0000i 1.5000 + 0.0000i %虛部都是0，說明是實(shí)數(shù)；%事實(shí)上，在Matlab中，由于數(shù)值計(jì)算的截?cái)嗾`差，對原復(fù)數(shù)做乘法后，答案的虛部可能有一極小的量。答案：w(n)=0,0

33、,4,3,2,1圖3.5.2 輸出的w(n)。第（3）小題>> Q=Convertor3(X);>> q=ExportIDFT(6,Q)q = Columns 1 through 4 1.5000 - 0.0000i -0.1667 - 0.2887i 0.7500 - 1.2990i 0.8333 - 0.0000i Columns 5 through 6 -0.5000 - 0.8660i 1.0833 - 1.8764i 這里的答案都是幅值、相位均非0的復(fù)數(shù)，而教材（實(shí)驗(yàn)指導(dǎo)第109頁）并未要求作圖，這里略去。 答案：q(n)=1.5，0.8333，1.0833-

34、1.8764i【分析】對原序列進(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ì)可以判斷，對原頻譜乘上旋轉(zhuǎn)因子之后進(jìn)行IDFT得到的y(n)，就是對原序列做圓周移位：。第（2）小題，由于對原頻譜取了實(shí)部，那么根據(jù)DFT的奇偶虛實(shí)性知，得到的w(n)也是實(shí)數(shù)的。第（3）小題，對原信號進(jìn)行了尺度變換（“抽取”），導(dǎo)致丟失了一些譜線，使得無法通過IDFT得到原來的序列x(n)。說明頻譜記錄了原有信號的信息，若頻譜發(fā)生變化，則對應(yīng)的時(shí)域信號也隨之改變。3.6已知信號，其中f1=4 Hz、f2

35、=4.02 Hz、f3=5 Hz，采用采樣頻率為20 Hz進(jìn)行采樣，求（1） 當(dāng)采樣長度N分別為512和2048情況下x(t)的幅度頻譜；（2） 當(dāng)采樣長度N為32，且增補(bǔ)N個(gè)零點(diǎn)、4N個(gè)零點(diǎn)、8N個(gè)零點(diǎn)、16N個(gè)零點(diǎn)情況下x(t)的幅度頻譜?！窘獯稹克悸罚翰蓸邮怯邢耷译x散的，用DFT（FFT算法）計(jì)算頻譜，以便得到離散的頻譜，并且具有較高速度。20 Hz對應(yīng)的采樣周期Ts=0.05s。M文件源代碼：i)采樣函數(shù)（其中Plus表示采樣后補(bǔ)零的個(gè)數(shù)）function xs=sampJune6(N,Plus)%This is a function applied to Problem 6%N:sa

36、mling points;Plus:quantity of additinal zero points.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);for p=(N+1):(N+Plus) sample(p)=0;%Add zero points.endxs=sample;%Returnendii)由于只要求顯示幅度頻譜，所以刪去myFFT函數(shù)中繪制相位頻譜的命令，使它的最后部分修改如下

37、：原來的：function myFFT(N,x)%This is a base-2 FFT function wrote by myself.Amp=abs(x);%Acquire the amplification.Phs=angle(x);%Acquire the phase angle (radian).subplot(1,2,1);stem(lov,Amp,'.');xlabel('FFT k');ylabel('FFT Amplification');hold on;%Plot the amplification.subplot(1,2

38、,2);stem(lov,Phs,'*');xlabel('FFT k');ylabel('FFT Phase Angle (radian)');hold off;end修改后的：function myFFT(N,x)%This is a base-2 FFT function wrote by myself.Amp=abs(x);%Acquire the amplification.stem(lov,Amp,'.');xlabel('FFT k');ylabel('FFT Amplification

39、9;);%Plot the amplification.end命令窗口中的運(yùn)行及其結(jié)果：第（1）小題>> x512=sampJune6(512,0);>> x2048=sampJune6(2048,0);>> myFFT(512,x512);>> myFFT(2048,x2048);圖3.6.1 進(jìn)行512點(diǎn)采樣得到的頻譜圖3.6.2 進(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,3

40、2*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);>> myFFT(32+32*8,x32p8N);>> myFFT(32+32*16,x32p16N);圖3.6.3 采樣N=32點(diǎn)，補(bǔ)零N=32

41、點(diǎn)，共64點(diǎn)的頻譜圖3.6.4 采樣N=32點(diǎn)，補(bǔ)零4N=128點(diǎn)，共160點(diǎn)的頻譜圖3.6.5 采樣N=32點(diǎn)，補(bǔ)零8N=32點(diǎn)，共288點(diǎn)的頻譜圖3.6.6 采樣N=32點(diǎn)，補(bǔ)零16N=32點(diǎn)，共544點(diǎn)的頻譜【分析】請注意，題目只要求繪制幅度頻譜。第（1）小題：首先，由于采樣時(shí)間都不是原有信號周期的整數(shù)倍，兩個(gè)采樣方式對應(yīng)的頻譜均發(fā)生了泄漏。不過由于2048點(diǎn)采樣對應(yīng)的采樣時(shí)間較長，它頻譜泄漏的程度比512點(diǎn)采樣輕。其次，由于20 Hz的2048點(diǎn)采樣的頻率分辨率為20/2048=0.0098 Hz < 0.2 Hz，因此放大頻譜圖之后我們可以看到4 Hz、4.02 Hz和5 Hz

42、對應(yīng)的譜線，而512點(diǎn)采樣的頻率分辨率為20/512=0.039 Hz > 0.2 Hz，因此4 Hz和4.02 Hz對應(yīng)的譜線無法區(qū)分。第（2）小題：首先，由于采樣時(shí)間都不是原有信號周期的整數(shù)倍，頻譜均發(fā)生了泄漏。而且由于采樣時(shí)間較短，頻譜泄漏比第（1）小題的兩個(gè)頻譜更加嚴(yán)重。其次，由于都是32點(diǎn)采樣，因此實(shí)際的頻率分辨率較低，無法區(qū)分4 Hz和4.02 Hz對應(yīng)的譜線。最后，雖然都是32點(diǎn)采樣，但由于補(bǔ)0個(gè)數(shù)的不同，各頻譜譜線間距各不相同。例如，補(bǔ)零最多的序列一共有544個(gè)數(shù)據(jù)點(diǎn)，因此譜線間距小。由此還可以得出結(jié)論：數(shù)據(jù)點(diǎn)個(gè)數(shù)越多，則頻譜越傾向于連續(xù)。可見，當(dāng)采樣時(shí)間不是原信號周期整

43、數(shù)倍而且采樣時(shí)間較短時(shí)，頻譜泄漏相當(dāng)嚴(yán)重的，所有的頻率上都有了幅值即能量，可見當(dāng)取樣信號的樣點(diǎn)數(shù)取得不夠時(shí)，原信號所攜帶的信息就不能被完全取得。而若將取樣信號補(bǔ)零，由圖可見信號的能量相應(yīng)的泄漏到了幾乎所有頻率上了，這樣所得的信號仍然嚴(yán)重失真，因此不能靠將信號補(bǔ)零這樣的方法來取得更精確的采樣信號。要想獲得不泄漏的頻譜，在采樣頻率不變的前提下，必須使采樣時(shí)間等于原信號周期的整數(shù)倍，或者盡量延長采樣時(shí)間以減少泄漏。四、實(shí)驗(yàn)體會(huì)4.1關(guān)于各個(gè)實(shí)驗(yàn)的分析，請見第3部分每道題的末尾。4.2在Matlab中，數(shù)組的序號是從1開始的，這與信號處理時(shí)通常的序號起點(diǎn)（0）不一致。我在編程充分注意到了這個(gè)問題。4.

44、3由于Matlab進(jìn)行數(shù)值計(jì)算的過程中存在截?cái)嗾`差，所以最后算得的值并不是準(zhǔn)確值。例如，對一個(gè)復(fù)數(shù)z，即使f(z)的虛部為零，但由于截?cái)嗾`差的存在（特別是z的虛部為無窮小數(shù)時(shí)），最終f(z)值的虛部可能是一個(gè)極小的非零值，從而在顯示時(shí)出現(xiàn)“零虛部”（例如，2+0.0000i )。4.4通過利用軟件對離散信號的各種變換DTFT、DFT以及其快速算法FFT進(jìn)行計(jì)算，使得在實(shí)驗(yàn)中比較難以實(shí)現(xiàn)的信號分析過程（離散信號的采集和顯示都是比較困難的）在計(jì)算機(jī)計(jì)算中實(shí)現(xiàn)，證明了理論的正確性，說明仿真計(jì)算是一種十分有效的輔助手段。4.5通過這次實(shí)驗(yàn)和上次實(shí)驗(yàn)信號的采集與恢復(fù)我知道了，要想盡量不失真地取得一個(gè)信號

45、的頻譜（低混疊、低泄漏），應(yīng)該盡量滿足以下條件：（1） 使用的開關(guān)函數(shù)要盡量接近理想沖激串；（2） 采樣頻率要高于原始信號的奈奎斯特頻率。對于頻譜不受限的信號，為了避免頻譜混疊，應(yīng)該使用低通濾波器進(jìn)行濾波；（3） 對于頻帶不受限的信號，抗混疊濾波器要盡量接近理想濾波器。（4） 采樣的持續(xù)時(shí)間最好能夠是原信號周期的整數(shù)倍，一避免頻譜泄漏。而當(dāng)不知道原信號的周期（或者周期不穩(wěn)定）時(shí)，就要通過延長采樣時(shí)間來盡量減少泄漏，從而突出原信號的譜線。（5） 當(dāng)信號混有白噪聲時(shí)，就更應(yīng)注意減少頻譜的泄漏和混疊，否則信號分析更加困難，甚至可能會(huì)使原信號被誤差“淹沒”。（6） 若原信號有多個(gè)頻率成分，應(yīng)該盡量提高

46、采樣的頻率分辨率，以區(qū)分出更細(xì)微的頻率差異。4.6在實(shí)驗(yàn)中，在計(jì)算2048點(diǎn)采樣時(shí)，初步體會(huì)到了FFT算法的優(yōu)越性。在我的計(jì)算機(jī)上，F(xiàn)FT算法的確比原始的DFT更快。不過由于采樣點(diǎn)數(shù)較少，這一差別僅限于幾秒鐘。在采樣點(diǎn)更多時(shí)，F(xiàn)FT在速度上的優(yōu)越性應(yīng)該能進(jìn)一步突出。4.7實(shí)驗(yàn)中遇到的問題及其解決：實(shí)驗(yàn)中有些M文件代碼總是出錯(cuò)。解決方法：重新檢查，在稿紙上演算，體會(huì)運(yùn)算過程。例如，Matlab數(shù)組序號起始位置（為1，而非C語言的0）的問題，就是這樣發(fā)現(xiàn)的。編程時(shí)對代碼不熟悉，使得思路比較混亂。解決方法：畫流程圖，理清思路。對比較復(fù)雜的程序尤其如此。感到大一時(shí)C語言教學(xué)并未強(qiáng)調(diào)流程圖，這是一個(gè)教學(xué)中值得改進(jìn)的地方。C語言中，比起掌握運(yùn)算符優(yōu)先級，對流程圖思想的培養(yǎng)顯得更為重要。附錄：附錄A值得指出的是，在Matlab中，數(shù)組的序號是從1開始的（