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

上傳人：奔*** IP屬地：河北更新時間：2024-03-02 格式：doc 頁數(shù)：27 大?。?54.00KB 人氣指數(shù)：12 舉報 版權(quán)申訴

已閱讀1頁，還剩26頁未讀，繼續(xù)免費(fèi)閱讀

版權(quán)說明：本文檔由用戶提供并上傳，收益歸屬內(nèi)容提供方，若內(nèi)容存在侵權(quán)，請進(jìn)行舉報或認(rèn)領(lǐng)

文檔簡介

1、　　實驗報告　　課程名稱：信號分析與處理指導(dǎo)老師：成績：__________________　　實驗名稱：離散傅里葉變換和快速傅里葉變換實驗類型：基礎(chǔ)實驗同組學(xué)生姓名： 　　第二次

2、實驗離散傅里葉變換和快速傅里葉變換　　一、實驗?zāi)康?lt;/b>　　1.1掌握離散傅里葉變換（DFT）的原理和實現(xiàn)；　　1.2掌握快速傅里葉變換（FFT）的原理和實現(xiàn)，掌握用FFT對連續(xù)信號和離散信號進(jìn)行譜分析的方法。　　1.3 會用Matlab

3、軟件進(jìn)行以上練習(xí)。　　二、實驗原理　　2.1關(guān)于DFT的相關(guān)知識　　序列x(n)的離散事件傅里葉變換（DTFT）表示為　　，　　如果x(n)為因果有限長序

4、列，n=0,1,...,N-1，則x(n)的DTFT表示為　　，　　x(n)的離散傅里葉變換（DFT）表達(dá)式為　　，　　序列的N點DFT是序列DTFT在頻率區(qū)間[0,2π]上的N點燈間隔采樣，采樣

5、間隔為2π/N。通過DFT，可以完成由一組有限個信號采樣值x(n)直接計算得到一組有限個頻譜采樣值X(k)。X(k)的幅度譜為，其中下標(biāo)R和I分別表示取實部、虛部的運(yùn)算。X(k)的相位譜為。　　離散傅里葉反變換（IDFT）定義為　　。　　2.2關(guān)于FFT的相關(guān)知識</p

6、>　　快速傅里葉變換（FFT）是DFT的快速算法，并不是一個新的映射。FFT利用了函數(shù)的周期性和對稱性以及一些特殊值來減少DFT的運(yùn)算量，可使DFT的運(yùn)算量下降幾個數(shù)量級，從而使數(shù)字信號處理的速度大大提高。　　若信號是連續(xù)信號，用FFT進(jìn)行譜分析時，首先必須對信號進(jìn)行采樣，使之變成離散信號，然后就可以用FFT來對連續(xù)信號進(jìn)行譜分析。為了滿足采樣定理，一般在采樣之前要設(shè)

7、置一個抗混疊低通濾波器，且抗混疊濾波器的截止頻率不得高于與采樣頻率的一半。　　比較DFT和IDFT的定義，兩者的區(qū)別僅在于指數(shù)因子的指數(shù)部分的符號差異和幅度尺度變換，因此可用FFT算法來計算IDFT。　　實驗內(nèi)容與相關(guān)分析（共6道）　　說明：為了便于老師查看，現(xiàn)將各題的內(nèi)容寫在這里——&l

8、t;p>　　題目按照3.1、3.2、...、3.6排列。每道題包含如下內(nèi)容：題干、解答（思路、M文件源代碼、命令窗口中的運(yùn)行及其結(jié)果）、分析。其中“命令窗口中的運(yùn)行及其結(jié)果”按照小題順序排列，各小題包含命令與結(jié)果（圖形或者序列）。　　3.1 求有限長離散時間信號x(n)的離散時間傅里葉變換（DTFT）X(ejΩ)并繪圖。　　已知；

9、（2）已知。　　【解答】　　思路：這是DTFT的變換，按照定義編寫DTFT的M文件即可?？紤]到自變量Ω是連續(xù)的，為了方便計算機(jī)計算，計算時只取三個周期[-2π,4π]中均勻的1000個點用于繪圖。　　理論計算的各序列DTFT表達(dá)式，請見本題的分析。<

