1、1 引言 12基于MATLAB勺FFT算法實現 21.1 系統總體流程圖 21.2 FFT 運算規律及編程思想 31.2.1 語音信號的采集 31.2.2 DIT-FFT 算法的基本原理 31.2.3 DIT-FFT 算法的運算規律及編程思想 53 Matlab 程序實現 104 系統人機對話界面 134.1 GUI 簡介 134.2 界面設計 134.3 運行調試 145 心得體會 16參考文獻 17附錄I 18附錄H 211 引言MATLA是矩陣實驗室(Matrix Laboratory )的簡稱，是美國 MathWorks公 司出品的商數學軟件，用于算法開發、 數據可視化、 數據分析以及

2、數值計算的高 級技術計算語言和交互式環境，主要包括MATLAB口 Simulink兩大部分。MATLAB勺應用范圍非常廣，包括信號和圖像處理、通訊、控制系統設計、 測試和測量、 財務建模和分析以及計算生物學等眾多應用領域。 附加的工具箱 (單 獨提供的專用MATLAB函數集)擴展了 MATLAB環境，以解決這些應用領域內特 定類型的問題。它以矩陣運算為基礎，把計算、可視化、程序設計融合在一個簡 單易用的交互式工作環境中， 是一款數據分析和處理功能都非常強大的工程適用 軟件。 它可以將聲音文件變換為離散的數據文件， 然后利用其強大的矩陣運算能 力處理數據，如數據濾波、傅立葉變換、時域和頻域分析、

3、聲音回放以及各種圖 的呈現等， 它的信號處理與分析工具箱位語音信號分析提供了十分豐富的功能函 數， 利用這些功能函數可以快捷而又方便的完成語音信號的處理和分析以及信號 的可視化。數字信號處理是MATLABt要應用的領域之一。對于有限長序列x(n),若要求其N點的傅里葉變換(DFT需要經過M次復 數乘法運算和N*(N-1)次復數加法運算。隨著N的增加，運算量將急劇增加，而 在實際問題中，N往往是較大的，如當N=1024時，完成復數乘法和復數加法的 次數分別為百萬以上，無論是用通用計算機還是用 DSPE片，都需要消耗大量的 時間和機器內存，不能滿足實時的要求。因此， DFT的這種運算只能進行理論上

4、 的計算，不適合對實時處理要求高的場合。因此，研究作為 DSP 的快速算法的 FFT是相當必要的，快速傅里葉變換(FFT)是為提高DFT運算速度而采用的一 種算法， 快速算法的種類很多，而且目前仍在改進和提高， 它是根據離散傅里葉 變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的?；?本學期所學的 DIT-FFT 的運算規律和編程思想以及Matlab 的學習和使用，本課設要求在Matlab 環境下編寫基2 DIT-FFT 算法實現對離散信號的快速傅里葉變換，再與Matlab 軟件自帶的FFT 函數實現對離散信號的傅里葉變換進行比較，如果得到的頻譜相同，那么我們編寫的程序就是正

5、確的。其中離散信號是通過PC自帶的錄音機錄制一段 wav語音信號，用Matlab采樣得到離散序列x1。如果 有能力可以選做系統人機對話界面。 用GUI界面完成人機交互方便使用的。 本課 程設計主要是對數字信號的分析。2基于MATLAB勺FFT算法實現2.1 系統總體流程圖本設計要求錄制一段個人自己的語音信號，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；在Matlab環境下編寫基2 DIT-FFT算法;利用自己編寫的算法對已采集的語音信號進行頻譜分析，并畫出語音信號的時域與頻譜圖，并與Matlab數字信號處理工具箱中的fft函數進行對比研究，驗證 自編算法的正確性。所以得到系統

6、總體流程圖如圖1所示。語音信號采集編寫fft程序，畫出信號頻譜圖完成信號時域圖兀 成 信 號 頻 率 響 應實 現 輸 入 信 號 的 倒 序實現一級中不同種蝶形算運實 現級 中 相 同 種 蝶 形 運Matlab自帶的FFT比較圖1系統總體流程圖2.2 FFT運算規律及編程思想2.2.1 語音信號的采集利用PM自帶的錄音機，錄制一段語音信號，保存格式為 wave的文件，并 將其保存在電腦中。在 MATLAB 中， fn=input(' Enter WAV filename:','s'); x,fs,nb=wavread(fn,n1 n2);用于讀取語音，采樣值

7、放在向量 x中，fs表示采樣頻率(Hz) , nb表示采樣位數。n1 n2表示讀取從n1點到n2點的值(若 只有一個n的點則表示讀取前n點的采樣值)。sound(x,fs,nb); 用于對聲音的 回放。向量x則就代表了一個信號(也即一個復雜的“函數表達式”)也就是說可以像處理一個信號表達式一樣處理這個聲音信號。采集到語音信號之后，需要對語音信號進行分析，如語音信號的時域分析、 頻譜分析、譜圖分析。2.2.2 DIT-FFT算法的基本原理快速傅里葉變換(FFT)是為提高DFT運算速度而采用的一種算法。 對一個有限長度序列x(n)的N點的DFT為：N 1X k xn 01 Nx n 一WNnk ,

8、 k 0, 1, 2,., N 1k Wn nk,n 0, 1, 2,., N 1所以，要求N點的DFT,需要2次的復數乘法運算，N*(N-1)次復數乘法運算算。隨著N的增加，運算量將急劇增加，而在實際問題中，N往往是較大的，如當N=1024時，完成復數乘法和復數加法的次數分別為百萬以上，無論是用通用 計算機還是用DSP芯片，都需要消耗大量的時間，不能滿足實時的要求不適合于對實時處理要求高的場合。為了能實時處理DFT,要想減少DFT的運算量可以有兩個途徑：第一是降N, N的值減小了，運算量就減少了；第二是利用旋轉 因子的周期性，對稱性和可約性。利用這兩個途徑實現DFT的快速傅里葉變換(FFT)

9、, FFT算法基本上可分為按時間抽取的 FFT算法(DIT-FFT)和按頻率抽 設序列x(n)的長度為N,且N滿足N=2,M為正整數。若N不能滿足上述關系， 可以將序列x(n)補零實現。按時間抽取基 2-FFT算法的基本思路是將 N點序列 按時間下標的奇偶分為兩個 N/2點序列，計算這兩個N/2點序列的N/2點DFI 計算量可減小約一半；每一個N/2點序列按照同樣的劃分原則，可以劃分為兩個 N/4點序列，最后，將原序列劃分為多個 2點序列，將計算量大大降低。取的FFT算法(DIF-FFT)旋轉因子的性質：(1)周期性(2)共腕對稱性(3)可約性本次課設要求用用基WW;n WWN啊*WNWN)n

