基于DSP的實數(shù)序列FFT算法探討

郭建濤，范春鳳

(信陽師范學(xué)院物理電子工程學(xué)院，河南信陽 464000）

“數(shù)字信號處理”是電氣信息類學(xué)生的一門專業(yè)基礎(chǔ)課程，該課程以算法為核心，其特點是理論性強(qiáng)，抽象概念多，對數(shù)學(xué)基礎(chǔ)要求高，傳統(tǒng)的以教師“單純授課”、學(xué)生死記公式的教育模式不僅教學(xué)效果差，而且使得學(xué)生學(xué)習(xí)興趣下降，以致今后不愿從事該領(lǐng)域的工作。而“DSP技術(shù)及應(yīng)用”課程重點偏重數(shù)字信號處理的實現(xiàn)，學(xué)生盡管有濃厚的興趣，但是由于前期的理論基礎(chǔ)有所欠缺，學(xué)習(xí)“DSP技術(shù)及應(yīng)用”課程時比較吃力，在理論深度、算法設(shè)計思想上不能融會貫通，難以達(dá)到培養(yǎng)創(chuàng)新能力、實踐能力的根本要求。

傅里葉變換是一種將信號從時域轉(zhuǎn)變?yōu)轭l域表示的變換手段，它在信號的頻譜成分分析中起著基礎(chǔ)性的作用。由于快速傅里葉變換（FFT）算法的出現(xiàn)，使得傅里葉變換在信號分析和系統(tǒng)設(shè)計中得到了廣泛的應(yīng)用。DSP技術(shù)由于其快速處理大數(shù)據(jù)量的能力在現(xiàn)代語音、圖像以及視頻應(yīng)用中得到了快速發(fā)展，基于DSP芯片實現(xiàn)FFT算法顯得日益重要。

本文從基本FFT算法理論出發(fā)，基于實際應(yīng)用中信號多為實數(shù)數(shù)據(jù)，給出了實數(shù)序列的FFT快速算法，并結(jié)合當(dāng)前通用的TITMS320C54系列DSP芯片，匯編編程實現(xiàn)FFT算法。經(jīng)過理論分析和實驗仿真，更加有助于加強(qiáng)對FFT算法內(nèi)容的掌握和理解，大大提高了數(shù)字信號處理技術(shù)的教學(xué)效果。

1 FFT算法基本原理

1.1 N點復(fù)數(shù)序列的FFT

給定N點復(fù)數(shù)序列d（n），0≤n≤N-1。為計算方便起見，這里設(shè)N=2M，M為正整數(shù)。

d（n）的DFT為[1-2]：

將d（n）按照n的奇偶分為兩組，即按n=2r及n=2r+1分為兩組：

1.2 2N點實數(shù)序列的FFT

對于2N 點實數(shù)序列 a（n），n=0，1，…，2N-1，按照如下方式構(gòu)成一N 點復(fù)數(shù)序列 d（n），n=0,1,…，N-1，a(2n)表示其實部，a(2n+1)是其虛部：

根據(jù)DFT的共軛對稱性，知道d（n）的實部的DFT對應(yīng)于D（k）的共軛偶對稱分量，d（n）的虛部的DFT對應(yīng)于D（k）的共軛奇對稱分量：

從而得到a（n）的DFT為：

將式（4）代入式（5），得到

從而，一個2N點的實數(shù)序列的傅里葉變換可以通過一個N點復(fù)數(shù)序列的FFT獲取。

2 FFT的DSP程序?qū)崿F(xiàn)

下面以2N=256點的實數(shù)序列為例說明FFT的DSP實現(xiàn)過程。

2.1 實數(shù)FFT運(yùn)算序列的存儲分配

程序代碼安排在0x3000h開始的存儲器中，其中0x3000h～0x3080h存放中斷向量表，F(xiàn)FT程序使用的正弦表和余弦表數(shù)據(jù)(.data段)安排在0xc00h開始的地方，長度為0x1400h（分別有256個值），變量(.bss段定義)存放在0x80h開始的地址中。設(shè)置256點實數(shù)FFT程序的數(shù)據(jù)緩沖區(qū)位于0x2300h～0x23ffh,FFT變換后功率譜的計算結(jié)果存放在0x2200h～0x22ffh中[3-4]。

2.2 基2實數(shù)FFT運(yùn)算的算法

該算法主要分為以下2個步驟進(jìn)行。

第一，輸入數(shù)據(jù)的組合和位排序。首先，原始輸入的2N=256個點的實數(shù)序列復(fù)制放到0x2300h開始的相鄰單元，當(dāng)成N=128點的復(fù)數(shù)序列d[n]，其中奇數(shù)地址是d[n]實部，偶數(shù)地址是d[n]的虛部。然后，依據(jù)圖1輸入數(shù)據(jù)的排列，把復(fù)數(shù)序列位倒序后，存儲在以0x2200h開始的數(shù)據(jù)處理緩沖區(qū)中。