10、/p>　　M文件源代碼（我的Matlab源文件不支持中文注釋，抱歉）：　　function DTFT(n1,n2,x)　　%This is a DTFT function for my experiment of Signal Processing & Analysis.　　w

11、=0:2*pi/1000:2*pi;%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);%

12、It is the definition of DTFT.　　end　　Amp=abs(X);%Acquire the amplification.　　Phs=angle(X);%Acquire the phase angle (radian).<

13、;p>　　subplot(1,2,1);　　plot(w,Amp,'r'); xlabel('\Omega');ylabel('Amplification');hold on;　　%Plot amplification on the left.　　subplo

14、t(1,2,2);　　plot(w,Phs,'b');xlabel('\Omega');ylabel('Phase Angle (radian)');hold off;　　%Plot phase angle on the right.　　end&l

15、t;/b>　　命令窗口中的運(yùn)行及其結(jié)果（理論計算的各序列DTFT表達(dá)式，請見本題的分析）：　　第（1）小題　　>> n=(-2:2);　　>> x=1.^n;

16、;　　>>DTFT(-2,2,x);　　第（2）小題　　>> n=(0:10);　　>> x=2.^n;　　>> DTFT(0,10,x);&

17、lt;b>　　【分析】　　對于第（1）小題，由于序列x(n)只在有限區(qū)間（-2，-1，-，1，2）上為1，所以是離散非周期的信號。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號。事實上，我們可計算出它的表達(dá)式：　　，可見幅度頻譜擁有主極大和次極大，兩個主極大間有|5-1|=4個極小，即有3個次級大。而對于它的相位頻譜，則是周期性地在-π、

18、0、π之間震蕩。　　對于第（2）小題，由于是離散非周期的信號。它的幅度頻譜相應(yīng)地應(yīng)該是周期連續(xù)信號。而它的表達(dá)式：，因此主極大之間只有|0-1|=1個極小，不存在次級大。而對于它的相位頻譜，則是在一個長為2π的周期內(nèi)有|11-1|=10次振蕩。　　而由DTFT的定義可知，頻譜都是以2π為周期向兩邊無限延伸的。由于DTFT是連續(xù)譜，對于計算機(jī)處理來說特別困

19、難，因此我們才需要離散信號的頻譜也離散，由此構(gòu)造出DFT（以及為加速計算DFT的FFT）。　　3.2已知有限長序列x(n)={8,7,9,5,1,7,9,5},試分別采用DFT和FFT求其離散傅里葉變換X(k)的幅度、相位圖。　　【解答】　　思路：按照定義編寫M文件即可。&l

20、t;/p>　　M文件源代碼：　　i) DFT函數(shù)：　　function DFT(N,x)　　%This is a DFT function for my experiment of Signal Process

21、ing & 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)*e

22、xp((-1)*j*2*k*pi/N*i);%It is the definition of DFT.　　end　　Amp=abs(X);%Acquire the amplification.　　Phs=angle(X);%Acquire the phase angle (r

23、adian).　　subplot(1,2,1);　　stem(k,Amp,'.',’MarkerSize’,18); xlabel('k');ylabel('Amplification');hold on;　　%Plot amplification on the l

24、eft.　　subplot(1,2,2);　　stem(k,Phs,'*');xlabel('k');ylabel('Phase Angle (radian)');hold off;　　%Plot phase angle on the right.

25、;　　end　　ii) 基2-FFT函數(shù)　　function myFFT(N,x)　　%This is a base-2 FFT function.　　lov=(0:N-1);<p&

26、gt;　　j1=0;　　for i=1:N %indexed addressing　　if i<j1+1　　temp=x(j1+1);　　x(j1+1)=x(i);&

27、lt;p>　　x(i)=temp;　　end　　k=N/2;　　while k<=j1　　j1=j1-k;

28、;　　k=k/2;　　end　　j1=j1+k;　　end　　digit=0;&l

29、t;/p>　　k=N;　　while k>1　　digit=digit+1;　　k=k/2;　　end&l

30、t;/b>　　n=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

31、.　　idx=1;　　for i=1:n　　idx1=idx;　　idx2=1;

