1、 DSP原理與應用教程課程設計 DSP原理與應用(yngyng)教程課程設計說明書課題(kt)：基于TMS320VC5416的簡易(jiny)信號分析儀設計專業： 電子信息工程 班級： 電子信息工程1102班 姓名： 陳瑋 學號： 3110209424 指導老師： 方衛東、李錦彬、朱悅涵 2014年 12 月 14日 - 12月 26日 摘要(zhiyo)系統(xtng)基于快速傅立葉變換(FFT)算法(sun f)，以TMS320VC5416為控制與數據處理核心，通過CCS開發環境以及MATLAB模擬實現對音頻信號頻率成分的分析。系統由電源模塊、音頻輸入輸出模塊、時鐘模塊等組成。通過A/DC

2、采樣并進行FFT變換，分析音頻信號的頻譜結構。關鍵詞： 音頻信號,CCS，FFT, 頻譜目錄(ml) TOC o 1-3 h z u HYPERLINK l _Toc407306637 摘要(zhiyo) PAGEREF _Toc407306637 h 2 HYPERLINK l _Toc407306638 一、設計任務(rn wu)與目標 PAGEREF _Toc407306638 h 4 HYPERLINK l _Toc407306639 1.1設計課題 PAGEREF _Toc407306639 h 4 HYPERLINK l _Toc407306640 1.2設計要求 PAGEREF

3、_Toc407306640 h 4 HYPERLINK l _Toc407306641 二、諧波的概念 PAGEREF _Toc407306641 h 4 HYPERLINK l _Toc407306642 三、諧波的分析方法 PAGEREF _Toc407306642 h 4 HYPERLINK l _Toc407306643 2.1模擬電路 PAGEREF _Toc407306643 h 4 HYPERLINK l _Toc407306644 2.2傅立葉變換 PAGEREF _Toc407306644 h 5 HYPERLINK l _Toc407306645 2.3小波變換 PAGER

4、EF _Toc407306645 h 5 HYPERLINK l _Toc407306646 四、DFT/FFT分析原理 PAGEREF _Toc407306646 h 5 HYPERLINK l _Toc407306647 4.1 DFT計算公式 PAGEREF _Toc407306647 h 5 HYPERLINK l _Toc407306648 4.2 N點DFT的計算量 PAGEREF _Toc407306648 h 6 HYPERLINK l _Toc407306649 4.3旋轉因子WN的特性 PAGEREF _Toc407306649 h 6 HYPERLINK l _Toc40

5、7306650 4.4基-2 FFT算法推導 PAGEREF _Toc407306650 h 6 HYPERLINK l _Toc407306651 4.5 N點基-2 FFT算法的計算量 PAGEREF _Toc407306651 h 10 HYPERLINK l _Toc407306652 4.6 N點基-2 FFT算法的實現方法 PAGEREF _Toc407306652 h 10 HYPERLINK l _Toc407306653 4.6.1對于輸入數據序列進行倒位序變換 PAGEREF _Toc407306653 h 10 HYPERLINK l _Toc407306654 4.6.

6、2蝶形運算的循環結構 PAGEREF _Toc407306654 h 11 HYPERLINK l _Toc407306655 4.6.3浮點到定點轉換需要注意的關鍵問題 PAGEREF _Toc407306655 h 11 HYPERLINK l _Toc407306656 4.6.4計算過程中的溢出問題 PAGEREF _Toc407306656 h 12 HYPERLINK l _Toc407306657 五、硬件系統 PAGEREF _Toc407306657 h 12 HYPERLINK l _Toc407306658 5.1硬件電路的總體框圖 PAGEREF _Toc4073066

7、58 h 12 HYPERLINK l _Toc407306659 5.2各子模塊的組成 PAGEREF _Toc407306659 h 13 HYPERLINK l _Toc407306660 5.2.1電源模塊 PAGEREF _Toc407306660 h 13 HYPERLINK l _Toc407306661 5.2.2時鐘模塊 PAGEREF _Toc407306661 h 14 HYPERLINK l _Toc407306662 5.2.3音頻輸入輸出模塊 PAGEREF _Toc407306662 h 15 HYPERLINK l _Toc407306663 5.2.4 TMS

