DFT在信號(hào)頻譜分析中的應(yīng)用

DFT在信號(hào)頻譜分析中的應(yīng)用目錄TOC\o"1-1"\h\z\uⅠ.設(shè)計(jì)題目 1Ⅱ.設(shè)計(jì)目的 1Ⅲ.設(shè)計(jì)原理 1Ⅳ.實(shí)現(xiàn)方法 1Ⅴ.設(shè)計(jì)內(nèi)容及結(jié)果 5Ⅵ.改進(jìn)及建議 11Ⅶ.思考題及解答 14Ⅷ.設(shè)計(jì)體會(huì)及心得 15Ⅸ.參考文獻(xiàn) 16Ⅰ.設(shè)計(jì)題目DFT在信號(hào)頻譜分析中的應(yīng)用Ⅱ.設(shè)計(jì)目的MatlabDFTDFT應(yīng)用，用DFT對(duì)序列進(jìn)行頻譜分析，了解DFTFFT的應(yīng)用。Ⅲ.設(shè)計(jì)原理DFT是一種時(shí)工程實(shí)際中，經(jīng)常遇到的連續(xù)信號(hào)Xa(t),其頻譜函數(shù)Xa(jW)也是連續(xù)函數(shù)。數(shù)字計(jì)算機(jī)難于處理，因而我們采用DFT來(lái)對(duì)連續(xù)時(shí)間信號(hào)的傅里葉變換進(jìn)行逼近，進(jìn)而分析連續(xù)時(shí)間信號(hào)的頻譜。Ⅳ.實(shí)現(xiàn)方法葉級(jí)數(shù)本質(zhì)是一樣的。（1）N的序列DFTDFT（2）WN(nk)的周期性2DFT（對(duì)稱性Wnk2N1

Wnk，N2WN1;周期性Wn(rNk)2

WnrNW

Wnk，r為任意整數(shù),WnrN1）N N N N N N離散傅里葉變換的推導(dǎo)：離散傅里葉級(jí)數(shù)定義為xp

(n)

1N1xN k0