10、)* W；n wN k)n* W；(n)*Wkn WCWkn WNk肝2的按時間抽取的FFT算法(DIT-FFT)實現FFT功能,按時間下標的奇偶將N點x(n)分別抽取組成兩個N/2點序列，分別記為x1(n) 和x2(n),將x(n)的DFTW化為x1(n)和x2(n)的DFT的計算。r 0,1,x(2r) xi(r)x(2r 1)x2(r)n 0,2,4kx nx n WNnkn 1,3,5N12x 2rW2rkr 0,1N12x 2r 1 Wn2r 1 kr 0,12rkxr Wnr 0,1N 122r 1 kx2r Wnr 0,1利用旋轉因子的可約性，即:W；*X kj2_2rk e N

11、N 12 .rkx r Wn r 02.2jNrkTN2 k WnrWn：x20r WN"X 1(k) W；X2(k), 0 k N 1用蝶形運算可表示為如圖2所示:(k) 人 一_ Y X (#)=國 + 際X" (k)M+幻=X") -際與圖2 DIT-FFT蝶形運算流圖符號以此類推，還可以把x1(n)和x2(n)按n值得奇偶分為兩個序列，這樣就達到了降N得目的，從而減少了運算量FFT對DFT的數學運算量改進:直接采用DFT進行計算，運算量為 M次復數乘法和N*(N-1)次復數乘法。當采用M次FFT時，由N=皆求得M=logN,運算流圖有M級蝶形，每一級都 由