8、320VC5416和CODEC接口 PAGEREF _Toc407306663 h 16 HYPERLINK l _Toc407306664 5.2.5 TMS320VC5416DSP最小系統 PAGEREF _Toc407306664 h 17 HYPERLINK l _Toc407306665 六、基于MATLAB的程序仿真 PAGEREF _Toc407306665 h 18 HYPERLINK l _Toc407306666 6.1程序流程圖 PAGEREF _Toc407306666 h 18 HYPERLINK l _Toc407306667 6.2 MATLAB構造音頻信號 PA

9、GEREF _Toc407306667 h 19 HYPERLINK l _Toc407306668 6.3仿真步驟 PAGEREF _Toc407306668 h 19 HYPERLINK l _Toc407306669 6.4仿真結果 PAGEREF _Toc407306669 h 20 HYPERLINK l _Toc407306670 6.5結果分析 PAGEREF _Toc407306670 h 20 HYPERLINK l _Toc407306671 七、課程設計總結 PAGEREF _Toc407306671 h 21 HYPERLINK l _Toc407306672 八、參考

10、文獻 PAGEREF _Toc407306672 h 22 HYPERLINK l _Toc407306673 附錄 PAGEREF _Toc407306673 h 22一、設計(shj)任務與目標1.1設計(shj)課題基于TMS320VC5416的簡易(jiny)信號分析儀設計1.2設計要求1.給出算法原理2.寫出應用軟件流程圖3.以TMS320VC5416 DSP為核心，設計一DSP應用系統，用DSP C語言和匯編混合編程的方法設計應用軟件，實現音頻信號功率譜分析，并驗證最終的分析結果。最后按要求寫出設計說明書，繪出設計的軟硬件系統.1.3設計目標 學習FFT算法在功率譜分析中的應用；以

11、TMS320VC5416 DSP為核心CPU，設計一小型DSP應用系統，編制軟件實現頻譜分析，在CCS上觀察分析的結果.二、諧波的概念任何周期性波形均可分解為一個基頻正弦波加上許多高次頻率的正弦波，高次頻率是基頻的整倍數（N，只能為整數），直流成分稱為0次諧波，基波稱為1次諧波，二次以上的波形稱為高次諧波，其中偶次頻率的波形稱為偶次諧波，奇次頻率的波形稱為奇次諧波。三、諧波的分析方法2.1.模擬電路 消除諧波的方法很多，即有主動型，又有被動型；既有無源的，也有有源的，還有混合型的，目前較為先進的是采用有源電力濾波器。但由于其檢測環節多采用模擬電路，因而造價較高，且由于模擬帶通濾波器對頻率和溫度

12、的變化非常敏感，故使其基波幅值誤差很難控制在10%以內，嚴重影響了有源濾波器的控制性能。近年來，人工神經 HYPERLINK /network/ 網絡的研究取得了較大進展，由于神經元有自適應和自學習能力，且結構簡單，輸入輸出關系明了，因此可用神經元替代自適應濾波器，再用一對與基波頻率相同，相位相差90度的正弦向量作為神經元的輸入。由神經元先得到基波電流，然后檢測出應補償的電流，從而完成諧波電流的檢測。但人工神經網絡的硬件目前還是一個比較薄弱的環節，限制了其應用范圍。2.2.傅立葉變換(binhun)利用傅立葉變換可在數字域進行(jnxng)諧波檢測，電力系統的諧波分析，目前大都是通過該方法實現

13、的，離散傅立葉變換所需要處理的是經過采樣和A/D轉換得到的數字信號，設待測信號為x(t)，采樣間隔為 t秒，采樣頻率 =1/ t滿足采樣定理，即 大于信號最高頻率分量的2倍，則采樣信號為x(n t)，并且采樣信號總是有限長度的，即n=0，1N-1。這相當于對無限長的信號做了截斷，因而造成了傅立葉變換的泄露現象，產生誤差。此外，對于離散傅立葉變換來說，如果不是整數周期采樣，那么即使信號只含有單一頻率，離散傅立葉變換也不可能求出信號的準確參數，因而出現柵欄效應。通過加窗可以減小泄露現象的影響。2.3.小波變換(binhun)小波變換已廣泛應用于信號分析、語音識別與合成、自動控制、圖象處理與分析等領

14、域。電力諧波是由各種頻率成分合成的、隨機的、出現和消失都非常突然的信號，在應用離散傅立葉變換進行處理受到局限的情況下，可充分發揮小波變換的優勢。即對諧波采樣離散后，利用小波變換對數字信號進行處理，從而實現對諧波的精確測定。小波可以看作是一個雙窗函數，對一信號進行小波變換相當于從這一時頻窗內的信息提取信號。對于檢測高頻信息，時窗變窄，可對信號的高頻分量做細致的觀測；對于分析低頻信息，這時時窗自動變寬，可對信號的低頻分量做概貌分析。所以小波變換具有自動“調焦”性。其次，小波變換是按頻帶而不是按頻點的方式處理頻域信息，因此信號頻率的微小波動不會對處理產生很大的影響，并不要求對信號進行整周期采樣。另外