32、;　　for j1=1:N/(2*n)　　r=(idx2-1)*2^(digit-mu);　　wn=exp(j*(-2)*pi*r/N);%It is the "circulating coefficients".　　temp=x(idx);　　x(i

33、dx)=temp+x(idx+dif)*wn;　　x(idx+dif)=temp-x(idx+dif)*wn;　　idx=idx+1;　　idx2=idx2+1;　　end　　idx=idx1

34、+2*dif;　　end　　n=n/2;　　end　　Amp=abs(x);%Acquire the amplification.&l

35、t;p>　　Phs=angle(x);%Acquire the phase angle (radian).　　subplot(1,2,1);　　stem(lov,Amp,'.',’MarkerSize’,18);　　xlabel('FFT k');ylabel('FF

36、T Amplification');hold on;　　%Plot the amplification.　　subplot(1,2,2);　　stem(lov,Phs,'*');xlabel('FFT k');ylabel('FFT Phase Angle (rad

37、ian)');hold off;　　end　　命令窗口中的運(yùn)行及其結(jié)果：　　DFT：　　>> x=[8,7,9,5,1,7,9,5];

38、　　>> DFT(8,x);　　FFT：　　>> x=[8,7,9,5,1,7,9,5];　　>> myFFT(8,x);　　【分析】&

39、lt;p>　　DFT是離散信號、離散頻譜之間的映射。在這里我們可以看到序列的頻譜也被離散化。事實上，我們可以循著DFT構(gòu)造的方法驗證這個頻譜：　　首先，將序列做N=8周期延拓，成為離散周期信號。然后利用DFS計算得到延拓后的頻譜：　　，從而取DFS的主值區(qū)間得到DFT，與圖一致。因此計算正確。　　而對于FF

40、T，我們可以看到它給出和DFT一樣的結(jié)果，說明了FFT算法就是DFT的一個等價形式。不過，由于序列不夠長，F(xiàn)FT在計算速度上的優(yōu)越性尚未凸顯。　　3.3已知連續(xù)時間信號x(t)=3cos8πt, X(ω)=，該信號從t=0開始以采樣周期Ts=0.1 s進(jìn)行采樣得到序列x(n)，試選擇合適的采樣點數(shù)，分別采用DFT和 FFT求其離散傅里葉變換X(k)的幅度、相位圖，并將結(jié)果與X(k)的幅度、相位圖

41、，并將結(jié)果與X(ω)相比較。　　【解答】　　思路：此題與下一題都是一樣的操作，可以在編程時統(tǒng)一用變量g（0或1）來控制是否有白噪聲。這里取g=0（無白噪聲）。　　另外，分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏的關(guān)系。<

42、p>　　M文件源代碼：　　i)采樣函數(shù)：　　function xs=sampJune3(N,Ts,g)　　%This is a function applied to Problem 3 & 4.&l

43、t;p>　　%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

44、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 consistent with our habit. Thus I made a adaptation i

45、n codes.　　%It is a sampling process with(g=1)/without(g=0) noise.　　end　　xs=sample;　　plot(n-1,sample,'.',’Mark

46、erSize’,18);xlabel('n');ylabel('value');hold off;　　% Must use (n-1), because in Matlab, index starts at 1.　　end　　ii)DFT和基2

47、-FFT函數(shù)的代碼，請見第3.2節(jié)。不需再新編一個。　　命令窗口中的運(yùn)行及其結(jié)果：　　12點采樣：　　>> xs=sampJune3(12,0.1,0);%末尾的0表示無噪聲。　　>> DFT(12,xs);&

48、lt;/p>　　>> myFFT(12,xs);　　20點采樣：　　>> xs=sampJune3(20,0.1,0);%末尾的0表示無噪聲。　　>> DFT(20,xs);<p&g

49、t;　　>> myFFT(20,xs);　　28點采樣：　　>> xs=sampJune3(28,0.1,0);%末尾的0表示無噪聲。　　>> DFT(28,xs);　　>> myFFT

50、(28,xs);　　【分析】　　分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏之間的關(guān)系?，F(xiàn)在與原信號頻譜比較后可以得出如下結(jié)論：　　原信號的頻譜是，在±8π上各有一強(qiáng)度為3π的譜線，在其余頻率上為0。

