基于ADSP-TS201的大點數(shù)FFT實現(xiàn)

穆文爭 王啟智 朱子平

(中國電子科技集團公司第三十八研究所 合肥 230088)

0 引言

在雷達信號處理中，F(xiàn)FT作為是頻域脈壓和譜分析的基礎，是較為常用的一種變換，一般由數(shù)字信號處理器完成。在工程應用中，當需要做FFT變換時，設計師可以調(diào)用開發(fā)環(huán)境自帶的庫函數(shù)，也可以將FFT做成公用基礎模塊，供設計者調(diào)用。此類模塊一般采用匯編語言設計，針對所使用處理器的結構特征進行優(yōu)化，具有較高的運算效率。但是當需要變換的FFT點數(shù)較大時，該類模塊不再適用，因為輸入數(shù)據(jù)和中間緩存都要占據(jù)很大的內(nèi)存空間，而一般DSP的內(nèi)存空間有限且被劃分為多個不連續(xù)的block，這就破壞了該類模塊的并行運算結構，降低了運算效率，如果采用與外部存儲器交換數(shù)據(jù)的辦法，由于與外存交換數(shù)據(jù)的速度較慢，也會降低運算效率。因此，在工程上需要做大點數(shù)FFT的場合，就要考慮如何實現(xiàn)的問題。

本文首先介紹了大點數(shù)FFT變換的數(shù)學原理，即將一維FFT拆成二維FFT實現(xiàn)的算法原理，然后以96K點為例，介紹其在通用處理器ADSP-TS201上的實現(xiàn)過程。

1 數(shù)學原理

設有限長序列x(n)的長度為N，則其DFT為

令N=N1N2，可以將x(n)分解為N2個長度為N1的序列，將這些序列用陣列x'表示，則

令n和k的序號映射定義如下:

則N點DFT可以表示為:

令

則G(n2，k1)為序列x(N2n1+n2)的N1點 DFT，即上式二維陣列N2行元素的DFT。計算陣列x'每一行N1點DFT就得到另一個矩陣:

矩陣的元素為復數(shù)G(n2，k1)。因n2行的數(shù)據(jù)在計算完x(N2n1+n2)的N1點DFT后不再需要，所以G(n2，k1)可以存儲在同一行。計算X(k)，還應乘以旋轉(zhuǎn)因子Wk1n2N形成新的陣列

這樣DFT變換結果就可由二維陣列得出。

由此可見，將一維FFT拆成二維FFT的實現(xiàn)過程就是將N點FFT分解為N2個N1點行變換與N1個N2點列變換的級聯(lián)［2］，經(jīng)過二維拆分以后，設計師就可以調(diào)用已有的FFT函數(shù)。

2 工程實現(xiàn)面臨的問題

除了輸入數(shù)據(jù)和中間緩存要占據(jù)很大的內(nèi)存空間外，做大點數(shù)FFT遇到的另外一個問題就是存儲旋轉(zhuǎn)因子。當變換的點數(shù)很大時，在FPGA或DSP中存儲這些旋轉(zhuǎn)因子需要很大的內(nèi)存空間，這在工程實現(xiàn)上不方便。下面說明如何對旋轉(zhuǎn)因子進行分解，以降低對內(nèi)存的需求量。

如果將N分解為N1×N2，k分解為N1×α+β，其中α為倍數(shù)，β為余數(shù)。則有

因為k＜N，即N1·α+β＜N1·N2，又β是k對 N1的余數(shù)，有0≤β≤N1-1 ，所以 N1·α ＜ N1·N2，即α＜N2。

以1M點FFT為例，N可以分解為N1×N2=1024×1024，這樣存儲的旋轉(zhuǎn)因子個數(shù)就由1M點減少為1024+1024=2048點。

3 96K點FFT工程實現(xiàn)

以上對大點數(shù)FFT的實現(xiàn)方法做了原理性介紹，下面通過一個具體的例子，說明其在通用處理器ADSP-TS201上的實現(xiàn)過程。

由上述一維FFT拆成二維FFT的算法原理，可以歸納出大點數(shù)FFT實現(xiàn)的基本步驟［3］:

a.將N點輸入數(shù)據(jù)表示成L行C列矩陣;

b.對矩陣所有列分別做L點FFT，得到G(l1，c);

c.將G(l1，c)乘以中間旋轉(zhuǎn)因子得到GG(l1，c);

d.對GG(l1，c)所有行做C點FFT，得到X(c+l1×C);

e.將X(c+l1×C)整序為X(k)，變換完成。

圖1為大點數(shù)FFT實現(xiàn)流程圖。

需要注意的是，數(shù)據(jù)輸入分解和因計算中間旋轉(zhuǎn)因子而進行的分解是兩種不同的分解，它們可以不同。例如本例中，數(shù)據(jù)輸入可以分解為96k=3×32k，計算中間旋轉(zhuǎn)因子可以分解為96k=256×384。在對96k=N1×N2分解時，需要遵循一個原則，就是其中一個因子要為2的整冪次，因為這樣便于利用ADSP-TS201的匯編指令FEXT快速提取α和β。

可知，旋轉(zhuǎn)因子只需存儲一半即可，另一半可由共軛對稱性得到;在存儲旋轉(zhuǎn)因子時，由于 β 的范圍是0～(N1-1)，所以這N1個旋轉(zhuǎn)因子要全部存儲。這樣需要存儲的旋轉(zhuǎn)因子個數(shù)就可以由(N1+N2)點減少為點，從而進一步減少了所需的內(nèi)存空間。有了這兩種旋轉(zhuǎn)因子，計算大點數(shù)FFT時任何一種旋轉(zhuǎn)因子都可以通過(9)式實時生成。

表1是96K點FFT的N1和N2不同分解方式對應的旋轉(zhuǎn)因子存儲量，可見256×384是存儲量最小的一種分解方式。

表1 不同分解方式對應的旋轉(zhuǎn)因子存儲量

由于所使用的處理器為ADSP-TS201，下面對該處理器的內(nèi)存加以說明。ADSP-TS201具有24Mbit內(nèi)存，被劃分為 6個 block，每個 block為4Mbit(128K ×32bit)，組織方式如下［4］。

在采用開發(fā)環(huán)境(VisualDSP++5.0)默認的鏈接描述文件時，block0用來存放代碼，block4的部分空間用作堆棧，這樣可供用戶使用的完整的128K字內(nèi)存塊就只有4塊，即 block2、block6、block8和block10。

處理96K點FFT時，由于一個block最多存儲64K點復數(shù)，所以需要1.5個block作為輸入數(shù)據(jù)緩存?？梢钥紤]將96k點復數(shù)分為三段、每段32K點進行存儲，其中一種存儲方式如圖3所示。

圖2 ADSP-TS201的內(nèi)存組織

圖3 96K點輸入數(shù)據(jù)的緩存方式

下面按照前述的實現(xiàn)步驟，對96K點復數(shù)進行FFT變換。

第一步:

對輸入的3段數(shù)據(jù)，先對每列做3點DFT，然后對做過DFT的數(shù)據(jù)乘以中間旋轉(zhuǎn)因子。需要注意的是，每列3點數(shù)據(jù)所乘的旋轉(zhuǎn)因子均不同，旋轉(zhuǎn)因子的上標為行號×列號。乘積結果覆蓋原數(shù)據(jù)見圖4。

圖4 每列3點DFT

每列旋轉(zhuǎn)因子的上標見表2。

表2 旋轉(zhuǎn)因子上標

第二步:

數(shù)據(jù)經(jīng)過第一步處理后，再對每一行做FFT，結果仍然覆蓋原數(shù)據(jù)，如圖5所示。

圖5 每行做FFT

第三步:

將結果重新排序。由于DSP沒有足夠的內(nèi)存用作緩存，所以重排只能分段進行，這樣重排就需要兩輪，重排示意圖見圖6。

圖6 數(shù)據(jù)重排

第一輪排序:

申請32K空間作為重排的緩存，對這32K緩存，只用到其中的24K。數(shù)據(jù)重排的程序由匯編指令設計，為了提高效率，存取數(shù)據(jù)同時進行。先取數(shù)據(jù)至寄存器，存的時候以需要的順序從寄存器讀出，這樣將數(shù)據(jù)存至緩存的時候，就已經(jīng)重排好了。存完24K數(shù)據(jù)(每段取8K)后，再將數(shù)據(jù)搬回，將原數(shù)據(jù)覆蓋，這一過程如圖7。

經(jīng)過第一輪重排后，數(shù)據(jù)分布規(guī)律如圖8所示。

經(jīng)過第一輪重排后，每段數(shù)據(jù)的首地址如圖9所示。

可見，經(jīng)過第一輪重排后，數(shù)據(jù)雖然得到了調(diào)整，但沒有達到圖6所示的順序，還要經(jīng)過第二輪重排，才能完成整個排序。

第二輪排序:

第二輪排序就是對圖8的12個數(shù)據(jù)段進行交換，直至達到圖6的順序。交換時以整個數(shù)據(jù)段進行，第二輪排序的交換順序如圖10所示。

第二輪排序時，每次交換的2個數(shù)據(jù)段的首地址見表3。

圖7 第一輪重排示意圖

圖8 第一輪重排后的數(shù)據(jù)分布

圖9 第一輪重排后每段數(shù)據(jù)的首地址

圖10 第二輪排序數(shù)據(jù)交換過程

表3 第二輪排序時每次交換的2個數(shù)據(jù)段的首地址

至此，經(jīng)過第二輪排序以后，數(shù)據(jù)在內(nèi)存中的順序就達到了期望的順序，即完成了整個排序，從而完成了整個FFT變換。

對于大點數(shù)IFFT，可以在FFT的基礎上實現(xiàn)，步驟如圖11所示。

圖11 大點數(shù)IFFT實現(xiàn)步驟

運行時間統(tǒng)計:在工程實現(xiàn)時，96K點復數(shù)FFT模塊由四個函數(shù)組成:a.對每列數(shù)據(jù)做3點DFT并乘以中間旋轉(zhuǎn)因子的函數(shù)DFT_3point_multMidT-wid.asm;b.對每行做32K點FFT的函數(shù)bw_cvfft_c.asm，該函數(shù)為通用的參數(shù)化FFT模塊;c.第一輪排序函數(shù)ReArrange_dat_1stRound.asm;d.第二輪排序函數(shù)ReArrange_dat_2ndRound.asm。當ADSPTS201的核時鐘運行于480MHz時，每個函數(shù)模塊的運行時間統(tǒng)計見表4。

表4 函數(shù)運行時間統(tǒng)計

可見運行完96K點復數(shù)FFT需要的時間為:1.66+1.6*3+0.35+0.41=7.22ms。經(jīng)測試，做完96K點頻域脈壓并對結果求模取對數(shù)，共計用時16.7ms?？梢姡@種處理大點數(shù)FFT的方法不僅解決了內(nèi)存需求量大的問題，而且函數(shù)模塊具有較好的實時性。

4 結束語

針對雷達信號處理遇到的大時寬脈壓問題，給出了大點數(shù)FFT的變換方法。該方法大大降低了因存儲旋轉(zhuǎn)因子對內(nèi)存的需求，在不與外部存儲器交換數(shù)據(jù)的情況下，能夠快速完成FFT變換。文章首先介紹了大點數(shù)FFT變換的數(shù)學原理，然后以96K點為例，闡述了其在通用處理器ADSP-TS201上的實現(xiàn)過程，為工程上實現(xiàn)其它點數(shù)的大點數(shù)FFT提供了有益的借鑒。

［1］丁玉美，高西全.數(shù)字信號處理［M］.西安:西安電子科技大學出版社，2000:68-77.

［2］蘇濤，莊德靖.大點數(shù)FFT算法的改進及其實現(xiàn)［J］. 現(xiàn)代雷達，2005，27(7):23-26.

［3］王曉君，龍騰，周希元.二維級聯(lián)流水結構大點數(shù)FFT運算器實現(xiàn)研究［J］.無線電工程，2010，40(11):19-22.

［4］ 劉書明，羅勇江.ADSP TS20XS系列DSP原理與應用設計［M］.北京:電子工業(yè)出版社，2007:115-118.