15、，由小波變換的時間局部可知，在信號的局部發生波動時，不會象傅立葉變換那樣把影響擴散到整個頻譜，而只改變當時一小段時間的頻譜分布，因此，采用小波變換可以跟蹤時變和暫態信號。四、DFT/FFT分析原理4.1 DFT計算公式其中x(n)表示輸入的離散數字信號序列，WN為旋轉因子，X(k)為輸入序列x(n)對應(duyng)的N個離散頻率點的相對幅度。一般情況下，假設x(n)來自于低通采樣，采樣頻率為fs，那么X(k)表示了從-fs/2率開始，頻率間隔為fs/N,到fs/2-fs/N截至的N個頻率點的相對幅度。因為DFT計算得到的一組離散頻率幅度值實際上是在頻率軸上從成周期變化的，即X(k+N)=X(

16、k)。因此任意取連續的N個點均可以表示DFT的計算效果，負頻率成分比較抽象，難于理解，根據X(k)的周期特性，于是我們又可以認為X(k)表示了從零頻率開始，頻率間隔為fs/N,到fs-fs/N截至的N個頻率點的相對幅度。4.2 N點DFT的計算(j sun)量根據(1)式給出的DFT計算公式，我們可以知道每計算一個頻率點X(k)均需要進行N次復數乘法(chngf)和N-1次復數加法，計算N各點的X(k)共需要N2次復數乘法和N*(N-1)次復數加法。當x(n)為實數的情況下，計算N點的DFT需要2*N2次實數乘法，2*N*(N-1)次實數加法。4.3旋轉因子WN的特性1.WN的對稱性 2.WN

17、的周期性3.WN的可約性根據以上這些性質，我們可以得到式(5)的一系列有用結果4.4基-2 FFT算法(sun f)推導假設采樣序列點數為N=2L，L為整數，如果不滿足這個條件可以人為地添加若干個0以使采樣序列點數滿足這一要求。首先我們將序列x(n)按照(nzho)奇偶分為兩組如下：于是(ysh)根據DFT計算公式(1)有：至此，我們將一個N點的DFT轉化為了式(7)的形式，此時k的取值為0到N-1,現在分為兩段來討論，當k為0N/2-1的時候，因為x1(r)，x2(r)為N/2點的序列，因此，此時式(7)可以寫為：而當 k取值為N/2N-1時，k用k+N/2取代，k取值為0N/2-1。對式(

18、7)化簡可得：綜合以上推導我們可以(ky)得到如下結論：一個N點的DFT變換過程可以用兩個N/2點的DFT變換過程來表示，其具體公式如式(10)所示DFT快速算法的迭代公式：上式中X(k)為偶數項分支(fnzh)的離散傅立葉變換，X(k)為奇數項分支的離散傅立葉變換。 式(10)的計算過程可以用圖1的蝶形算法流圖直觀地表示出來。圖1 時間(shjin)抽取法蝶形運算流圖在圖1中，輸入為兩個N/2點的DFT輸出為一個N點的DFT結果，輸入輸出點數一致。運用這種表示方法，8點的DFT可以用圖2來表示：圖2 8點DFT的4點分解(fnji)根據公式(10)，一個N點的DFT可以由兩個N/2點的DFT

19、運算構成，再結合(jih)圖1的蝶形信號流圖可以得到圖2的8點DFT的第一次分解。該分解可以用以下幾個步驟來描述： 1.將N點的輸入序列(xli)按奇偶分為2組分別為N/2點的序列 2.分別對1中的每組序列進行DFT變換得到兩組點數為N/2的DFT變換值X1和X2 3.按照蝶形信號流圖將2的結果組合為一個N點的DFT變換結果根據式(10)我們可以對圖2中的4點DFT進一步分解，得到圖3的結果，分解步驟和前面一致。圖3 8點DFT的2點分解最后對2點DFT進一步分解得到最終的8點FFT信號計算流圖：圖4 8點DFT的全分解(fnji)從圖2到圖4的過程中關于旋轉系數的變化規律需要說明一下。看起來