51、　　可見原信號被0.1 s采樣周期的采樣信號離散化之后，譜線以20π為周期重復(fù)，并且只在（20k±8）π (k為整數(shù)）處非0。那么，在20點DFT（采樣時間原信號周期的整數(shù)倍）中，只有第8根、第12根譜線非0。而在12點、28點DFT中，由于采樣時間不是原信號周期的整數(shù)倍，譜線將向兩邊泄漏。　　不過，對比12點采樣和28點采樣，我們還可以發(fā)現(xiàn)，28點采樣頻譜的主譜線高度是次譜線高度的4

52、倍，兒12點采樣頻譜的主譜線高度是次譜線高度的3倍?？梢?，在無法保證采樣時間是信號周期整數(shù)倍的情況下，增加采樣時間有助于減輕頻譜泄漏的程度。　　3.4對第3步中所述連續(xù)時間信號疊加高斯白噪聲信號，重復(fù)第3步過程。　　【解答】　　思路：此題與上一題都是一樣的操作，可以在編程時統(tǒng)

53、一用變量g（0或1）來控制是否有白噪聲。這里取g=1（有白噪聲）。　　另外，仍然分別取12點、20點、28點采樣，以考察采樣長度的選擇與頻譜是否泄漏的關(guān)系。　　M文件源代碼：　　不需要再新編程序?？梢灾苯右蒙厦娴暮瘮?shù)：　　sampJ

54、une3(N,Ts,g)，取g=1，以體現(xiàn)存在白噪聲　　DFT(N,x)　　myFFT(N,x)　　命令窗口中的運(yùn)行及其結(jié)果：　　12點采樣：　　>

55、> xs=sampJune3(12,0.1,1);%末尾的1表示有噪聲。　　>> DFT(12,xs);　　>> myFFT(12,xs);　　20點采樣：　　>> xs=sampJune3(

56、20,0.1,1);%末尾的1表示有噪聲。　　>> DFT(20,xs);　　>> myFFT(20,xs);　　28點采樣：　　>> xs=sampJune3(28,0.1,0);%末尾的1表示有

57、噪聲。　　>> DFT(28,xs);　　>> myFFT(28,xs);　　【分析】　　依然分別取12點、20點、28點采樣。仍然與原信號的頻譜（圖3.3.10）比較，可以得到結(jié)論：&

58、lt;p>　　由于疊加了噪聲，所以頻譜都受到了一定的干擾。由于白噪聲在各個頻率的功率相等，因此頻譜上各處的干擾也是均勻隨機(jī)的。　　不過，通過對比我們可以發(fā)現(xiàn)，20點采樣（無噪聲時不發(fā)生泄漏的采樣方法）在存在噪聲時，仍然可以明顯區(qū)分出原信號的譜線。第二好的是28點采樣，因為采樣時間較長，即使存在頻譜泄漏也能較好地區(qū)分原信號的譜線。而最差的是12點采樣，由于噪聲的存在和嚴(yán)重的頻譜泄漏，它的次譜

59、線與主譜線的高度相差不大，使原信號不明顯。　　3.5已知序列，X(k)是x(n)的6點DFT，設(shè)。　　若有限長序列y(n)的6點DFT是，求y(n)。　　若有限長序列w(n)的6點DFT W(k)是的實部，求w(n)。　　若有限長序列q(n)的3點DFT是，k=0,1,2

60、，求q(n)。　　【解答】　　思路：這是對DFT進(jìn)行變換后求IDFT。考慮到IDFT和DFT定義的對稱性，可以在DFT的基礎(chǔ)上略加調(diào)整既可用于計算。　　首先，∵，　　∴它的6點采樣是序列是。值得指

61、出的是，在Matlab中，數(shù)組的序號是從1開始的（而在信號分析中習(xí)慣從0開始），不過我在上面編程時已考慮到這一情況，具體可見實驗報告最后的“附錄”。　　首先生成x(n)的6點DFT，再按照各小題分別轉(zhuǎn)換，最后求相應(yīng)的IDFT。　　M文件源代碼：　　i) 輸出x(n)的6點DF