程序設(shè)計中，使用位倒序?qū)ぶ贩绞娇梢源蟠筇岣叱绦驁?zhí)行的速度和使用存儲器的效率。在這種尋址方式中，AR0存放的整數(shù)是實數(shù)數(shù)據(jù)點數(shù)的一半（此處僅僅完成位倒序），即128。輸入緩沖指針AR3進(jìn)行加1減1、位倒序緩沖指針AR2連續(xù)加2次，完成一個復(fù)數(shù)數(shù)據(jù)倒序。然后按照位倒序?qū)ぶ贩绞叫薷腁R3，即，指令：MAR AR3+0B，就可以得到下一數(shù)據(jù)地址。重復(fù)上述過程128次，即可完成輸入數(shù)據(jù)位的排序工作；這一點通過設(shè)置塊重復(fù)計數(shù)器BRC由RPTB塊重復(fù)指令完成[5]。

第二，N點復(fù)數(shù)FFT運(yùn)算。在數(shù)據(jù)處理器里進(jìn)行N點復(fù)數(shù)FFT運(yùn)算。由于在FFT運(yùn)算中要用到旋轉(zhuǎn)因子WN，它是一個復(fù)數(shù)。我們把它分為正弦和余弦部分，用Q15格式將它們存儲在兩個分離的表中。FFT變換關(guān)鍵運(yùn)算主要包括蝶形結(jié)個數(shù)和單個蝶形結(jié)計算兩個部分。

(1)級、組和蝶形結(jié)的個數(shù)變量及其調(diào)整

N點復(fù)數(shù)序列FFT運(yùn)算的分成M級，每一級分別進(jìn)行2，22，…，2M點的FFT運(yùn)算，此時每一級的組數(shù)分別為2M-1，2M-2，…，1，每一級各組中的蝶形結(jié)的數(shù)目相同，不同級的蝶形結(jié)數(shù)目分布為20，21，…，2M-1（FFT運(yùn)算點數(shù)的一半）。因此，算法首先制定級、組和蝶形結(jié)三個計數(shù)器，通過三層循環(huán)分別完成各級、各組合單個蝶形結(jié)的運(yùn)算。

計算時，由2個數(shù)據(jù)指針PX、QX分別指向參加蝶形運(yùn)算的第一個數(shù)據(jù)和第二個數(shù)據(jù)。其中，第二個數(shù)據(jù)需要乘以旋轉(zhuǎn)因子（通過指向正余弦表的指針AR4、AR5完成定位）。

(2)蝶形結(jié)運(yùn)算

設(shè)定參加蝶形結(jié)運(yùn)算的兩個輸入數(shù)據(jù)分別為P、Q，包括數(shù)據(jù)的實部和虛部，用（PR、PI）和（QR、QI）表示；輸出分別 P′，Q′，旋轉(zhuǎn)因子用 W 表示，相應(yīng)的實部和虛部表示用其后添加R和I加以明示，見圖1。

圖1 單個蝶形運(yùn)算

對于旋轉(zhuǎn)因子W，正余弦表格給出的僅僅是正余弦值，用Q15格式將它們存儲在兩個分離的表中。每個表中有128項，對應(yīng)0°~180°。因為采用循環(huán)尋址地址來對表尋址，128=27<28，因此每張表排隊的開始地址就必須是8個LSB位為0的地址。按照系統(tǒng)的存儲器配置，把正弦表第一項“sine_table”放在0x0C00的位置，余弦表第一項“cos-table”放在0x0D00的位置。

在計算蝶形輸出時相應(yīng)的實部是用減法，虛部用加法，即

因此，蝶形結(jié)第一個輸出數(shù)據(jù)的實部表示為PR′＝PR＋QR×WR+QI×WI，第一個輸出數(shù)據(jù)的虛部表示為PI′=PI+QI×WR-QR×WI。同樣可得到蝶形結(jié)第二個輸出數(shù)據(jù)的表示形式，這里不再給出[6]。相應(yīng)的關(guān)鍵匯編指令如下：