12、N/2個蝶形運算構成，這樣每一級蝶形運算都需要 N/2次復數乘法和N次復數 加法。M級運算共需要復數乘法次數為 C=N/2*M,復數加法次數為C=N*M當N值較大時，FFT減少運算量的特點表現的越明顯。2.2.3 DIT-FFT 算法的運算規律及編程思想為了編寫 DIT-FFT 算法的運算程序， 首先要分析其運算規律， 總結編程思想 并繪出程序框圖。1. 原位計算對N 2M點的FFT共進行M級運算，每級由N/2個蝶形運算組成。在同一級 中，每個蝶的輸入數據只對本蝶有用，且輸出節點與輸入節點在同一水平線上，這就意味著每算完一個蝶后， 所得數據可立即存入原輸入數據所占用的數組元素 ( 存儲單元 )

13、 ，這種原位( 址 ) 計算的方法可節省大量內存。2. 蝶形運算實現FFT運算的核心是蝶形運算，找出蝶形運算的規律是編程的基礎。蝶形 運算是分級進行的； 每級的蝶形運算可以按旋轉因子的指數大小排序進行； 如果 指數大小一樣則可從上往下依次蝶算。對 N 2M點的FFT共有M級運算，用L表 示從左到右的運算級數(L=1, 2,M)。第L級共有B 2L1個不同指數的旋轉 因子，用R表示這些不同指數旋轉因子從上到下的順序(R=0, 1,，B-1) 0第R個旋轉因子的指數P 2MLR,旋轉因子指數為P的第一個蝶的第一節點標號k 從R開始，由于本級中旋轉因子指數相同的蝶共有2ML個，且這些蝶的相鄰間距為

14、2L ， 故旋 轉 因子指數為 P 的最后 一個蝶 的 第 一節 點標號 k 為 ：MLLL(21) 2 R N 2 R ，本級中各蝶的第二個節點與第一個節點都相距B 點。應用原位計算，蝶形運算可表示成如下形式：AL(J)= AL-1(J)+ AL-1 (J+B)* WNPAL(J+B)= AL-1(J)- AL-1(J+B)* WNP總結上述運算規律， 可采用如下運算方法進行DIT-FFT 運算。 首先讀入數據，根據數據長度確定運算級數 M,運算總點數N 2M ,不足補0處理。然后對讀入 數據進行數據倒序操作。數據倒序后從第1級開始逐級進行，共進行 M級運算。在進行第 L 級運算時，先算出該

15、級不同旋轉因子的個數B 2L1 (也是該級中各個蝶形運算兩輸入數據的間距)，再從R=0開始按序計算，直到R=B-1結束。每個 R對應的旋轉因子指數P 2MLR，旋轉因子指數相同的蝶從上往下依次逐個運算， 各個蝶的第一節點標號k都是從R開始，以2L為步長，到N 2L R(可簡取極值N-2）結束?？紤]到蝶形運算有兩個輸出，且都要用到本級的兩個輸入數據，故第 一個輸出計算完畢后，輸出數據不能立即存入輸入地址，要等到第二個輸出計算 調用輸入數據完畢后才能覆蓋。這樣數據倒序后的運算可用三重循環程序實現 整個蝶形運算流程圖如圖3所示。3. 序列倒序為了保證運算輸出的X(k)按順序排列，要求序列x(n)倒序

16、輸入，即在運算 前要先對輸入的序列進行位序顛倒。如果總點數為n 2M的x(n)的順序數是用M位二進制數表示，則倒序數只需將順序數的二進制位倒置即可，按照這一規律用硬件電路和匯編語言很容易產生倒序數。但用MATLA舒高級語言實現倒序時，直接倒置二進制數位的方法不可取，還須找出產生倒序的十進制規律。 將十進制順序數用I表示，與之對應的二進制數用 舊表示。十進制倒序數用J表示，與 之對應的二進制數用JB表示。JB是舊的位倒置結果，十進制順序數I增加1, 相當于 舊最低位加1且逢2向高位進1,即相當于JB最高位加1且逢2向低位 進10JB的變化規律反映到J的變化分二種情況：如果JB的最高位是0(J n

17、/2), 則直接由加1(J J n/2)得到下一個倒序值；如果JB的最高位是1(J n/2),則 要先將最高位變0(J J n/2),再在次高位加1(J J n/4)。但次高位加1時， 同樣要判斷0、1值，如果是0 (J N/4),則直接加1(J J N/4),否則要先將 次高位變0(J J n/4),再判斷下一位。依此類推，直到完成最高位加1,逢2向右進位的運算。利用這一算法可按順序數I的遞增順序，依次求得與之對應的 倒序數Jo為了節省內存，數據倒序可原址進行，當I = J時不需要交換，當I* J時需要交換數據。另外，為了避免再次調換前面已經調換過的一對數據，只對 I <J的情況進行數