62、T的函數(shù)：　　function X = ExportDFT(N,x)　　%This is a DFT function that exports the solution to vector X.　　k=(0:N-1); %Define var

63、iable 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

64、 the definition of DFT.　　end　　end　　ii)算第（1）小題的Y(k)的函數(shù)：　　function Y = Convertor1(X)<

65、p>　　%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-

66、1),because in Matlab index starts at 1.　　end　　end　　iii)算第（2）小題的W(k)的函數(shù)：　　function W = Convertor2(X

67、)　　%This is a mathematical convertor for the subproblem (2).　　W=real(X);%Acquire the real part of X.　　end　　iv)算第（3）小題

68、的Q(k)的函數(shù)：　　function Q = Convertor3(X)　　%This is a mathematical convertor for the subproblem (3).　　Q=zeros(3);% Detailed explanation goes here&l

69、t;p>　　for tmp=1:3　　Q(tmp)=X(2*tmp);　　end　　end　　v)輸出IDFT的函數(shù)：　　function x =

70、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.<

71、/p>　　for k=0:N-1　　x=x+X(k+1)*exp(j*2*k*pi/N*n);　　end　　x=x/N;　　a=real(x);%We MUST use

72、 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);xla

73、bel('n');ylabel('Solution');　　end　　命令窗口中的運(yùn)行及其結(jié)果：　　>> x=[4,3,2,1,0,0];　　>> X=ExportDFT(6,x

74、);　　第（1）小題　　>> Y=Convertor1(X);　　>> y=ExportIDFT(6,Y)　　y =　　Colum

75、ns 1 through 4　　0.0000 + 0.0000i 0.0000 + 0.0000i 4.0000 + 0.0000i 3.0000 + 0.0000i　　%虛部都是0，說明是實數(shù)　　Columns 5 through 6　　2.0000 + 0.

76、0000i 1.0000 - 0.0000i %虛部都是0，說明是實數(shù)　　%事實上，在Matlab中，由于數(shù)值計算的截斷誤差，對原復(fù)數(shù)做乘法后，答案的虛部可能有一極小的量。　　第（2）小題　　>> W=Convertor2(X);&l

77、t;p>　　>> 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&l

78、t;/p>　　%虛部都是0，說明是實數(shù)　　Columns 5 through 6　　1.0000 + 0.0000i 1.5000 + 0.0000i %虛部都是0，說明是實數(shù)；　　%事實上，在Matlab中，由于數(shù)值計算的截斷誤差，對原復(fù)數(shù)做乘法后，答案的虛部可能有一極小的量。&

79、lt;/p>　　第（3）小題　　>> Q=Convertor3(X);　　>> q=ExportIDFT(6,Q)　　q =　　Columns

80、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

81、　　這里的答案都是幅值、相位均非0的復(fù)數(shù)，而教材（實驗指導(dǎo)第109頁）并未要求作圖，這里略去。　　答案：q(n)={1.5，-0.1667-0.2887i，0.7500-1.2990i，0.8333，-0.5000-0.8660i，1.0833-1.8764i}　　【分析】&

82、lt;p>　　對原序列進(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）小題，由于對原頻譜取了實部，那么根據(jù)DFT的奇偶虛實性知，

83、得到的w(n)也是實數(shù)的。　　第（3）小題，對原信號進(jìn)行了尺度變換（“抽取”），導(dǎo)致丟失了一些譜線，使得無法通過IDFT得到原來的序列x(n)。說明頻譜記錄了原有信號的信息，若頻譜發(fā)生變化，則對應(yīng)的時域信號也隨之改變。　　3.6已知信號，其中f1=4 Hz、f2=4.02 Hz、f3=5 Hz，采用采樣頻率為20 Hz進(jìn)行采樣，求

84、　　當(dāng)采樣長度N分別為512和2048情況下x(t)的幅度頻譜；　　當(dāng)采樣長度N為32，且增補(bǔ)N個零點、4N個零點、8N個零點、16N個零點情況下x(t)的幅度頻譜。　　【解答】　　思路：采樣是有限且離散的，用DFT（FFT算法）計算頻譜，以便得到離散的頻譜

85、，并且具有較高速度。　　20 Hz對應(yīng)的采樣周期Ts=0.05s。　　M文件源代碼：　　i)采樣函數(shù)（其中Plus表示采樣后補(bǔ)零的個數(shù)）　　function xs=sampJune6(N,Plus)<p&