其中，第二個數(shù)據(jù)與系數(shù)相乘后，所得結(jié)果是32位數(shù)，與第一個數(shù)據(jù)的實部相加時，需要將該實部數(shù)據(jù)左移16位后執(zhí)行；并行指令可以在一個指令周期內(nèi)，將B高端移位ASM，存于AR2所指定的存儲單元中，并執(zhí)行將AR2所指的單元中的值左移16位后減去累加器B的值；在實際編程時，為了避免溢出，對每一步的輸出右移一位（相當(dāng)于將結(jié)果除以2，相當(dāng)于把的比例因子分散到各級運(yùn)算之中，可以在保持輸出信號方差不變、輸出信號最大幅度等于輸入值N倍的情況下，減小輸出噪聲電平，增加輸入信號幅度）。應(yīng)該指出，在FFT的各級運(yùn)算過程中，相應(yīng)的中間結(jié)果、FFT輸出和整理后的結(jié)果都存儲在該緩沖區(qū)中，從而實現(xiàn)原位計算。

第三，產(chǎn)生輸出結(jié)果。這里將D（k）重新表示為實部和虛部的形式[7]：

代入式(6)可以得到2N點實數(shù)序列的DFT表示式：

A（k）=AR（k）+jAI（k）

從而有

其中，RP（k）和 IP（k）分別表示 D（k）的偶對稱部分的實部和虛部，RM（k）和 IM（k）分別表示 D（k）的奇對稱部分的實部和虛部，即

編程時首先實現(xiàn)輸出復(fù)數(shù)序列的奇偶分解，求解RP(k)、IP(k)、RM(k)和IM(k)。這一點通過4個指針分別指向R(k)、R(N-k)、I(k)和I(N-k)，然后完成式(10)的計算。程序最后按照原始序列的DFT對上述結(jié)果按照式(9)計算，并以自然順序存儲在整個4N長的數(shù)據(jù)緩沖區(qū)內(nèi)。

3 實驗結(jié)果與分析

給定包含兩個頻率分別為5Hz、10Hz正弦信號，采樣頻率200Hz，取數(shù)據(jù)長度為256。將Matlab生成的信號源文件轉(zhuǎn)換為CCS能夠識別的.dat格式，這里命名為input.dat文件。在CCS中，首先設(shè)置探針斷點Toggle Probe point，從File菜單項中選擇File I/O，在出現(xiàn)的窗口中選擇輸入文件input.dat。單擊Add Probe Point，出現(xiàn)探針斷點窗口下，將輸入文件與設(shè)置的探測點連接起來。根據(jù)公式(3)將a(n)組合為復(fù)數(shù)序列d(n)，經(jīng)過128點的FFT后得到其傅里葉變換D(k)，并順序存放在0x2200h-0x23ffh的緩沖區(qū)中，前半部分為d(n)、后半部分為D(k)，如圖2。由于復(fù)數(shù)的實部和虛部分開存放，因此一共占用512個字長。圖3是256點實數(shù)序列a(n)的傅里葉變換A(k)，頻譜關(guān)于零頻對稱，相應(yīng)的功率譜顯示在圖4中，與A(k)不同的是，功率譜僅有256點。為了比較起見，圖5給出了利用Matlab軟件計算得到的實數(shù)序列a(n)的功率譜結(jié)果，與圖4相比較，可知DSP算法匯編程序給出了信號正確的傅里葉變換結(jié)果。

圖2 D(k)和輸入信號

圖3 實數(shù)序列的FFT

圖4 信號功率譜（ccs）

圖5 信號功率譜（Matlab）

4 結(jié)論

2N點實數(shù)序列可以用N點FFT算法實現(xiàn)，實驗證明，在TMS320VC5416定點DSP中運(yùn)行結(jié)果正確。同樣該程序經(jīng)過簡單的參數(shù)修改可以適用于其它點數(shù)的信號傅里葉變換，并應(yīng)用于實際系統(tǒng)中。這種從理論到實現(xiàn)的流程教學(xué)方法，不僅增強(qiáng)了實際的教學(xué)效果，也有助于培養(yǎng)學(xué)生的工程實踐和創(chuàng)新能力。

[1]吳鎮(zhèn)揚(yáng).數(shù)字信號處理[M].北京:高等教育出版社,2004.

[2]陳佩青.數(shù)字信號處理教程[M].北京:清華大學(xué)出版社,2007.

[3]趙宏怡.DSP技術(shù)與應(yīng)用實例[M].北京:電子工業(yè)出版社,2008.

[4]陳金鷹.DSP技術(shù)及應(yīng)用[M].北京:機(jī)械工業(yè)出版社,2011.

[5]劉和平.TMS320C240xDSP C語言開發(fā)應(yīng)用[M].北京:北京航空航天大學(xué)出版社,2003.

[6]邱立存,聞武,劉海英.TMS320C54X系列DSP上FFT運(yùn)算的實現(xiàn)[J].微計算機(jī)信息,2005(7):136-137.

[7]肖宛昂.嵌入式系統(tǒng)中FFT算法研究[J].單片機(jī)與嵌入式系統(tǒng)應(yīng)用,2003(4):68-69.