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

1、實驗報告課程名稱： 信號分析與處理 指導老師： 成績：_實驗名稱：離散傅里葉變換和快速傅里葉變換 實驗類型： 基礎實驗 同組學生姓名： 第二次實驗 離散傅里葉變換和快速傅里葉變換裝 訂 線一、實驗目的1.1掌握離散傅里葉變換（DFT）的原理和實現；1.2掌握快速傅里葉變換（FFT）的原理和實現，掌握用FFT對連續(xù)信號和離散信號進行譜分析的方法。1.3 會用Matlab軟件進行以上練習。二、實驗原理2.1關于DFT的相關知識序列x(n)的離散事件傅里葉變換（DTFT）表示為，如果x(n)為因果有限長序列，n=0,1,.,N-1，則x(n)的DTFT表示為，x(n)的離散傅里葉變換（DFT）表達式

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

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

4、間信號x(n)的離散時間傅里葉變換（DTFT）X(ej)并繪圖。（1） 已知；（2）已知?！窘獯稹克悸罚哼@是DTFT的變換，按照定義編寫DTFT的M文件即可?？紤]到自變量是連續(xù)的，為了方便計算機計算，計算時只取三個周期-2,4中均勻的1000個點用于繪圖。理論計算的各序列DTFT表達式，請見本題的分析。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命令窗口中的運行及其結果（理論計算的各序列DTFT表達式，請見本題的分析）：第（1）小題>> n=(-2:2);>&g

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

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

9、文件即可。M文件源代碼：i) DFT函數：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函數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命令窗口中的運行及其結果：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是離散信號、離散頻譜之間的映射。在這里我們可以看到序列的頻譜也被離散化。事實上，我們可以循著DFT構造的方法驗證這個頻譜：首先，將序列做N=8周期延拓，成為離散周期信號。然后利用DFS計算得到延拓后的頻譜：，從而取DFS的主值區(qū)間得到DFT，與圖一致。因此計算正確。而對于FFT，我們可以看到它給出和DFT一樣的結果，說明了FFT算法就是DFT的一個等價形式。不過，由于序列不夠長，FFT在計算速度上的優(yōu)越性尚未凸顯。3.3已知連續(xù)

16、時間信號x(t)=3cos8t, X()=，該信號從t=0開始以采樣周期Ts=0.1 s進行采樣得到序列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)%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函數的代碼，請見第3.2

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

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

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

22、是信號周期整數倍的情況下，增加采樣時間有助于減輕頻譜泄漏的程度。3.4對第3步中所述連續(xù)時間信號疊加高斯白噪聲信號，重復第3步過程?！窘獯稹克悸罚捍祟}與上一題都是一樣的操作，可以在編程時統一用變量g（0或1）來控制是否有白噪聲。這里取g=1（有白噪聲）。 另外，仍然分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏的關系。M文件源代碼：不需要再新編程序?？梢灾苯右蒙厦娴暮瘮担簊ampJune3(N,Ts,g)，取g=1，以體現存在白噪聲DFT(N,x)myFFT(N,x)命令窗口中的運行及其結果：12點采樣：>> xs=sampJune3(12,0.1,1);%

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

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

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

26、基礎上略加調整既可用于計算。首先，它的6點采樣是序列是。值得指出的是，在Matlab中，數組的序號是從1開始的（而在信號分析中習慣從0開始），不過我在上面編程時已考慮到這一情況，具體可見實驗報告最后的“附錄”。 首先生成x(n)的6點DFT，再按照各小題分別轉換，最后求相應的IDFT。M文件源代碼：i) 輸出x(n)的6點DFT的函數：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)的函數：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)的函數：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)的函數：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的函數：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命令窗口中的運行及其結果：>> 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，說明是實數 Columns 5 through 6 2.0000 + 0.0000i 1.0000 - 0.0000i %虛部都是0，說明是實數%事實上，在Matlab中，由于數值計算的截斷誤差，對原復數做乘法后，答案的虛部可能有一極小的量。答案：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，說明是實數 Columns 5 through 6 1.0000 + 0.0000i 1.5000 + 0.0000i %虛部都是0，說明是實數；%事實上，在Matlab中，由于數值計算的截斷誤差，對原復數做乘法后，答案的虛部可能有一極小的量。答案：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的復數，而教材（實驗指導第109頁）并未要求作圖，這里略去。 答案：q(n)=1.5，-0.1667-0.2887

34、i，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的性質可以判斷，對原頻譜乘上旋轉因子之后進行IDFT得到的y(n)，就是對原序列做圓周移位：。第（2）小題，由于對原頻譜取了實部，那么根據DFT的奇偶虛實性知，得到的w(n)也是實數的。第（3）小題，對原信號進行了尺度變換（“抽取”），導致丟失了一些譜線，使得無法通過IDFT得到原來的序列x(n)。說明頻譜記錄了原有信

35、號的信息，若頻譜發(fā)生變化，則對應的時域信號也隨之改變。3.6已知信號，其中f1=4 Hz、f2=4.02 Hz、f3=5 Hz，采用采樣頻率為20 Hz進行采樣，求（1） 當采樣長度N分別為512和2048情況下x(t)的幅度頻譜；（2） 當采樣長度N為32，且增補N個零點、4N個零點、8N個零點、16N個零點情況下x(t)的幅度頻譜?！窘獯稹克悸罚翰蓸邮怯邢耷译x散的，用DFT（FFT算法）計算頻譜，以便得到離散的頻譜，并且具有較高速度。20 Hz對應的采樣周期Ts=0.05s。M文件源代碼：i)采樣函數（其中Plus表示采樣后補零的個數）function xs=sampJune6(N,Plu

36、s)%This is a function applied to Problem 6%N:samling 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;%Returnendi