20、似乎向前推一級，在奇數分組部分的旋轉系數因子增量(zn lin)似乎就要變大，其實不是這樣。事實上奇數分組部分的旋轉因子指數每次增量固定為1，只是因為每向前推進一次，該分組序列的數據個數變少了，為了統一使用以原數據N為基的旋轉因子就進行了變換導致的。每一次分組奇數部分的系數WN,這里的N均為本次分組前的序列點數。以上邊的8點DFT為例，第一次分組N=8,第二次分組N為4，為了統一根據式(4)進行了變換將N變為了8，但指數相應的需要乘以2。4.5 N點基-2 FFT算法(sun f)的計算量從圖4可以看到N點DFT的FFT變換可以轉為log2(N)級級聯的蝶形運算，每一級均包含有N/2次蝶形計算

21、。而每一個蝶形運算包含了1次復數乘法，2次復數加法。因此N點FFT計算的總計算量為：復數乘法N/2log2(N) 復數加法Nlog2(N)。假設被采樣的序列為實數序列，那么也只有第一級的計算為實數與復數的混合計算，經過一次迭代后來的計算均變為復數計算，在這一點上和直接的DFT計算不一致。因此對于輸入序列是復數還是實數對FFT算法的效率影響較小。一次復數乘法包含了4次實數乘法，2次實數加法，一次復數加法包含了2次復數加法。因此對于N點的FFT計算需要總共的實數乘法數量為：2Nlog2(N)；總的復數加法次數為：2xNxlog2(N)。4.6 N點基-2 FFT算法的實現方法從圖4我們可以總結出對

22、于點數為N=2L的DFT快速計算方法的流程：4.6.1對于輸入數據序列進行倒位序變換該變換的目的是使輸出能夠得到X(0)X(N-1)的順序序列，同樣以8點DFT為例，該變換將順序輸入序列x(0)x(7)變為如圖4的x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其實現方法是：假設順序輸入序列一次村在A(0)A(N-1)的數組元素中，首先我們將數組下標進行二進制化(例：對于點數為8的序列只需要LOG2(8) = 3位二進制序列表示，序號6就表示為110)。二進制化以后就是將二進制序列進行倒位，倒位的過程就是將原序列從右到左書寫一次構成新的序列，例如序號為6的二

23、進制表示為110，倒位后變為了011，即使十進制的3。第三步就是將倒位前和倒位后的序號對應的數據互換。依然以序號6為例，其互換過程如下：Temp = A(6); A(6) = A(3); A(3) = A(6);實際上考慮(kol)到執行效率，如果對于每一次輸入的數據都需要這個處理過程是非常浪費時間的。我們可以采用指向指針的指針來實現該過程，或者是采用指針數組來實現該過程。4.6.2蝶形運算(yn sun)的循環結構從圖4中我們可以看到對于點數為N = 2L的fft運算，可以分解為L階蝶形圖級聯，每一階蝶形圖內又分為M個蝶形組，每個蝶形組內包含K個蝶形。根據這一點我們就可以構造三階循環來實現蝶

24、形運算。編程過程需要注意旋轉因子(ynz)與蝶形階數和蝶形分組內的蝶形個數存在關聯。4.6.3浮點到定點轉換需要注意的關鍵問題上邊的分析都是基于浮點運算來得到的結論，事實上大多數嵌入式系統對浮點運算支持甚微，因此在嵌入式系統中進行離散傅里葉變換一般都應該采用定點方式。對于簡單的DFT運算從浮點到定點顯得非常容易。根據式(1)，假設輸入x(n)是經過AD采樣的數字序列，AD位數為12位，則輸入信號范圍為04096。為了進行定點運算我們將旋轉因子實部虛部同時擴大212倍，取整數部分代表旋轉因子。之后，我們可以按照(1)式計算，得到的結果與原結果成比例關系，新的結果比原結果的212倍。但是，對于使用

25、蝶形運算的fft我們不能采用這種簡單的放大旋轉因子轉為整數計算的方式。因為fft是一個非對稱迭代過程，假設我們對旋轉因子進行了放大，根據蝶形流圖我們可以發現其最終的結果是，不同的輸入被放大了不同的倍數，對于第一個輸入x(0)永遠也不會放大。舉一個更加形象的例子，還是以圖4為例。從圖中可以看出右側的X(0)可以直接用下式表示： 可以通過簡單的右移位運算將之前的放大(fngd)抵消，而右移位又代替了 從上式我們(w men)可以看到不同輸入項所乘的旋轉因子個數(注意這里是個數，就算是wn0，也被考慮進去了，因為在沒有放大(fngd)時wn0等于1，放大后所有旋轉因子指數模均不為1，因此需要考慮)。

