數(shù)字信號除處理課程設(shè)計-DFT在信號頻譜分析中的應(yīng)用

1、湘南學(xué)院課程設(shè)計課程名稱 數(shù)字信號處理 系 別： 計算機科學(xué)系 專業(yè)班級： 通信一班 學(xué) 號： 04 06 02 03 10 36 姓 名： 胡金霞 、肖雅青 、許芬、李真真 、薛明 、蔣小松 題 目： DFT在信號頻譜分析中的應(yīng)用 完成日期： 2021年 12月 23日 指導(dǎo)老師： 樊洪斌 2010年 12月 23日 目 錄1、設(shè)計題目·····················

2、3;·····················32、設(shè)計目的···························

3、83;···············33、設(shè)計原理·································&#

4、183;·········34、實現(xiàn)方法·······································&

5、#183;···35、設(shè)計內(nèi)容及結(jié)果·····································66、改良建議······&#

6、183;···································127、思考題及解答············

7、3;·························158、設(shè)計體會·······················&#

8、183;··················159、參考文獻(xiàn)······························

9、············16.設(shè)計題目 DFT在信號頻譜分析中的應(yīng)用.設(shè)計目的 掌握離散傅里葉變換的有關(guān)性質(zhì)，利用Matlab實現(xiàn)DFT變換。了解DFT應(yīng)用，用DFT對序列進行頻譜分析，了解DFT算法存在的問題及改良方法。學(xué)習(xí)并掌握FFT的應(yīng)用。.設(shè)計原理所謂信號的頻譜分析就是計算信號的傅里葉變換。連續(xù)信號與系統(tǒng)的傅里葉分析顯然不便于直接用計算機進行計算，使其應(yīng)用受到限制，而DFT是一種時域和頻域均離散化的變換，適合數(shù)值運算，成為分析離散信號和系統(tǒng)的有力工具。 工程實際中，經(jīng)

10、常遇到的連續(xù)信號Xa(t),其頻譜函數(shù)Xa(jW)也是連續(xù)函數(shù)。數(shù)字計算機難于處理，因而我們采用DFT來對連續(xù)時間信號的傅里葉變換進行逼近，進而分析連續(xù)時間信號的頻譜。.實現(xiàn)方法 離散傅里葉變換是有限長序列的傅里葉變換，它相當(dāng)于把信號的傅里葉變換進行等頻率間隔采樣，并且有限長序列的離散傅里葉變換和周期序列的離散傅里葉級數(shù)本質(zhì)是一樣的?？焖俑道锶~變換FFT并不是一種新的變換，它是離散傅里葉變換的一種快速算法，并且主要是基于這樣的思路而開展起來的：1把長度為N的序列的DFT逐次分解成長度較短的序列的DFT來計算。2利用WN(nk)的周期性和對稱性，在DFT運算中適當(dāng)?shù)姆诸悾蕴岣哌\算速度。對稱性，

11、;周期性，r為任意整數(shù) 離散傅里葉變換的推導(dǎo)：離散傅里葉級數(shù)定義為 1-1將上式兩端乘以并對n在0N-1求和可得 因為 所以 這樣用k代替m得1-2令那么1-2成為DFS 1-31-1成為IDFS 1-4式1-3、1-4式構(gòu)成周期序列傅里葉級數(shù)變換關(guān)系。其中都是周期為N的周期序列，DFS·表示離散傅里葉級數(shù)正變換，IDFS·表示離散傅里葉級數(shù)反變換。習(xí)慣上，對于長為N的周期序列，把0nN-1區(qū)間稱為主值區(qū)，把稱為的主值序列，同樣也稱為的主值序列。由于，對于周期序列僅有N個獨立樣值，對于任何一個周期進行研究就可以得到它的全部信息。在主值區(qū)研究與是等價的，因此在主值區(qū)計算DFS