86、gt;　　%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

87、;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+

88、1):(N+Plus)　　sample(p)=0;%Add zero points.　　end　　xs=sample;%Return　　end　　ii)由

89、于只要求顯示幅度頻譜，所以刪去myFFT函數(shù)中繪制相位頻譜的命令，使它的最后部分修改如下：　　原來的：　　function myFFT(N,x)　　%This is a base-2 FFT function wrote by myself.<p

90、>　　...　　Amp=abs(x);%Acquire the amplification.　　Phs=angle(x);%Acquire the phase angle (radian).　　subplot(1,2,1);

91、stem(lov,Amp,'.');xlabel('FFT k');ylabel('FFT Amplification');hold on;　　%Plot the amplification.　　subplot(1,2,2);　　stem(lov,Phs,'

92、*');xlabel('FFT k');ylabel('FFT Phase Angle (radian)');hold off;　　end　　修改后的：　　function myFFT(N,x)&l

93、t;/p>　　%This is a base-2 FFT function wrote by myself.　　...　　Amp=abs(x);%Acquire the amplification.　　stem(lov,Amp,'.'

94、;);xlabel('FFT k');ylabel('FFT Amplification');　　%Plot the amplification.　　end　　命令窗口中的運(yùn)行及其結(jié)果：

95、;　　第（1）小題　　>> x512=sampJune6(512,0);　　>> x2048=sampJune6(2048,0);　　>> myFFT(512,x512);　　>> myFFT(2048,

96、x2048);　　第（2）小題　　>> x32p1N=sampJune6(32,32*1);%32點采樣，補(bǔ)零N=32個，共64個數(shù)據(jù)點　　>> x32p4N=sampJune6(32,32*4);%32點采樣，補(bǔ)零4N=128個，共160個數(shù)

97、據(jù)點　　>> x32p8N=sampJune6(32,32*8);%32點采樣，補(bǔ)零8N=256個，共288個數(shù)據(jù)點　　>> x32p16N=sampJune6(32,32*16);%32點采樣，補(bǔ)零16N=128個，共544個數(shù)據(jù)點　　>> myFFT(32+32*1,

98、x32p1N);　　>> myFFT(32+32*4,x32p4N);　　>> myFFT(32+32*8,x32p8N);　　>> myFFT(32+32*16,x32p16N);　　【分析】請注意，題目只要求繪制幅度頻譜。</

99、p>　　第（1）小題：　　首先，由于采樣時間都不是原有信號周期的整數(shù)倍，兩個采樣方式對應(yīng)的頻譜均發(fā)生了泄漏。不過由于2048點采樣對應(yīng)的采樣時間較長，它頻譜泄漏的程度比512點采樣輕。　　其次，由于20 Hz的2048點采樣的頻率分辨率為20/2048=0.0098 Hz < 0.2

100、 Hz，因此放大頻譜圖之后我們可以看到4 Hz、4.02 Hz和5 Hz對應(yīng)的譜線，而512點采樣的頻率分辨率為20/512=0.039 Hz > 0.2 Hz，因此4 Hz和4.02 Hz對應(yīng)的譜線無法區(qū)分。　　第（2）小題：　　首先，由于采樣時間都不是原有信號周期的整數(shù)倍，頻譜均發(fā)生了泄漏。而且由于采樣時間

101、較短，頻譜泄漏比第（1）小題的兩個頻譜更加嚴(yán)重。　　其次，由于都是32點采樣，因此實際的頻率分辨率較低，無法區(qū)分4 Hz和4.02 Hz對應(yīng)的譜線。　　最后，雖然都是32點采樣，但由于補(bǔ)0個數(shù)的不同，各頻譜譜線間距各不相同。例如，補(bǔ)零最多的序列一共有544個數(shù)據(jù)點，因此譜線間距小。由此還可以得出結(jié)論：數(shù)據(jù)點個數(shù)越多，則頻譜越傾向于連續(xù)。

102、;　　可見，當(dāng)采樣時間不是原信號周期整數(shù)倍而且采樣時間較短時，頻譜泄漏相當(dāng)嚴(yán)重的，所有的頻率上都有了幅值即能量，可見當(dāng)取樣信號的樣點數(shù)取得不夠時，原信號所攜帶的信息就不能被完全取得。而若將取樣信號補(bǔ)零，由圖可見信號的能量相應(yīng)的泄漏到了幾乎所有頻率上了，這樣所得的信號仍然嚴(yán)重失真，因此不能靠將信號補(bǔ)零這樣的方法來取得更精確的采樣信號。要想獲得不泄漏的頻譜，在采樣頻率不變的前提下，必須使采樣時間等于原信號周期的整數(shù)倍，

103、或者盡量延長采樣時間以減少泄漏。　　四、實驗體會　　4.1關(guān)于各個實驗的分析，請見第3部分每道題的末尾。　　4.2在Matlab中，數(shù)組的序號是從1開始的，這與信號處理時通常的序號起點（0）不一致。我在編程充分注意到了這個問題。　　4

104、.3由于Matlab進(jìn)行數(shù)值計算的過程中存在截斷誤差，所以最后算得的值并不是準(zhǔn)確值。例如，對一個復(fù)數(shù)z，即使f(z)的虛部為零，但由于截斷誤差的存在（特別是z的虛部為無窮小數(shù)時），最終f(z)值的虛部可能是一個極小的非零值，從而在顯示時出現(xiàn)“零虛部”（例如，2+0.0000i )。　　4.4通過利用軟件對離散信號的各種變換DTFT、DFT以及其快速算法FFT進(jìn)行計算，使得在實驗中比較難以實現(xiàn)的信