26、這就導致輸入不平衡，運算結果不正確。經查閱相關資料，比較妥善的做法是，首先對所有旋轉因子都放大2Q倍，Q必須要大于等于L，以保證不同旋轉因子的差異化。旋轉因子放大，為了保證其模為1，在每一次蝶形運算的乘積運算中我們需要將結果右移Q位來抵消這個放大，從而得到正確的結果。之所以采用放大倍數必須是2的整數次冪的原因也在于此，我們之后可以通過簡單的右移位運算將之前的放大抵消，而右移位又代替了除法運算，大大節省了時間。4.6.4計算過程中的溢出問題最后需要注意的一個問題就是計算過程中的溢出問題。在實際應用中，AD雖然有12位的位寬，但是采樣得到的信號可能較小，例如可能在08之間波動，也就是說實際可能只有

27、3位的情況。這種情況下為了在計算過程中不丟失信息，一般都需要先將輸入數據左移P位進行放大處理，數據放大可能會導致溢出，從而使計算錯誤，而溢出的極限情況是這樣：假設我們數據位寬為D位(不包括符號位)，AD采樣位數B位，數字放大倍數P位，旋轉因此放大倍數Q位，FFT級聯運算帶來的最大累加倍數L位。我們得到:假設AD位寬12，數據位寬32，符號位1位，因此有效位寬31位，采樣點數N，那么我們可以得到log2(N)+P+Q=L可以得到放大倍數P VECT PAGE 0.sysregs: BIOSREGS PAGE 1 .trcinit: EPROG PAGE 0 .gblinit: EPROG PAG

28、E 0 frt: EPROG PAGE 0 .text: EPROG PAGE 0/* 已初始化段，包含所有的可執行代碼、字符串和編譯器生成的常數 */ .cinit: EPROG PAGE 0/* 已初始化段，包含初始化變量和常數表 */ .pinit: EPROG PAGE 0/* 已初始化段，包含運行時調用的全局(qunj)目標構造器列表 */ .sysinit: EPROG PAGE 0 .bss: IDATA PAGE 1/* 未初始化段 */ .far: IDATA PAGE 1 .const: IDATA PAGE 1/* 已初始化段，包含(bohn)由C限定詞const定義的字

29、符串常量和數據 */ .switch: IDATA PAGE 1/* 已初始化段 */ .sysmem: IDATA PAGE 1/* 未初始化段，為動態(dngti)存儲器函數malloc等分配存儲器空間 */ .cio: IDATA PAGE 1 .MEM$obj: IDATA PAGE 1 .sysheap: IDATA PAGE 1 .stack: IDATA PAGE 1/* 未初始化段，為系統堆棧分配存儲器 */ fft.c#include #include #include #include #include fft.h#include sin.incconst float da

30、taR256=/y=100*sin(2*pi/N)*2*n)+200*sin(2*pi/N)*4*n)+300*sin(2*pi/N)*6*n)0,68.529,135.91,201,262.72,320.05,372.05,417.89,456.85,488.37,511.99,527.43,534.57,533.42,524.17,507.16,482.84,451.84,414.88,372.78,326.45,276.88,225.08,172.11,119,66.802,16.49,-31.003,-74.826,-114.22,-148.55,-177.27,-200,-216.4

31、8,-226.58,-230.34,-227.9,-219.56,-205.74,-186.96,-163.84,-137.1,-107.51,-75.891,-43.102,-10,22.563,53.768,82.843,109.08,131.86,150.65,165.02,174.66,179.4,179.16,174.01,164.13,149.82,131.49,109.64,84.883,57.869,29.322,7.3479e-014,-29.322,-57.869,-84.883,-109.64,-131.49,-149.82,-164.13,-174.01,-179.16

32、,-179.4,-174.66,-165.02,-150.65,-131.86,-109.08,-82.843,-53.768,-22.563,10,43.102,75.891,107.51,137.1,163.84,186.96,205.74,219.56,227.9,230.34,226.58,216.48,200,177.27,148.55,114.22,74.826,31.003,-16.49,-66.802,-119,-172.11,-225.08,-276.88,-326.45,-372.78,-414.88,-451.84,-482.84,-507.16,-524.17,-533