37、i)由于只要求顯示幅度頻譜，所以刪去myFFT函數中繪制相位頻譜的命令，使它的最后部分修改如下：原來的：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

38、9;);hold on;%Plot the amplification.subplot(1,2,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(&#

39、39;FFT k');ylabel('FFT Amplification');%Plot the amplification.end命令窗口中的運行及其結果：第（1）小題>> x512=sampJune6(512,0);>> x2048=sampJune6(2048,0);>> myFFT(512,x512);>> myFFT(2048,x2048);圖3.6.1 進行512點采樣得到的頻譜圖3.6.2 進行2048點采樣得到的頻譜第（2）小題>> x32p1N=sampJune6(32,32*1);%32點采

40、樣，補零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);>>

41、; myFFT(32+32*16,x32p16N);圖3.6.3 采樣N=32點，補零N=32點，共64點的頻譜圖3.6.4 采樣N=32點，補零4N=128點，共160點的頻譜圖3.6.5 采樣N=32點，補零8N=32點，共288點的頻譜圖3.6.6 采樣N=32點，補零16N=32點，共544點的頻譜【分析】請注意，題目只要求繪制幅度頻譜。第（1）小題：首先，由于采樣時間都不是原有信號周期的整數倍，兩個采樣方式對應的頻譜均發(fā)生了泄漏。不過由于2048點采樣對應的采樣時間較長，它頻譜泄漏的程度比512點采樣輕。其次，由于20 Hz的2048點采樣的頻率分辨率為20/2048=0.0098

42、Hz < 0.2 Hz，因此放大頻譜圖之后我們可以看到4 Hz、4.02 Hz和5 Hz對應的譜線，而512點采樣的頻率分辨率為20/512=0.039 Hz > 0.2 Hz，因此4 Hz和4.02 Hz對應的譜線無法區(qū)分。第（2）小題：首先，由于采樣時間都不是原有信號周期的整數倍，頻譜均發(fā)生了泄漏。而且由于采樣時間較短，頻譜泄漏比第（1）小題的兩個頻譜更加嚴重。其次，由于都是32點采樣，因此實際的頻率分辨率較低，無法區(qū)分4 Hz和4.02 Hz對應的譜線。最后，雖然都是32點采樣，但由于補0個數的不同，各頻譜譜線間距各不相同。例如，補零最多的序列一共有544個數據點，因此譜線間

43、距小。由此還可以得出結論：數據點個數越多，則頻譜越傾向于連續(xù)?？梢?，當采樣時間不是原信號周期整數倍而且采樣時間較短時，頻譜泄漏相當嚴重的，所有的頻率上都有了幅值即能量，可見當取樣信號的樣點數取得不夠時，原信號所攜帶的信息就不能被完全取得。而若將取樣信號補零，由圖可見信號的能量相應的泄漏到了幾乎所有頻率上了，這樣所得的信號仍然嚴重失真，因此不能靠將信號補零這樣的方法來取得更精確的采樣信號。要想獲得不泄漏的頻譜，在采樣頻率不變的前提下，必須使采樣時間等于原信號周期的整數倍，或者盡量延長采樣時間以減少泄漏。四、實驗體會4.1關于各個實驗的分析，請見第3部分每道題的末尾。4.2在Matlab中，數組的

44、序號是從1開始的，這與信號處理時通常的序號起點（0）不一致。我在編程充分注意到了這個問題。4.3由于Matlab進行數值計算的過程中存在截斷誤差，所以最后算得的值并不是準確值。例如，對一個復數z，即使f(z)的虛部為零，但由于截斷誤差的存在（特別是z的虛部為無窮小數時），最終f(z)值的虛部可能是一個極小的非零值，從而在顯示時出現“零虛部”（例如，2+0.0000i )。4.4通過利用軟件對離散信號的各種變換DTFT、DFT以及其快速算法FFT進行計算，使得在實驗中比較難以實現的信號分析過程（離散信號的采集和顯示都是比較困難的）在計算機計算中實現，證明了理論的正確性，說明仿真計算是一種十分有效

45、的輔助手段。4.5通過這次實驗和上次實驗信號的采集與恢復我知道了，要想盡量不失真地取得一個信號的頻譜（低混疊、低泄漏），應該盡量滿足以下條件：（1） 使用的開關函數要盡量接近理想沖激串；（2） 采樣頻率要高于原始信號的奈奎斯特頻率。對于頻譜不受限的信號，為了避免頻譜混疊，應該使用低通濾波器進行濾波；（3） 對于頻帶不受限的信號，抗混疊濾波器要盡量接近理想濾波器。（4） 采樣的持續(xù)時間最好能夠是原信號周期的整數倍，一避免頻譜泄漏。而當不知道原信號的周期（或者周期不穩(wěn)定）時，就要通過延長采樣時間來盡量減少泄漏，從而突出原信號的譜線。（5） 當信號混有白噪聲時，就更應注意減少頻譜的泄漏和混疊，否則信

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