105、號分析過程（離散信號的采集和顯示都是比較困難的）在計算機(jī)計算中實現(xiàn)，證明了理論的正確性，說明仿真計算是一種十分有效的輔助手段。　　4.5通過這次實驗和上次實驗《信號的采集與恢復(fù)》我知道了，要想盡量不失真地取得一個信號的頻譜（低混疊、低泄漏），應(yīng)該盡量滿足以下條件：　　使用的開關(guān)函數(shù)要盡量接近理想沖激串；　　采樣頻

106、率要高于原始信號的奈奎斯特頻率。對于頻譜不受限的信號，為了避免頻譜混疊，應(yīng)該使用低通濾波器進(jìn)行濾波；　　對于頻帶不受限的信號，抗混疊濾波器要盡量接近理想濾波器。　　采樣的持續(xù)時間最好能夠是原信號周期的整數(shù)倍，一避免頻譜泄漏。而當(dāng)不知道原信號的周期（或者周期不穩(wěn)定）時，就要通過延長采樣時間來盡量減少泄漏，從而突出原信號的譜線。<

107、;p>　　當(dāng)信號混有白噪聲時，就更應(yīng)注意減少頻譜的泄漏和混疊，否則信號分析更加困難，甚至可能會使原信號被誤差“淹沒”。　　若原信號有多個頻率成分，應(yīng)該盡量提高采樣的頻率分辨率，以區(qū)分出更細(xì)微的頻率差異。　　4.6在實驗中，在計算2048點采樣時，初步體會到了FFT算法的優(yōu)越性。在我的計算機(jī)上，F(xiàn)FT算法的確比原始的DFT更快。不過由于采樣點數(shù)較少，

108、這一差別僅限于幾秒鐘。在采樣點更多時，F(xiàn)FT在速度上的優(yōu)越性應(yīng)該能進(jìn)一步突出。　　4.7實驗中遇到的問題及其解決：　　實驗中有些M文件代碼總是出錯。解決方法：重新檢查，在稿紙上演算，體會運(yùn)算過程。例如，Matlab數(shù)組序號起始位置（為1，而非C語言的0）的問題，就是這樣發(fā)現(xiàn)的。　　編程時對代碼不熟悉，使得思路比較