33、.42,-534.57,-527.43,-511.99,-488.37,-456.85,-417.89,-372.05,-320.05,-262.72,-201,-135.91,-68.529,1.6111e-012,68.529,135.91,201,262.72,320.05,372.05,417.89,456.85,488.37,511.99,527.43,534.57,533.42,524.17,507.16,482.84,451.84,414.88,372.78,326.45,276.88,225.08,172.11,119,66.802,16.49,-31.003,-74.826,

34、-114.22,-148.55,-177.27,-200,-216.48,-226.58,-230.34,-227.9,-219.56,-205.74,-186.96,-163.84,-137.1,-107.51,-75.891,-43.102,-10,22.563,53.768,82.843,109.08,131.86,150.65,165.02,174.66,179.4,179.16,174.01,164.13,149.82,131.49,109.64,84.883,57.869,29.322,-3.1247e-013,-29.322,-57.869,-84.883,-109.64,-13

35、1.49,-149.82,-164.13,-174.01,-179.16,-179.4,-174.66,-165.02,-150.65,-131.86,-109.08,-82.843,-53.768,-22.563,10,43.102,75.891,107.51,137.1,163.84,186.96,205.74,219.56,227.9,230.34,226.58,216.48,200,177.27,148.55,114.22,74.826,31.003,-16.49,-66.802,-119,-172.11,-225.08,-276.88,-326.45,-372.78,-414.88,

36、-451.84,-482.84,-507.16,-524.17,-533.42,-534.57,-527.43,-511.99,-488.37,-456.85,-417.89,-372.05,-320.05,-262.72,-201,-135.91,-68.529;float re_inputK_FFT_SIZE,im_inputK_FFT_SIZE;/ FFT/DFT運算(yn sun)輸入float am_outputK_FFT_SIZE;/ FFT/DFT輸出(shch)諧波幅值float re_outputK_FFT_SIZE,im_outputK_FFT_SIZE;/ DFT運算(y

37、n sun)輸出諧波的實部和虛部void clr_ram(void) unsigned int i;for (i=0;iK_FFT_SIZE;i+) am_outputi=0.0;re_inputi=im_inputi=0.0;re_outputi=im_outputi=0.0;/ 用于DFTvoid put_in(unsigned int N) int i;for(i=0;iN;i+) re_inputi=dataR256i/1.0;/ 實部im_inputi=0.0;/ 虛部/ 位倒序(do x)運算void re_bit(unsigned int N)int i,j,M,s;float

38、temp;M=N/2;i=0;for(j=1;j=s) i=i-s;s=s/2;i=i+s;if(ji) temp=re_inputj;re_inputj=re_inputi;re_inputi=temp;temp=im_inputj;im_inputj=im_inputi;im_inputi=temp;/ 蝶形運算(yn sun)void butterfly(unsigned int l) int i,j,r;/分別(fnbi)是級號、蝶群號、運算蝶號int m1,m2,m3,m4;int N,k1,k2;float u,v;int x;N=1l;for(i=1;il+1;i+)/第i級/m