12、和DFT是相等的，所以DFT計算公式形式與DFS根本相同。其關(guān)系為 所以離散傅里葉正變換 0kN-1離散傅里葉變換DFT定義：設(shè)有限長序列x (n) 長為N0nN-1，其離散傅里葉變換是一個長為N的頻率有限長序列0kN-1，其正變換為 0kN-1 離散傅里葉變換的實質(zhì)是：把有限長序列當(dāng)做周期序列的主值序列進行DFS變換，x(n)、X(k)的長度均為N，都是N個獨立值，因此二者具有的信息量是相等的。x(n)可以唯一確定X(k),X(k)可以唯一確定x(n)。雖然離散傅里葉變換是兩個有限長序列之間的變化，但它們是利用DFS關(guān)系推導(dǎo)出來的，因而隱含著周期性。構(gòu)造離散傅里葉變換的Matlab實現(xiàn)程序如

13、下： functionXk=dft(xn,N) n=0:1:N-1; k=n; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.nk; Xk=xn*WNnk快速傅里葉變換FFT并不是與DFT不同的另外一種變換，而是為了減少DFT計算次數(shù)的一種快速有效的算法 共軛對稱性：設(shè)有限長序列的長度為N，以N為周期的周期延拓列為 周期序列的共軛對稱分量和共軛反對稱分量分別為 (1-5) (1-6)同樣可以證明，它們滿足 1-7 1-8 那么有限長序列的圓周共軛對稱分量和圓周共軛反對稱分量分別定義為： 1-9 1-10由于滿足 故 1-11顯然，長度為N的有限長序列可以分解為

14、圓周共軛對稱分量和圓周共軛反對稱分量之和，和的長度皆為N。利用有限長序列與周期序列的共軛對稱分量和反對稱分量的關(guān)系式1-9和式1-10，以及式1-11可以推導(dǎo)出DFT的一系列的對稱性質(zhì)（1） DFT 式中表示的共軛復(fù)序列。證明：DFT 又因為 所以DFT（2） 復(fù)序列實部的DFT等于DFT的圓周共軛對稱局部，即 DFT證明：DFTDFT=DFT+DFT=利用DFT的對稱性可求得的DFT:設(shè) 那么DFT因為 所以DFTDFT=.設(shè)計內(nèi)容及結(jié)果1. 用MATLAB語言編寫計算序列x(n)的N點DFT的m函數(shù)文件。并與MATLAB中的內(nèi)部函數(shù)文件作比擬。 functionXk=dft(xn,N) n

15、=0:1:N-1; k=n; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.nk; Xk=xn*WNnkfunction varargout = fft(varargin)if nargout = 0 builtin('fft', varargin:);else varargout1:nargout = builtin('fft', varargin:);end運算量估計：對于N=點序列進行時間抽選奇偶分解FFT計算，需分M級，每級計算N/2個蝶。每一級需N/2次復(fù)乘、N次復(fù)加，因此總共需要進行：復(fù)乘： 復(fù)加：直接計算N點的DF

16、T，需要次復(fù)乘、N(N-1)次復(fù)加。N值越大，時間抽選奇偶分解FFT算法越優(yōu)越。例如當(dāng)N=2048點時，時間抽選奇偶分解FFT算法比直接計算DFT速度快300多倍可以用一下Matlab程序比擬DFT和FFT的運算時間N=2048;M=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;運行結(jié)果：Time_fft =Time_dft =由此可見FFT算法比直接計算DFT速度快得多2. 對離散確定信號 作如下譜分析：截取使