109、混亂。解決方法：畫流程圖，理清思路。對比較復(fù)雜的程序尤其如此。感到大一時C語言教學(xué)并未強(qiáng)調(diào)流程圖，這是一個教學(xué)中值得改進(jìn)的地方。C語言中，比起掌握運(yùn)算符優(yōu)先級，對流程圖思想的培養(yǎng)顯得更為重要。　　附錄：　　附錄A值得指出的是，在Matlab中，數(shù)組的序號是從1開始的（而在信號分析中習(xí)慣從0開始），我在上面編程時已考

110、慮到這一情況。例如，下面可以驗證DFT函數(shù)的正確性（順便也可以驗證FFT函數(shù)）。　　>> n=1:8;　　>> xn=6*cos(pi*(n-1)/4);　　>> DFT(8,xn);　　附錄B出

111、于興趣，我們還可以看看DFT和myFFT的運(yùn)行時間。　　為了不計繪圖，刪去DFT和myFFT的繪圖命令。通過Matlab命令窗口查看運(yùn)行時間。　?。?）16（24）點采樣：　　>> tic;n=(0:15);xs=sin(0.1*2*pi*n)+randn(size(n));DFT(16,xs);to

112、c;　　Elapsed time is 0. seconds.　　>> tic;n=(0:15);xs=sin(0.1*2*pi*n)+randn(size(n));myFFT(16,xs);toc;　　Elapsed time is 0. seconds.

113、?。?）2048（211）點采樣：　　>> tic;n=(0:2047);xs=sin(0.1*2*pi*n)+randn(size(n));DFT(2048,xs);toc;　　Elapsed time is 0. seconds.　　>> tic;n=(0:2047);xs=sin(

114、0.1*2*pi*n)+randn(size(n));myFFT(2048,xs);toc;　　Elapsed time is 0. seconds.　?。?）8192（213）點采樣：　　>> tic;n=(0:8191);xs=sin(0.1*2*pi*n)+randn(size(n));DFT

115、(8192,xs);toc;　　Elapsed time is 8. seconds.　　>> tic;n=(0:8191);xs=sin(0.1*2*pi*n)+randn(size(n));myFFT(8192,xs);toc;　　Elapsed time is 0. seconds.</

116、p>　?。?）65536（216）點采樣：　　>> tic;n=(0:65535);xs=sin(0.1*2*pi*n)+randn(size(n));DFT(65536,xs);toc;　　Elapsed time is 2192. seconds.　　>>

117、tic;n=(0:65535);xs=sin(0.1*2*pi*n)+randn(size(n));myFFT(65536,xs);toc;　　Elapsed time is 4. seconds.　　測試環(huán)境配置：Matlab 2013a，Windows 7 Professional，Intel Centrino2雙核32位1.40 GHz，內(nèi)存3.0

118、 G。　　可見，在采樣點數(shù)較小時，可能由于反復(fù)賦值、變址等運(yùn)算，F(xiàn)FT算法的運(yùn)行速度較慢。　　而在采樣點數(shù)較大時，F(xiàn)FT算法在速度上的優(yōu)越性凸顯出來。例如在8192點采樣時，F(xiàn)FT算法運(yùn)行時間僅約為原始DFT算法的5.9%；而在65536點采樣時，原始DFT算法運(yùn)行時間達(dá)到了近37分鐘，而FFT算法運(yùn)行時間僅約4.4 s，為前者的0.2%，優(yōu)勢相當(dāng)明顯。

眾賞文庫> 全部分類> 畢業(yè)設(shè)計

溫馨提示

1. 本站所有資源如無特殊說明，都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
2. 本站的文檔不包含任何第三方提供的附件圖紙等，如果需要附件，請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
3. 本站RAR壓縮包中若帶圖紙，網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽，若沒有圖紙預(yù)覽就沒有圖紙。
4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
5. 眾賞文庫僅提供信息存儲空間，僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理，對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯，并不能對任何下載內(nèi)容負(fù)責(zé)。
6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容，請與我們聯(lián)系，我們立即糾正。
7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

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

文檔簡介

溫馨提示

最新文檔

評論

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

文檔簡介

溫馨提示

最新文檔

評論

免費(fèi)下載