39、1=pow(2,i-1);/m1=1,2,4,8,16,32,64,128,.m1=1(i-1);m2=2*m1;/第i級的蝶群長度m3=N/m2;/第i級包含的蝶形單元，即“群”個數for(j=0;jm3;j+)/第j號蝶群m4=j*m2;/m4=j*(2i) 蝶群址for(r=0;rm1;r+)/第r號運算蝶k1 = (m4+r);/s=j*(2i)+r 蝶址k2 = (k1+m1);/u=Rek2*cos(2*p*r/m2)+Imk2*sin(2*p*r/m2);/v=Imk2*cos(2*p*r/m2)-Rek2*sin(2*p*r/m2);x = (1(l-i);x = x*r;u

40、= (re_inputk2*cos_tab256x)+(im_inputk2*sin_tab256x);v= (im_inputk2*cos_tab256x)-(re_inputk2*sin_tab256x);re_inputk2 = (re_inputk1-u);im_inputk2 = (im_inputk1-v);re_inputk1 = (re_inputk1+u);im_inputk1 = (im_inputk1+v);/ FFT運算(yn sun)/ input: n=采樣(ci yn)點數N（必須是2的l次冪）void fft(unsigned int n,unsigned in

41、t l) re_bit(n);/ 位倒序(do x)butterfly(l);/ DFT運算，計算結果為X(0),X(1),.X(N-1)/ sign=1:正變換;sign=-1:反變換void dft(int sign) int i,k;double c,q,w,s;q=(2.0*PI)/K_FFT_SIZE;for (k=0; kK_FFT_SIZE; k+) w=k*q;re_outputk=im_outputk=0.0;for (i=0; iK_FFT_SIZE; i+) unsigned int l;/d=i*w;/c=cos(d);/s=sin(d)*sign;l=(k*i)%K_

42、FFT_SIZE;c=cos_tab256l;s=sin_tab256l;re_outputk += (c*re_inputi+s*im_inputi);/ 計算(j sun)X(k)的實部im_outputk += (c*im_inputi-s*re_inputi);/ 計算(j sun)X(k)的虛部if (sign=-1) c=1.0/K_FFT_SIZE;for (k=0; kK_FFT_SIZE; k+) re_outputk=c*re_outputk;im_outputk=c*im_outputk;/ 計算(j sun)單次X(k)，即只計算某一次諧波的值/ input:k=諧波次

43、數void dft_k(unsigned int k) int i;double c,s;re_outputk=im_outputk=0.0;for (i=0; iK_FFT_SIZE; i+) unsigned int l;l=(k*i)%K_FFT_SIZE;c=cos_tab256l;s=sin_tab256l;re_outputk += (c*re_inputi+s*im_inputi);/ 計算X(k)的實部 im_outputk += (c*im_inputi-s*re_inputi);/ 計算X(k)的虛部/ 求FFT運算結果諧波的幅值/ input:n=采樣點數=諧波個數/ o

44、utput:am_output=諧波的幅值void get_am(unsigned int N) unsigned int i;for (i=0; iN; i+)am_outputi=sqrt(re_inputi*re_inputi+im_inputi*im_inputi);am_outputi /= (float)K_FFT_SIZE;am_outputi *= 2;if (am_outputi ZERO_VALVE)/ 屏蔽計算帶來的0的誤差am_outputi = 0.0;void get_am_dft(void) unsigned int i;for (i=0; iK_FFT_SIZE

45、; i+)am_outputi=sqrt(re_outputi*re_outputi+im_outputi*im_outputi);am_outputi /= (float)K_FFT_SIZE;am_outputi *= 2;if (am_outputi ZERO_VALVE)am_outputi = 0.0;#defineEXE_FFT0/ FFT運算(yn sun)#defineEXE_DFT1/ DFT計算(j sun)所有諧波#defineEXE_DFT_K2/ DFT計算指定(zhdng)次諧波#defineEXE_MODEEXE_FFTvoid main() int m;clr_

46、ram();/ 清內存put_in(K_FFT_SIZE);/ 偽A/D采樣#if EXE_MODE=EXE_FFT/ 快速傅里葉變換FFTfft(K_FFT_SIZE, K_LOGN);/ 得X(k),k=0,1,.,N-1get_am(K_FFT_SIZE);/ 求諧波幅值/ 快速傅里葉逆變換IFFTfor (m=0; mK_FFT_SIZE; m+)/ 求X(k)共軛復數im_inputm=-im_inputm;fft(K_FFT_SIZE, K_LOGN);for (m=0; mK_FFT_SIZE; m+)/ 最后取共軛并乘以1/N得到序列x(n),n=0,1,.,N-1im_inp

47、utm=-im_inputm/K_FFT_SIZE;re_inputm/=K_FFT_SIZE;if (abs(re_inputm) ZERO_VALVE)re_inputm = 0.0;if (abs(im_inputm) ZERO_VALVE)im_inputm = 0.0;asm( NOP);#endif#if EXE_MODE=EXE_DFTdft(1);get_am_dft();asm( NOP);#endif#if EXE_MODE=EXE_DFT_Kdft_k(5);get_am_dft();asm( NOP);#endifasm( NOP);while(1);基于MATLAB音

48、頻(ynpn)構造代碼function st=mstgN=256Fs=2800;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/256*2;fc2=Fs/256*4;fc3=Fs/256*6;xt1=100*sin(2*pi*fc1*t);xt2=200*sin(2*pi*fc2*t);xt3=300*sin(2*pi*fc3*t);st=xt1+xt2+xt3;fxt=fft(st,N);dlmwrite(dsp.txt,st);subplot(3,1,1)plot(t,st);grid;xlabel(t/s);ylabel(s(t);ax

49、is(0,Tp,min(st),max(st);title(a)st的波形(b xn);subplot(3,1,2)stem(f,abs(fxt)/max(abs(fxt),.);grid;title(b)st的頻譜);axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel(幅度(fd);基于MATLAB的SIN/COS的查表法代碼fidc=fopen(sin.inc,wt);fprintf(fidc,%s,const float cos_tab256=);fprintf(fidc,%s,1.0000,);for n=1:255a(n)=cos(2*pi/255)*1*n

50、);fprintf(fidc,%.4f,a(n);endfprintf(fidc,%s,);fprintf(fidc,nnn);fprintf(fidc,%s,const float sin_tab256=);fprintf(fidc,%s,0,);for n=1:255a(n)=sin(2*pi/255)*1*n);fprintf(fidc,%.4f,a(n);endfprintf(fidc,%s,);fclose(fidc);sin.inc: const floatcos_tab256=1.0000,0.9997,0.9988,0.9973,0.9951,0.9924,0.9891,0.9

51、852,0.9806,0.9755,0.9698,0.9635,0.9566,0.9491,0.9411,0.9325,0.9233,0.9135,0.9032,0.8924,0.8810,0.8691,0.8566,0.8437,0.8302,0.8162,0.8017,0.7867,0.7713,0.7554,0.7390,0.7222,0.7049,0.6872,0.6691,0.6506,0.6317,0.6124,0.5928,0.5727,0.5524,0.5317,0.5106,0.4893,0.4677,0.4457,0.4235,0.4011,0.3784,0.3555,0.

52、3324,0.3090,0.2855,0.2618,0.2379,0.2139,0.1898,0.1656,0.1412,0.1168,0.0923,0.0677,0.0431,0.0185,-0.0062,-0.0308,-0.0554,-0.0800,-0.1045,-0.1290,-0.1534,-0.1777,-0.2019,-0.2260,-0.2499,-0.2737,-0.2973,-0.3207,-0.3439,-0.3670,-0.3898,-0.4124,-0.4347,-0.4567,-0.4785,-0.5000,-0.5212,-0.5421,-0.5626,-0.5

53、828,-0.6026,-0.6221,-0.6412,-0.6599,-0.6782,-0.6961,-0.7136,-0.7307,-0.7473,-0.7634,-0.7791,-0.7943,-0.8090,-0.8233,-0.8370,-0.8502,-0.8629,-0.8751,-0.8868,-0.8979,-0.9085,-0.9185,-0.9280,-0.9369,-0.9452,-0.9529,-0.9601,-0.9667,-0.9727,-0.9781,-0.9830,-0.9872,-0.9908,-0.9939,-0.9963,-0.9981,-0.9993,

54、-0.9999,-0.9999,-0.9993,-0.9981,-0.9963,-0.9939,-0.9908,-0.9872,-0.9830,-0.9781,-0.9727,-0.9667,-0.9601,-0.9529,-0.9452,-0.9369,-0.9280,-0.9185,-0.9085,-0.8979,-0.8868,-0.8751,-0.8629,-0.8502,-0.8370,-0.8233,-0.8090,-0.7943,-0.7791,-0.7634,-0.7473,-0.7307,-0.7136,-0.6961,-0.6782,-0.6599,-0.6412,-0.6

55、221,-0.6026,-0.5828,-0.5626,-0.5421,-0.5212,-0.5000,-0.4785,-0.4567,-0.4347,-0.4124,-0.3898,-0.3670,-0.3439,-0.3207,-0.2973,-0.2737,-0.2499,-0.2260,-0.2019,-0.1777,-0.1534,-0.1290,-0.1045,-0.0800,-0.0554,-0.0308,-0.0062,0.0185,0.0431,0.0677,0.0923,0.1168,0.1412,0.1656,0.1898,0.2139,0.2379,0.2618,0.2

56、855,0.3090,0.3324,0.3555,0.3784,0.4011,0.4235,0.4457,0.4677,0.4893,0.5106,0.5317,0.5524,0.5727,0.5928,0.6124,0.6317,0.6506,0.6691,0.6872,0.7049,0.7222,0.7390,0.7554,0.7713,0.7867,0.8017,0.8162,0.8302,0.8437,0.8566,0.8691,0.8810,0.8924,0.9032,0.9135,0.9233,0.9325,0.9411,0.9491,0.9566,0.9635,0.9698,0.9755,0.9806,0.9852,0.9891,0.9924,0.9951,0.9973,0.9988,0.9997,1.0000,const float sin_tab256=0,0.0246,0.0493,0.0739,0.0984,0.1229,0.1473,0.1716,0.1958,0.2199,0.243