17、成為有限長序列N()，(長度N自己選)寫程序計算出的N點DFT ,畫出時域序列圖xnn和相應(yīng)的幅頻圖。解：1求x(n)的前10點數(shù)據(jù)對應(yīng)的X(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('時域序列圖xn');xlabel('n');axis(0,10,-2.5,2.5);w=2*pi*(0:1:2047)/2048;Xw=xn*exp(-j*n'*w

18、);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('頻率單位：pi');axis(0,1,0,10);x(n)的前10點數(shù)據(jù)對應(yīng)的x(n)、X(ejw)、X(k)如圖1-1所示。 圖1-1 x(n)的前10點數(shù)據(jù)對應(yīng)的x(n)、X

19、(ejw)、X(k)這兩個頻率分量。2將x(n)補零至100點，求N=100點的X(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('時域序列圖x1');xlabel('n');axis(0,100,-2.5,2.5);w=2*pi*(0:2047)/2048;X1=x1*exp(-j*n1'*w

20、);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)補零至100點對應(yīng)的x(n)、X(ejw)、X(k)如圖1-2所示。圖1-

21、2 x(n)補零至100點對應(yīng)的x(n)、X(ejw)、X(k)這兩個頻率分量。這說明，補零僅僅是提高了計算分辨率，得到的是高密度頻譜，而得不到高分辨率譜。3求x(n)的前100點數(shù)據(jù)，求N=100點的X(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');title('時域序列圖xn');xlabel('n');axis(0,100,-2.5,2.5);w=2*pi*(

22、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);100點x(n)對應(yīng)的x(n)、X(

23、ejw)、X(k)如圖1-3所示圖1-3 100點x(n)對應(yīng)的x(n)、X(ejw)、X(k)這兩個頻率分量。這說明通過增加數(shù)據(jù)的記錄長度Tp來提高物理分辨率可以得到分辨率譜。.改良及建議由圖1-1、圖1-2、圖1-3可以看出。當(dāng)取0<=n<=9時，從相應(yīng)的圖中幾乎無法看出有關(guān)信號頻譜的信號；將x(n)補90個零點后作N=100點的DFT，從相應(yīng)的Xk圖中可以看出，這時的譜線相當(dāng)密，故稱為高密度譜線圖，但是從中很難看出信號的頻譜局部；對x(n)加長取樣數(shù)據(jù)，得到長度為N=100的序列，此時相應(yīng)的X(k)圖中可以清晰地看到信號的頻譜成分，這稱為高分辨頻譜。為了得到更高的分辨率，增加

24、N點的取值進行改良取x(n)的前128點數(shù)據(jù)，求N=128點的X(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)

25、;title('幅頻特性曲線X(ejw)');xlabel('w');axis(0,1,0,70);subplot(3,1,3);k1=0:1:63;w1=2*pi/128*k1;stem(w1/pi,abs(Xk(1:1:64),'.k');title('頻域序列圖Xk');xlabel('頻率單位：pi');axis(0,1,0,70);128點x(n)對應(yīng)的x(n)、X(ejw)、X(k)如圖1-4所示 圖1-4 128點x(n)對應(yīng)的x(n)、X(ejw)、X(k) 128點x(n)的數(shù)據(jù)對應(yīng)的x(n)、X

26、(ejw)、X(k)如圖1-4所示。這兩個頻率分量，但還呈現(xiàn)頻譜泄露。分辨率提高了，但仍出現(xiàn)了頻譜泄露現(xiàn)象，故要求N取值為周期序列的整數(shù)倍。取x(n)的前150點數(shù)據(jù)，求N=150點的X(ejw)、X(k)。N=150;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*e

27、xp(-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(Xk(1:1:75),'.k');title('頻域序列圖Xk');xlabel('頻率單位：pi');axis(0,1,0,80);15點x(n)對應(yīng)的x(n)、X(ejw)、X(k)如圖1-5所示 圖1-5 1

28、50點x(n)對應(yīng)的x(n)、X(ejw)、X(k) 150點x(n)的數(shù)據(jù)對應(yīng)的x(n)、X(ejw)、X(k)如圖1-4所示。這兩個頻率分量，到達(dá)了很好的效果。.思考題及解答(1)比照設(shè)計內(nèi)容2中(1)(2) (3)的圖，說明補零DFT的作用。 由圖1-1、1-2、1-3可知DFT是有限長序列的頻譜等間隔采樣所得到的樣本值，這就相當(dāng)于透過一個柵欄去觀察原來信號的頻譜，因此必然有一些地方被柵欄所遮擋，這些被遮擋的局部就是未被采樣到的局部，這種現(xiàn)象稱為柵欄效應(yīng)。如下列圖由于柵欄效應(yīng)總是存在的，因而可能會使信號頻率中某些較大的頻率分量由于被“遮擋二無法得到反映。此時，通常在有限長序列的尾部增補假

29、設(shè)干個零值，介意改變原序列的長度。這樣加長的序列作DFT時，由于點數(shù)增加就相當(dāng)于調(diào)整了原來柵欄的間隙即間隔頻率，可以使得原來的不到反映的那些較大的頻率分量落在采樣點上而得到反映。但要注意，由于柵欄效應(yīng)，使得被分析的頻譜變得較為稀疏，為此，在采樣樣本序列x(n)后補零，在數(shù)據(jù)長度不變的情況下，可以改變頻譜的頻率取樣密度，得到高密度頻譜。(2)解釋設(shè)計內(nèi)容3中圖和圖有什么區(qū)別？補零DFT能否提高信號的頻譜分辨率，說明提高頻譜密度、頻譜分辨率的措施各是什么？ 圖1-2在圖1-1的根底上提高了頻譜密度和計算分辨率，但是頻譜的包絡(luò)沒有發(fā)生變化所以圖和圖沒有區(qū)別。 補零作用不能提高信號的頻譜分辨率，因在x

30、(n)后面補零并沒有增加新的信息量，改善的僅是柵欄效應(yīng)，所以補零是不能提高頻率分辨率的，即得不到高分辨率譜。這說明，補零僅僅是提高了計算分辨率，得到的是高密度頻譜，而要得到高分辨率譜，那么要通過增加數(shù)據(jù)的記錄長度來提高物理分辨率。.設(shè)計體會及心得 計算機是進行數(shù)字信號處理的主要工具，計算機只能處理有限長序列，這就決定了有限長序列處理在數(shù)字信號處理中的重要地位。離散傅里葉變換建立了有限長序列與其近似頻譜之間的聯(lián)系，在理論上具有重要意義。離散傅里葉變換DFT在數(shù)字通信、語音處理、圖像處理、譜估計、仿真、系統(tǒng)分析、雷達(dá)、光學(xué)影像、地震等各個領(lǐng)域得到廣泛應(yīng)用，但是這都是以卷積和相關(guān)運算，對連續(xù)信號和序

31、列進行譜分析為根底的。 通過該課程設(shè)計，我們受益匪淺，對DFT在進行頻譜的分析上有了根深刻的理解和掌握。DFT實現(xiàn)了頻域采樣，同時DFT存在快速算法FFT，所以在實際應(yīng)用中，可以利用計算機，用DFT來逼近連續(xù)時間信號的傅里葉變換，進而分析連續(xù)時間信號頻譜。同時知道了補零點的作用，其僅僅是提高了計算分辨率，得到的是高密度頻譜，并不能得到高分辨率譜，要提高頻率分辨率，那么要通過增加數(shù)據(jù)記錄長度來提高物理分辨率。 在編程實現(xiàn)中，遇到了一些問題，為此我們翻閱一些了參考書，并通過討論一一解決。期間我們不僅學(xué)到了許多課本上的知識，還有課本以外的內(nèi)容，學(xué)到了許多課本上所沒提到的東西，這些東西都讓我們耳目一新

32、，開闊了視野，拓寬了知識面。從以前僅僅掌握離散傅里葉變換的概念，到現(xiàn)在漸漸領(lǐng)悟到離散傅里葉變換的一些實際應(yīng)用，更明白它在實際設(shè)計中的作用，從理論到實踐的逐步過渡，增了動手能力。知道了到團隊精神的重要性，大家互相討論，分工合作，享受了合作的樂趣。.參考文獻(xiàn)1余成波，陶紅艷.數(shù)字信號處理及MATLAB實現(xiàn)第二版.北京：清華大學(xué)出版社.2021.1.p98-123.2王艷芬，王剛.數(shù)字信號處理原理及實現(xiàn).北京：清華大學(xué)出版社.2021.3.p96-105.3從玉良，王宏志.數(shù)字信號處理原理及其MATLAB實現(xiàn)第2版.北京：電子工業(yè)出版社.2021.7.p63-100.目錄一、設(shè)計題目二、設(shè)計內(nèi)容題目

33、一：用窗函數(shù)法設(shè)計FIR數(shù)字低通濾波器一設(shè)計目的二設(shè)計任務(wù)與要求三設(shè)計原理四實現(xiàn)方法五設(shè)計結(jié)果六思考題三、設(shè)計心得體會四、參考文獻(xiàn)一、設(shè)計題目FIR數(shù)字低通濾波器二、設(shè)計內(nèi)容題目一：用窗函數(shù)法設(shè)計FIR數(shù)字低通濾波器一設(shè)計目的1. 熟悉設(shè)計線性相位數(shù)字濾波器的一般步驟。2. 掌握用窗函數(shù)法設(shè)計FIR數(shù)字濾波器的原理和方法。3. 熟悉各種窗函數(shù)的作用以及各種窗函數(shù)對濾波特性的影響。4. 學(xué)會根據(jù)指標(biāo)要求選取適宜的窗函數(shù)。二設(shè)計任務(wù)與要求1. 學(xué)會計算濾波器各項性能指標(biāo)及如何來滿足給定的指標(biāo)要求。2. 用MATLAB語言編程實現(xiàn)給定指標(biāo)要求的濾波器的設(shè)計。3. 熟悉MATLAB語言，獨立編寫程序。

34、4. 設(shè)計低通FIR濾波器的指標(biāo)：通帶最大波動 ，阻帶最小衰減 ，三設(shè)計原理1.FIR數(shù)字濾波器的獨特優(yōu)點是容易得到嚴(yán)格的線性相位，其極點都位于原點，所以總是穩(wěn)定的。FIR主要采用非遞歸結(jié)構(gòu)，而且可以利用FFT計算以提高運算效率，因此，F(xiàn)IR濾波器得到了越來越廣泛的應(yīng)用，它的設(shè)計方法主要有窗函數(shù)法和頻率響應(yīng)法。(1)給出希望設(shè)計的濾波器的頻率響應(yīng)函數(shù)H(ejw)(2)根據(jù)允許的過渡帶寬度和阻帶衰減，來確定N 值窗函數(shù)旁瓣峰值幅度/db主瓣過渡帶寬度最小阻帶衰減/db矩形窗-134pi/N-21Hamming窗-418pi/N-53Hanning窗-318pi/N-44Kaiser窗-5710p

35、i/N-80(3)求出Hd(n). (4)將Hd(n)與窗函數(shù)相乘得到fir數(shù)字濾波器的單位取樣響應(yīng)h(n);h(n)=hd(n)*w(n)(5)計算fir數(shù)字濾波器的頻率響應(yīng)，并檢驗是否到達(dá)所要求的指標(biāo)(6)由Hejw計算幅頻響應(yīng)H(W)和相頻響應(yīng)(w)窗函數(shù)的傅式變換W(ej)的主瓣決定了H(ej)過渡帶寬。W(ej)的旁瓣大小和多少決定了H(ej)在通帶和阻帶范圍內(nèi)波動幅度，常用的幾種窗函數(shù)有：· 矩形窗 w(n)=RN(n)； · Hanning窗 ； · Hamming窗 ； · Blackmen窗 ； · Kaiser窗 。 式中I

36、o(x)為零階貝塞爾函數(shù)。四實現(xiàn)方法用Matlab設(shè)計觀察窗函數(shù)相關(guān)波形1、 編寫計算N階差分方程所描述系統(tǒng)頻響函數(shù) 。functionH=fr(b,a,w);m=0:length(b)-1;l=0:length(a)-1;num=b*exp(-j*m'*w);den=a*exp(-j*l'*w);H=num./den2、 編寫計算理想低通濾波器單位抽樣響應(yīng)hd(n)的m函數(shù)文件idealm。functionhd=ideal(wc,N)q=(N-1)/2;n=0:N-1;m=n-q+eps;hd=sin(wc*m)./(pi*m)3、 編寫m程序文件，通過調(diào)用idealm和frm文件，計算你設(shè)計的實際低通FIR濾波器的單位抽樣響應(yīng)h(n)和頻率響應(yīng)，打印在頻率區(qū)間O，上的幅頻響應(yīng)特性曲線，幅度用分貝表示。wp=0.2*pi;ws=0.3*pi;width=ws-wpN=ceil(6