18、據交換即可實現數據倒序操作。圖 3中數據倒序的程序流程圖 如圖4所示。例如，N=8時，序列倒序結果如表1所示。表1碼位倒序(N=8)00 0 00 0 0010 0 11 0 0420 1 00 1 0230 1 11 1 0641 0 00 0 1151 0 11 0 1561 1 00 1 1371 1 11 1 17圖4倒序流程圖3 Matlab程序實現用ringin.wav作為例子，運行調試程序。1.程序運行開始時，要求輸入采樣頻率fs (fs=1024),所需要變換的起點N1 (N1=100O和終點N2 (N2=509§以及要采樣白語音文件(ringin.wav );其 中

19、N1和N2的差值必須在語音信號長度范圍內， 但不能太小，否則聽不到較為清 晰的語音。具輸入窗口顯示如圖 5所示。圖5程序運行開始的輸入窗口將fs=1024, N1=100Q N2=5095和語音信號ringin.wav 輸入后得到采樣后 的語音信號x1的時域波形，如圖6所示，其頻率響應特性如圖7所示，Matlab 自帶的FFT函數實現的x1的頻譜與編寫的FFT程序實現的x1的頻譜的比較如圖 8所示。語音信號時域波形0.80.60.40.20-0.2-0.4-0.6-0.8 匚1LJ050010001500200025003000350040004500圖6采樣后的語音信號的時域波形60x1的頻

20、率響應特性50403020100.511.522.5頻率3.5圖7 x1的頻率響應特性Matlab自帶的FFT函數實現的x1的頻譜200150 -k 100 一-Y|50 -I ,-JhLI11 JH 10仁；. N'一 2、:-中, 一050010001500200025003000350040004500KK圖8 FFT函數實現與編寫的 FFT程序實現的比較取不同的點數進行FFT變換，經觀察，編寫FFT程序得到的語音信號的頻譜 圖與Matlab中自帶的FFT函數得到的語音信號頻譜圖總是基本一致，但是如果 輸入的N1和N2差值加1不是2的整數次幕就會有細微差別，這是因為編寫的快 速傅

21、里葉變換計算信號的 N點傅里葉變換要求N為2的整數次幕，不夠的話信 號將會被補零后運算，即參與運算的信號已經不同（差別很?。?，所以FFT運算后的結果也不盡相同，所示頻譜圖自然有細微的差別。由用MATLAB自帶FFT函數實現的頻譜圖與用 MATLA踹寫的FFT程序實現 的頻譜圖相比較，可知，兩個算法計算后的結果幾乎相同，驗證了自編算法的正 確性。如果改變輸入的值和語音信號，那么得到不一樣的時域波形，頻率響應，通 過FFT得到頻譜都會不同。需要注意的是自由輸入的語音信號不同，其長度也不同，所以根據要輸入的語音信號輸入 N1, N2的值。4系統人機對話界面4.1 GUI簡介圖形用戶界面(GUI),是

22、一種提供人機交互的工具和方法。GUI是包含圖形對象，如窗口、圖標、菜單和文本等圖文并茂的用戶界面。4.2 界面設計用MATLA圖形用戶界面開發環境設計 GUI點的一般步驟是：1 .進行界面設計。2 .設計控件屬性。3 .進行M語言編程。以本設計要求為例介紹。第一步，該選擇本圖形用戶界面需要的控件：三個文本框(Edit Text) 用來輸入N1, N2和fs。三個靜態文本(static text)用來提示我們需要人為輸入的內容。兩個推按鈕(Push button),用來運行和退出。四個軸對象(axes)用來顯示信號時域波形，頻譜圖，和兩種不同FFT實現的頻譜 圖。完成人機設計界面如圖9所示：圖9

23、完成人機設計界面第二步，設置控件屬性：雙擊組件可以設置文本框，推按鈕的屬性，如顯示大小，名稱和默認值等。第三步，編寫回調函數。組件事件的發生是通過回調函數進行工作的。控件設置完成后保存，然后運行GUI （操彳乍為ctrl+T ）,就會進入editor窗口，加入各個控件功能的函數代碼。 完成后保存即可。第四步，運行GUI。運行editor窗口的程序后，會彈出已經激活的人機對話界面。系統人機對話界面如圖10所示。圖10系統人機對話界面4.3運行調試運行GUI,已經彈出圖10的系統人機對話界面，根據提示“Input N1:"，“Input N2:"和“Input fs :&quo

24、t;分別輸入N1, N2和fs ,然后點擊“運行”推按鈕就會得 到本設計的要求。運行結果如圖11所示。如果運行結束，可以通過點擊“退出” 推按鈕退出該人機對話界面，返回 Matlab。匚慶圖11運行結果5 心得體會本次實習的主要內容是通過用 Matlab 實現FFT 的設計，可以實現對一段自己錄制的語音信號進行分析，并畫出采樣信號的時域與頻域圖。把自己編寫的FFT 算法與 Matlab 自帶 FFT 算法進行比較。程序運行調試時，自己選擇輸入要采樣的語音信號， 采樣頻率以及要變換的范圍， 可以實現對不同信號的信號采樣和進行不同點的FFT運算。在之前數字信號處理的學習以及完成實驗的過程中，已經使

25、用過Matlab ，對其有了一些基礎的了解和認識， 通過這次的課程設計使我進一步了解了信號的產生， 采樣及頻譜分析的方法， 以及其中產生信號和繪制信號的基本命令和一些基礎編程語言。 讓我感受到只有在了解課本知識的前提下， 才能更好的應用這個工具，并且熟練的應用 Matlab 也可以很好的加深我對課程的理解，方便我的思維。這次課程設計使我了解了 Matlab 的使用方法，提高了自己的分析和動手實踐能力。同時我相信，進一步加強對MATLAB勺學習與研究對我今后的學習將會起到很大的幫助。這次的課程設計是對本學期所學知識的一次重要鞏固， 使得在課堂上掌握的知識得到了真正的運用。 在學習的過程中和同學討

26、論， 更明白了理論知識與實踐的聯系。 書到用時方恨少， 有些知識學會是一回事， 掌握是一回事， 但應用起來，確實不是那么簡單的，需要很多知識的融會貫通。程序運行調試初期， 曾經多次出現錯誤、不能產生圖形等問題，但在我翻閱資料認真改正及老師同學的幫助下基本功能還是完成了， 經過 1 個星期的上機實習，程序已得到一些完善，能完成基本的要求的功能。最后經過努力，又深入學習了圖形用戶界面（GUI）,完成了選做要求的人機對話界面。學習就是一個了解， 疑惑，進而解惑的過程，這次實習就是提供了這樣一個發現自己知識漏洞，與同學老師探討進行解惑的的機會。通過這次課程設計實習， 我更深刻的了解了 Matlab 的

27、運用， 重新復習了 FFT中的重要的序列倒序和蝶形變換的程序， 對課本上的知識有了更深的理解， 使我對數字信號處理有了系統的認知。在這里特別感謝董老師和李老師， 他們給了我們很大的發揮空間， 讓我們真正自己動手真正掌握了知識，感謝他們細心指導。也非常感謝我的同學， 他們解開了我在實習中出現的諸多知識死角，謝謝大家！參考文獻1范壽康DSP技術與DSP芯片.北京：電子工業出版社2 程佩青 . 數字信號處理教程. 北京：清華大學出版社出版， 20013 高西全 丁玉美等 . 數字信號處理. 北京：電子工業出版社， 20094李勇徐震.MATLAB輔助現代工程數字信號處理.西安電子科技大學出版社5 陳

28、杰 .Matlab 寶典 . 電子工業出版社6蘇金明張蓮花劉波.MATLAB工具箱應用，電子工業出版社附錄Ifs=input(' 輸入采樣頻率fs='); % 語音信號采樣頻率為 fsN1=input(' 輸入所需變換的起點N1=');N2=input(' 輸入所需變換的終點N2=');fn=input(' Enter WAV filename:','s'); %獲取一個 *.wav 的文件x,fs,nb=wavread(fn,N1 N2); % 讀取語音信號的數據sound(x,fs,nb); % 播放語音信號%

29、n=N2-N1+1;臉語音信號文件較大時用這兩條%x1=reshape(x,1,2*n);% 語句替換 x1=x'x1=x'y1=fft(x1);figure(1)plot(x1) % 做原始語音信號的時域圖形title(' 語音信號時域波形 ')xlabel('n');ylabel(' 幅值 ');M=nextpow2(x1); % 求x的長度對應的2的最低幕次mN=2AM;if length(x1)<Nx1=x1,zeros(1,N-length(x1);%若x的長度不是2的幕，補零到2的整數冪end%數據倒序操作J=0

30、;%給倒序數賦初值for I=0:N-1;% 按序交換數據和算倒序數if I<J;% 條件判斷及數據交換T=x1(I+1);x1(I+1)=x1(J+1);x1(J+1)=T;end% 算下一個倒序數K=N/2;while J>=K;J=J-K;K=K/2;endJ=J+K;end %x1;y=x1; %將x倒序排列作為y的初始值WN=exp(-i*2*pi/N);for L=1:MB=2AL/2;% 第L級中，每個蝶形的兩個輸入數據相距B個點，每級有Bt不同的旋轉因子p=J*2A(M-L);WNp=WNAp;for k=J+1:2AL:N %kp=k+B; %t=y(kp)*WN

31、p; % y(kp)=y(k)-t; %for J=0:B-1 % J 代表了不同的旋轉因子本次蝶形運算的跨越間隔為 2AL蝶形運算的兩個因子對應單元下標的關系蝶形運算的乘積項蝶形運算， 注意必須先進行減法運算，然后進行加法運算，否則要使用中間變量來傳遞y(k)y(k)=y(k)+t; %蝶形運算endendend%yfigure(2)x1,w1=freqz(x1,1); % 繪制原始語音信號的頻率圖plot(w1,abs(x1);title('x1 的頻率響應特性')xlabel(' 頻率');ylabel(' 幅度');figure(3)su

32、bplot(2,1,1);plot(abs(y1) %Matlab 自帶的FFT®數實現的語音信號的FF砌譜圖title('Matlab 自帶的FFT®數實現的x1的頻譜')xlabel('K');ylabel('Y1(k)');subplot(2,1,2);plot(abs(y); % 編寫的FFTS序實現的語音信號的FFT1譜圖 title(' 編寫的FFTS序實現的的x1的頻譜') xlabel('K');附錄nfunction varargout = gaoli208(varargin)

33、% GAOLI208 M-file for gaoli208.fig% GAOLI208, by itself, creates a new GAOLI208or raises the existing% singleton*.% H = GAOLI208 returns the handle to a new GAOLI208 or the handleto% the existing singleton*.% GAOLI208('CALLBACK',hObject,eventData,handles,.) calls the localfunction named CALL

34、BACKin GAOLI208.M with the given inputarguments.%GAOLI208('Property','Value',.)creates a newGAOLI208or raisesthe% existing singleton*. Starting from the left, property value pairs are% applied to the GUI before gaoli208_OpeningFunction gets called.An% unrecognized property name or in

35、valid value makes property application% stop. All inputs are passed to gaoli208_OpeningFcn via varargin.% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one% instance to run (singleton)".% See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the respons

36、e to help gaoli208% Last Modified by GUIDE v2.5 11-Jan-2010 14:14:19 % Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, .'gui_Singleton', gui_Singleton, .'gui_OpeningFcn', gaoli208_OpeningFcn, .'gui_OutputFcn', gaoli20

37、8_OutputFcn, .'gui_LayoutFcn', , .'gui_Callback', );if nargin && ischar(varargin1)gui_State.gui_Callback = str2func(varargin1);endif nargoutvarargout1:nargout = gui_mainfcn(gui_State, varargin:);elsegui_mainfcn(gui_State, varargin:);end% End initialization code - DO NOT EDIT%

38、 - Executes just before gaoli208 is made visible.function gaoli208_OpeningFcn(hObject, eventdata, handles, varargin)% This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user

39、data (see GUIDATA)% varargin command line arguments to gaoli208 (see VARARGIN)% Choose default command line output for gaoli208handles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes gaoli208 wait for user response (see UIRESUME) % uiwait(handles.figure1);% - Outp

40、uts from this function are returned to the command line. function varargout = gaoli208_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT);% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structu

41、re with handles and user data (see GUIDATA)% Get default command line output from handles structure varargout1 = handles.output;function N1_Callback(hObject, eventdata, handles)% hObject handle to N1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with

42、 handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of N1 as text% str2double(get(hObject,'String') returns contents of N1 as a double% - Executes during object creation, after setting all properties.function N1_CreateFcn(hObject, eventdata, handles)%

43、hObject handle to N1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Ba

44、ckgroundColor'),get(0,'defaultUicontrolBackgroundColor')set(hObject,'BackgroundColor','white');end function N2_Callback(hObject, eventdata, handles)% hObject handle to N2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with

45、handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of N2 as text% str2double(get(hObject,'String') returns contents of N2 as a double% - Executes during object creation, after setting all properties.function N2_CreateFcn(hObject, eventdata, handles)% h

46、Object handle to N2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Bac

47、kgroundColor'),get(0,'defaultUicontrolBackgroundColor')set(hObject,'BackgroundColor','white');end function fs_Callback(hObject, eventdata, handles)% hObject handle to fs (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with h

48、andles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of fs as text% str2double(get(hObject,'String') returns contents of fs as a double% - Executes during object creation, after setting all properties. function fs_CreateFcn(hObject, eventdata, handles)% h

49、Object handle to fs (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB%handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.ifispc&& isequal(get(hObject,'Bac

50、kgroundColor'),get(0,'defaultUicontrolBackgroundColor') set(hObject,'BackgroundColor','white');end% - Executes on button press in begin.function begin_Callback(hObject, eventdata, handles)% hObject handle to begin (see GCBO)% eventdata reserved - to be defined in a future

51、 version of MATLAB% handles structure with handles and user data (see GUIDATA)N1=str2double(get(handles.N1,'String');N2=str2double(get(handles.N2,'String');fs=str2double(get(handles.fs,'String');x,fs,bits=wavread('D:Matlabworkringin.wav',N1 N2);sound(x,fs,bits); % 播放語

52、音信號%n=N2-N1+1;臉語音信號文件較大時用這兩條%x1=reshape(x,1,2*n);% 語句替換 x1=x'x1=x'y1=fft(x1);axes(handles.tu1);plot(x1) % 做原始語音信號的時域圖形title(' 語音信號時域波形 ')xlabel('n');ylabel(' 幅值 ');M=nextpow2(x1); % 求x的長度對應的2的最低幕次mN=2AM;if length(x1)<Nx1=x1,zeros(1,N-length(x1);%若x的長度不是2的幕，補零到2的整數冪

53、end%數據倒序操作J=0;%給倒序數賦初值for I=0:N-1;% 按序交換數據和算倒序數if I<J;% 條件判斷及數據交換T=x1(I+1);x1(I+1)=x1(J+1);x1(J+1)=T;end% 算下一個倒序數K=N/2;while J>=K;J=J-K;K=K/2;endJ=J+K;end%x1;y=x1; %將x倒序排列作為y的初始值WN=exp(-i*2*pi/N);for L=1:MB=2AL/2;% 第L級中，每個蝶形的兩個輸入數據相距B個點，每級有Bt不同的旋轉因子for J=0:B-1 % J 代表了不同的旋轉因子WNp=WNAp;for k=J+1:

54、2AL:N %kp=k+B; %t=y(kp)*WNp; % y(kp)=y(k)-t; %p=J*2A(M-L);本次蝶形運算的跨越間隔為 2AL蝶形運算的兩個因子對應單元下標的關系蝶形運算的乘積項蝶形運算， 注意必須先進行減法運算，然后進行加法運算，否則要使用中間變量來傳遞y(k)y(k)=y(k)+t; %蝶形運算endendend%yaxes(handles.tu2);x1,w1=freqz(x1,1); % 繪制原始語音信號的頻率圖plot(w1,abs(x1);title('x1 的頻率響應特性')xlabel(' 頻率');ylabel(' 幅度');axes(handles.tu3);plot(abs(y1) %Matlab自帶的FFT®數實現的語音信號的FF砌譜圖title('Matlab自帶的FFT®數實現的x1的頻譜')xlabel('K');ylabel('Y1(k)');axes(handles.tu4);plot(abs(y); % 編寫的FFTS序實現的語音信號的FF砌譜