(k

jN

（1-1）N將上式兩端乘以ej2πnmN

并對(duì) n 在

0~N-1 N1

(n)ej2πnm

1NN1

(k)ej(km

N1

(k)

N1ejπn(km)n0

Np1 N1

Nn0k0

pN11-N

Nj2π(km)N

pk0

1 N n0 因?yàn)? NNn0N

jn(kmN

N NNN 1-ej(km)NN

1km0 kmN1xpn0

(n)ej

N1k0

(k(km) p

N1xp Nn0N

ej

用k代m

(k)P

Nn0

x(n)ep

jN

（1-2）令WN

ej2π則（1-2）DFS

x(n)Xp

(k)N1xp n0

(n)WN

（1-3） 1

N1 （1-1）IDFS

(k)xp

(n)

Nn0

X (knkp N

（1-4）式1-3（1-4）式構(gòu)成周期序列傅里葉級(jí)數(shù)變換關(guān)系。其中xp

(n)、X

(k)p都是周期為N的周期序列，DFS[·]表示離散傅里葉級(jí)數(shù)正變換，IDFS[·]表示離散傅里葉級(jí)數(shù)反變換。習(xí)慣上，對(duì)于長(zhǎng)為N的周期序列，把0nN-1區(qū)間稱為主值區(qū)，把x(0)~xp

(N1)稱為xp

(n

(0)~p

(N1)p

(k)的p主值序列。由于x(n)xp

(n)RN

(nxp

(nN個(gè)獨(dú)立樣值，對(duì)于任何一個(gè)周期進(jìn)行研究就可以得到它的全部信息。在主值區(qū)研究xp

(nx(n)是等價(jià)的，因此在主值區(qū)計(jì)算DFSDFT是相等的，所以DFT計(jì)算公式形式與DFS基本相同。其關(guān)系為x(n)xp所以離散傅里葉正變換

(n)RN

(n) X(k)

(k)Rp

(k)XkDFTxnN1xnWnkN

0kN-1n02離散傅里葉變換DF）定義：x(n)（0nN-N的頻率有限長(zhǎng)序列（0kN-，其正變換為XkDFTxnN1xnWnkN

W0kN-1 （

N

j）Nn0離散傅里葉變換的實(shí)質(zhì)是：把有限長(zhǎng)序列當(dāng)做周期序列的主值序列進(jìn)行DFSX(k)N，都是Nx(n)X(k),X(k)x(n)。DFS導(dǎo)出來(lái)的，因而隱含著周期性。構(gòu)造離散傅里葉變換的Matlab實(shí)現(xiàn)程序如下：function[Xk]=dft(xn,N)n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnkDFT不同的另外一種變換，而是為了減少DFT計(jì)算次數(shù)的一種快速有效的算法共軛對(duì)稱性：~設(shè)有限長(zhǎng)序列x(n)的長(zhǎng)度為N，以N為周期的周期延拓列為x(n)x((n))~N~ x(n)的共軛對(duì)稱分量(n)~

(n)分別為e o~ 1~ ~* 1 2 x(n) x(n)x(n) x((n))2 e

x*((Nn))N

(1-5)~ 1~ ~* 1 2 x(n) x(n)x(n) x((n))2 o

x*((Nn))N

(1-6)~ ~* ~ ~xe

(n)xe

(n) （1-7）xo

(n)x

*(n) （1-8）o則有限長(zhǎng)序列x(n)的圓周共軛對(duì)稱分量xep

(nxop

(n)分別定義為：3x (n)~ep

(n)RN

(n)1[x((n))2

x*((Nn)) ]RN N

(n) （1-9）~ 1x (n)xop

(n)RN

(n) [x((n))2

x*((Nn)) ]RN N

(n) （1-10）由于滿足~(n)~(n)~ (n) 故e o~ ~ ~x(n)x(n)RN

(n)[xe

(n)x(n)]RN

(n)xep

(n)xop

(n) Nx(n)xep

(n)和圓周共軛反對(duì)稱分量xop

(n)之和，xep

(n)xop

N。利用有限長(zhǎng)序列與周期序列的共軛對(duì)稱分量和反對(duì)稱分量的關(guān)系式1-9）和式1-1，以及式（1-11）可以推導(dǎo)出DFT的一系列的對(duì)稱性質(zhì)DFT[x*(n)]X*(k)X*(nK) x*(nx(n的共軛復(fù)序列。證明：DFT [x*(n)]

N

x*

N

x(n)WnkX*(k) 又因?yàn)閚0

Nn0

NN1 *W

ej

)nNej2πn

1 DFT[x*(n)]

x(n)W(Nk)n X*(Nk)N Nn0NDFT等于DFT的圓周共軛對(duì)稱部分，即N1DFT{Re[x(n)]}Xep證明：

(k) [X(k)X*(Nk2DFT {Re[x(n)]} DFT {1[x(n)x*(n)]} = 1 {DFT [x(n)] +DFT [x*(n)] }=1[X(k)X*(Nk)]X2

2 2(k)利用DFT的對(duì)稱性可求得cosn的DFT:0x(n)cos0

njsin0

ne則DFT[x(n)] X(k)

N1 jenen

nk

1ejoNWNkN

1

joN因?yàn)閏osnRe[x(n)]0所以

n0

N 1ejoWN

1ejoWkNDFT[cos

n]DFT{Re[x(n)]}X0

(k)X(k)X*(Nk)=24[1ejoN 1ejoWk

1cos N11ejoWk1ejoN]N2

kcos

W

cos( N1)N N12WkN

0cos0

N 0W2kNⅤ.設(shè)計(jì)內(nèi)容及結(jié)果MATLABx(n)NDFTmdft.mMATLAB中的內(nèi)部函數(shù)文件fft.m作比較。解：x(n)NDFT的mfunction[Xk]=dft(xn,N)n=[0:1:N-1];k=n;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnkMatlab中的內(nèi)部函數(shù)文件fft.mfunction[varargout]fft(varargin)ifnargout==0builtin('fft',else[varargout{1:nargout}]=builtin('fft',end運(yùn)算量估計(jì)：N=2MFFTMN/2個(gè)蝶。每一級(jí)需N/2次復(fù)乘、N次復(fù)加，因此總共需要進(jìn)行：2

MN2

log2

N 復(fù)加：NMNlog N2直接計(jì)算N點(diǎn)的DFT，需要N2次復(fù)乘、N(N-1)次復(fù)加。N值越大，時(shí)間抽選奇FFT算法越優(yōu)越。例如當(dāng)N=2048點(diǎn)時(shí)，時(shí)間抽選奇偶分解FFTDFT300多倍可以用一下Matlab程序比較DFT和FFT的運(yùn)算時(shí)間N=2048;5M=11;x=[1:M,zeros(1,N-M)];t=cputime;y1=fft(x,N);Time_fft=cputime-tt1=cputime;y2=dft(x,N);Time_dft=cputime-t1t2=cputime;運(yùn)行結(jié)果：Time_fft=0.0469Time_dft=15.2031由此可見(jiàn)FFT算法比直接計(jì)算DFT速度快得多對(duì)離散確定信號(hào)x(n)作如下譜分析：x(n)x(n)N(0nN-1N自己選)寫(xiě)程序計(jì)算出x(n)的N點(diǎn)DFTX(k) ,畫(huà)出時(shí)域序列圖xn～n和相應(yīng)的幅頻圖X(k)~k。解：x(n)10X(ejw)、X(k)MATLAB程序如下：N=10;n=0:1:N-1;xn=cos(0.48*pi*n)+cos(0.52*pi*n);Xk=dft(xn,N);subplot(3,1,1)stem(n,xn,'.k');title('時(shí)域序列圖xn');6xlabel('n');axis([0,10,-2.5,2.5]);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-j*n'*w);subplot(3,1,2);plot(w/pi,abs(Xw));title('幅頻特性曲線X(ejw)');xlabel('w');axis([0,1,0,10]);subplot(3,1,3)k1=0:1:9;w1=2*pi/10*k1;stem(w1/pi,abs(Xk),'.k');title('Xk');xlabel('頻率（axis([0,1,0,10]);x(n)的前10點(diǎn)數(shù)據(jù)對(duì)應(yīng)的x(n)、X(ejw)、X(k)如圖1-1所示。圖1-1 x(n)的前10點(diǎn)數(shù)據(jù)對(duì)應(yīng)的x(n)、X(ejw)、X(k)7由圖可見(jiàn)，由于截?cái)嗪瘮?shù)的頻譜混疊作用，X(k)不能正確分辨w1=0.48π、w2=0.52π這兩個(gè)頻率分量。x(n)100N=100X(ejw)、X(k)MATLAB主要程序如下：N=10;n=0:N-1;xn=cos(0.48*pi*n)+cos(0.52*pi*n);N1=100;n1=0:N1-1;x1=[xn(1:10)zeros(1,90)];subplot(3,1,1)stem(n1,x1,'.k');title('時(shí)域序列圖x1');xlabel('n');axis([0,100,-2.5,2.5]);w=2*pi*(0:2047)/2048;X1=x1*exp(-j*n1'*w);subplot(3,1,2);plot(w/pi,abs(X1));title('幅頻特性曲線X(ejw)');xlabel('w');axis([0,1,0,10]);subplot(3,1,3)Xk=dft(x1,N1);k1=0:1:49;w1=2*pi/100*k1;stem(w1/pi,abs(Xk(1:1:50)),'.k');title('頻域序列圖Xk');xlabel('頻率（單位：pi）');axis([0,1,0,10]);x(n)補(bǔ)零至100點(diǎn)對(duì)應(yīng)的x(n)、X(ejw)、X(k)如圖1-2所示。8圖1-2 x(n)補(bǔ)零至100點(diǎn)對(duì)應(yīng)的x(n)、X(ejw)、X(k)x(n)100x(n)X(ejw)X(k)1-2100X(k)的密度，截?cái)嗪瘮?shù)的頻譜混疊作用沒(méi)有改變，這X(k)w1=0.48πw2=0.52π這兩個(gè)頻率分量。x(n)100N=100X(ejw)、X(k)MATLAB主要程序如下：N=100;n=0:1:N-1;xn=cos(0.48*pi*n)+cos(0.52*pi*n);Xk=dft(xn,N);subplot(3,1,1)stem(n,xn,'.k');9title('時(shí)域序列圖xn');xlabel('n');axis([0,100,-2.5,2.5]);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-j*n'*w);subplot(3,1,2);plot(w/pi,abs(Xw));title('幅頻特性曲線X(ejw)');xlabel('w');axis([0,1,0,50]);subplot(3,1,3);k1=0:1:49;w1=2*pi/100*k1;stem(w1/pi,abs(Xk(1:1:50)),'.k');title('頻域序列圖Xk');xlabel('頻率（單位：pi）');axis([0,1,0,50]);100x(n)x(n)、X(ejw)、X(k)1-3所示圖1-3100點(diǎn)x(n)對(duì)應(yīng)的x(n)、X(ejw)、X(k)10由圖可見(jiàn)，截?cái)嗪瘮?shù)的加寬且為周期序列的整數(shù)倍，改變了頻譜混疊作用，X(k)w1=0.48πw2=0.52π這兩個(gè)頻率分量。Tp來(lái)提高物理分辨率可以得到分辨率譜。Ⅵ.改進(jìn)及建議1-11-21-30<=n<=9時(shí)，從相應(yīng)的圖中幾乎x(n)90N=100DFT，從X（k）圖中可以看出，這時(shí)的譜線相當(dāng)密，故稱為高密度譜線圖，但是x(n)N=100的序X(k)Nx(n)128N=128X(ejw)、X(k)。MATLAB主要程序如下：N=128;n=0:1:N-1;xn=cos(0.48*pi*n)+cos(0.52*pi*n);Xk=dft(xn,N);subplot(3,1,1)stem(n,xn,'.k');title('xn');xlabel('n');axis([0,100,-2.5,2.5]);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-j*n'*w);subplot(3,1,2);plot(w/pi,abs(Xw));title('幅頻特性曲線X(ejw)');xlabel('w');axis([0,1,0,70]);subplot(3,1,3);k1=0:1:63;w1=2*pi/128*k1;11stem(w1/pi,abs(Xk(1:1:64)),'.k');title('頻域序列圖Xk');xlabel('頻率（單位：pi）');axis([0,1,0,70]);128x(n)x(n)、X(ejw)、X(k)1-4所示圖1-4128點(diǎn)x(n)對(duì)應(yīng)的x(n)、X(ejw)、X(k)128x(n)x(n)、X(ejw)、X(k)1-4所示。由圖可見(jiàn)，截?cái)嗪瘮?shù)雖然進(jìn)一步加寬，但不是周期序列的整數(shù)倍，所以盡管X(k)能正確分辨w1=0.48π、w2=0.52π這兩個(gè)頻率分量，但還呈現(xiàn)頻譜泄露。分辨率提高了，但仍出現(xiàn)了頻譜泄露現(xiàn)象，故要求N取值為周期序列的整數(shù)倍。取x(n)的前150點(diǎn)數(shù)據(jù)，求N=150點(diǎn)的X(ejw)、X(k)。N=150;n=0:1:N-1;12xn=cos(0.48*pi*n)+cos(0.52*pi*n);Xk=dft(xn,N);subplot(3,1,1)stem(n,xn,'.k');title('時(shí)域序列圖xn');xlabel('n');axis([0,100,-2.5,2.5]);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-j*n'*w);subplot(3,1,2);plot(w/pi,abs(Xw));title('幅頻特性曲線X(ejw)');xlabel('w');axis([0,1,0,80]);subplot(3,1,3);k1=0:1:74;w1=2*pi/150*k1;stem(w1/pi,abs(X