利用對稱性加速實(shí)序列ＦＦＴ的方法及其ＦＰＧＡ實(shí)現(xiàn)
2009-01-01 00:00:00鄧宏貴郭晟偉
計(jì)算機(jī)應(yīng)用研究 2009年1期
(中南大學(xué) 物理學(xué)院 電信系， 長沙 410083)
摘 要：針對工程實(shí)踐中傅里葉變換的輸入序列一般為實(shí)序列的情況，充分利用FFT（快速傅里葉變換）奇偶虛實(shí)的對稱性質(zhì)，提出了一種實(shí)序列FFT的加速算法。將2N點(diǎn)的實(shí)序列DFT轉(zhuǎn)換為N點(diǎn)的復(fù)序列DFT，并行計(jì)算使運(yùn)算量明顯減少；并給出了基于FPGA的硬件實(shí)現(xiàn)方法。
關(guān)鍵詞：快速傅里葉變換； 實(shí)序列； 現(xiàn)場可編程門陣列
中圖分類號：TP301.6 文獻(xiàn)標(biāo)志碼：A
文章編號：10013695(2009)01009202
Using symmetry accelerating real FFT and its FPGA implementation
DENG Honggui， GUO Shengwei
(Dept. of Electronics Information， School of Physics Science Technology， Central South University， Changsha 410083， China)
Abstract:Regarding the input of FFT being generally real sequence in engineering practices，making full use of symmetrical property，this papar proposed a improvement algorithm for real FFT. Real FFT computation of 2N points were transformed to complex FFT computation of N points. The amount of calculation could significantly reduce while parallel calculating. And gave a hardware implementation on FPGA.
Key words：fast Fourier transform(FFT); real sequence; FPGA
離散傅里葉變換（DFT）是數(shù)字信號處理領(lǐng)域的一個(gè)重要分析工具，其作用是將信號從時(shí)域轉(zhuǎn)換到頻域，因而DFT是進(jìn)行濾波、頻譜細(xì)化、譜估計(jì)等頻域處理的基礎(chǔ)，在聲學(xué)、電學(xué)、圖像等工程實(shí)踐領(lǐng)域有著不可替代的作用[1~3]。Cooley和Tukey于1965年提出了傅里葉變換的快速算法之后，使N點(diǎn)DFT的乘法計(jì)算量由N2次降為（N/2） log2 N次，極大促進(jìn)了DFT的應(yīng)用。然而CooleyTukey算法只考慮了輸入序列為復(fù)序列的情況，若為實(shí)序列，將虛部補(bǔ)零當(dāng)做復(fù)序列處理，這樣必然影響算法效率[4，5]。實(shí)際上，工程實(shí)踐中輸入序列通常都是實(shí)序列[6，7]。本文充分利用傅里葉變換的對稱性和周期性，提出了一種針對實(shí)序列的改進(jìn)算法，通過將原2N點(diǎn)實(shí)序列奇偶分離，構(gòu)造一個(gè)N點(diǎn)復(fù)序列，把一個(gè)2N點(diǎn)的實(shí)序列DFT計(jì)算轉(zhuǎn)換為一個(gè)N點(diǎn)的復(fù)序列DFT計(jì)算，然后將N點(diǎn)復(fù)序列的輸出序列進(jìn)行適當(dāng)?shù)倪\(yùn)算組合，獲得原2N點(diǎn)實(shí)序列的DFT輸出序列。這樣在利用FPGA并行計(jì)算的時(shí)候，使計(jì)算量減少了將近一半，運(yùn)算效率提高了將近一倍，更進(jìn)一步地滿足了信號實(shí)時(shí)處理的要求。
1 實(shí)序列FFT傳統(tǒng)算法原理
DFT是將時(shí)域序列轉(zhuǎn)換為長度相同的頻域序列，從理論上講，默認(rèn)輸入序列為復(fù)序列。長度為N點(diǎn)的有限長序列x(n)的DFT定義為
x(k)=∑N-1n=0x(n)WnkN k=0，1，…，N-1(1)
其中WN=e-j(2π/N）為旋轉(zhuǎn)因子。
輸入序列x(n)當(dāng)做復(fù)序列處理，若實(shí)際應(yīng)用中為實(shí)序列，則將虛部補(bǔ)零。設(shè)x(n)是一2N點(diǎn)的實(shí)序列，將其進(jìn)行奇偶分離，其偶序號依次構(gòu)成的新序列為g(n)，奇序號依次構(gòu)成的新序列為h(n)，則g(n)、h(n)仍為N點(diǎn)實(shí)序列。記g(n)、h(n)的傅里葉變換分別為G(k)、H(k)。根據(jù)CooleyTukey時(shí)間抽取基2算法：
X(k)=G(k)+WKNH(K) k=0，1，…，N-1(2)
X(k+N)=G(k)-WKNH(k) k=0，1，…，N-1(3)
在計(jì)算N點(diǎn)實(shí)序列g(shù)(n)、h(n)的傅里葉變換G(k)、H(k)時(shí)，繼續(xù)按照上面的方法進(jìn)行拆分，直至兩點(diǎn)的DFT為止。在計(jì)算時(shí)，則從兩點(diǎn)的DFT開始，逐級計(jì)算出原始序列的DFT。該算法乘法計(jì)算量為N log2 2N次（2N點(diǎn)序列）。
2 實(shí)序列FFT加速算法原理
通過觀察式(2)(3)可知，經(jīng)過一次分解之后，需計(jì)算兩個(gè)N點(diǎn)實(shí)序列g(shù)(n)、h(n)的傅里葉變換G(k)、H(k)?？梢猿浞掷酶道锶~變換的周期性和對稱性，將兩個(gè)N點(diǎn)的實(shí)序列的DFT轉(zhuǎn)換為一個(gè)N點(diǎn)復(fù)序列DFT一次性求得。以下給出該算法的原理及理論依據(jù)。
構(gòu)造y(n)=g(n)+i h(n)，則y(n)為一N點(diǎn)復(fù)序列。y(n)的傅里葉變換Y(k)為
Y(k)=∑N-1n=0(g(n)+i h(n))WnkN k=0，1，…，N-1(4)
用N-k代換k得：
Y(N-k)=∑N-1n=0(g(n)+i h(n))Wn(N-k)N k=0，1，…，N-1(5)
由旋轉(zhuǎn)因子的周期性可知:Wn(N-k)N=W-nkN
式(5)轉(zhuǎn)換為
Y(N-k)=∑N-1n=0(g(n)+i h(n))W-nkN k=0，1，…，N-1(6)
記G(k)的實(shí)部和虛部分別為Gr(k)和Gi(k)，H(k)的實(shí)部和虛部分別為Hr(k)和Hi(k)。
G(k)=∑N-1n=0g(n)WnkN=∑N-1n=0g(n)(cos θ+i sin θ)
k=0，1，…，N-1;θ=-2πnk/N（7）
H(k)=∑N-1n=0h(n)WnkN=∑N-1n=0h(n)(
cos θ+i sin θ)
k=0，1，…，N-1;θ=-2πnk/N（8）
因?yàn)間(n)為實(shí)序列，Gr(k)=∑N-1n=0g(n)cos θ；Gi(k)=∑N-1n=0g(n)sin θ。式(4)+(6)得：
Y(k)+Y(N-k)=2∑N-1n=0(g(n)+i h(n))cos θ θ=-2πnk/N(9)
兩邊取實(shí)部得：
Yr(k)+Yr(N-k)=2∑N-1n=0g(n)cos θ=2Gr(k)(10)
可得：Gr(k)=(Yr(k)+Yr(N-k))/2(11)
用類似的方法可以得到其他幾個(gè)方程：
Gi(k)=(Yi(k)-Yi(N-k))/2(12)
Hr(k)=(Yi(k)+Yi(N-k))/2(13)
Hi(k)=-(Yr(k)-Yr(N-k))/2(14)
因而可以由Y(k)分離出G(k)、H(k)：
G(k)=(Yr(k)+Yr(N-k))/2+i(Yi(k)-Yi(N-k))/2(15)
H(k)=(Yi(k)+Yi(N-k))/2-i(Yr(k)-Yr(N-k))/2(16)
這樣只要求出N點(diǎn)復(fù)序列y(n)的傅里葉變換Y(k)，即可分離出g(n)、h(n)的傅里葉變換G(k)、H(k)；再代入式（2）（3），即可得到2N點(diǎn)的實(shí)序列x(n)的傅里葉變換X(k)。而使用CooleyTukey算法求N點(diǎn)復(fù)序列的DFT時(shí)，復(fù)雜度與求N點(diǎn)實(shí)序列DFT基本相當(dāng)，遠(yuǎn)低于原2N點(diǎn)實(shí)序列DFT的復(fù)雜度，運(yùn)算效率顯著提高。
3 實(shí)序列FFT算法FPGA實(shí)現(xiàn)
FFT算法的硬件實(shí)現(xiàn)一般有DSP和FPGA兩種方案。DSP速度較慢，接口不靈活，而且沒有FFT運(yùn)算所需要的巨量存儲(chǔ)器，需外置特定的接口控制芯片和RAM，限制了運(yùn)算速度；但DSP開發(fā)相對簡單、技術(shù)成熟、開發(fā)費(fèi)用相對較低，目前大部分FFT硬件仍是用DSP來實(shí)現(xiàn)。近年來，隨著微電子工藝和FPGA技術(shù)的成熟，F(xiàn)PGA的速度和容量迅速提高，且內(nèi)嵌硬件乘法器。由于其體積、速度、靈活性等各種性能都優(yōu)于DSP，成為了FFT算法硬件實(shí)現(xiàn)的理想選擇。本文的實(shí)序列FFT加速算法就是用Altera公司Stratix系列FPGA器件EP1S20F484C5實(shí)現(xiàn)的。
3.1 FFT實(shí)序列處理器硬件框圖
基于FPGA設(shè)計(jì)的FFT處理器主要由輸入數(shù)據(jù)存儲(chǔ)器1、輸出數(shù)據(jù)存儲(chǔ)器2、運(yùn)算模塊、控制模塊、ROM因子表等構(gòu)成，如圖1所示。此外，由于本設(shè)計(jì)是針對實(shí)輸入序列，必須先將2N點(diǎn)的實(shí)序列轉(zhuǎn)換為N點(diǎn)的復(fù)序列，再經(jīng)FFT變換，最后需要進(jìn)行奇偶分離還原出原2N點(diǎn)實(shí)序列的傅里葉變換。
3.2 并行處理器及同址運(yùn)算
根據(jù)CooleyTukey時(shí)間抽取基2算法，計(jì)算一個(gè)N點(diǎn)的FFT，需要log2 N級碟形運(yùn)算，每級由N/2個(gè)碟形單元構(gòu)成。各級順序執(zhí)行，而同一級內(nèi)的N/2個(gè)碟形單元相互獨(dú)立，可以并行計(jì)算。這樣整個(gè)設(shè)計(jì)需要N/2個(gè)碟形單元，既不浪費(fèi)FPGA的內(nèi)部資源，又以并行方式提高了運(yùn)算速度。同時(shí)，由于每一個(gè)碟形單元的輸入不再參與別的碟形單元運(yùn)算，基于這個(gè)特點(diǎn)可以進(jìn)行同址運(yùn)算，即將其輸出數(shù)據(jù)仍放到輸入數(shù)據(jù)的存儲(chǔ)單元中，既節(jié)省了存儲(chǔ)單元，又便于尋址。
3.3 W因子的生成
在FFT中，旋轉(zhuǎn)因子WnkN=cos(2πnk/N)-i sin(2πnk/N)。 產(chǎn)生旋轉(zhuǎn)因子一般有兩種方案：a)直接計(jì)算產(chǎn)生，每次調(diào)用都臨時(shí)計(jì)算旋轉(zhuǎn)因子的值，其優(yōu)點(diǎn)是節(jié)省ROM存儲(chǔ)單元;b)將W因子制表。由于WnkN只有N個(gè)獨(dú)立的值，可以預(yù)先計(jì)算出WnkN的值，存于FPGA的ROM中，調(diào)用時(shí)可直接查表得到，這樣提高了運(yùn)算速度。在一些信號實(shí)時(shí)處理系統(tǒng)中，F(xiàn)FT的運(yùn)算速度是最重要的。而FPGA中的ROM能夠提供足夠的存儲(chǔ)單元來存放W因子表。
3.4 復(fù)數(shù)乘法單元及其改進(jìn)
輸入數(shù)據(jù)x=xr+j xi與旋轉(zhuǎn)因子W=cos θ+j sin θ相乘的公式為
yr=xrcos θ-xisin θ (17)
yi=xicos θ+xrsin θ(18)
一般而言，在一個(gè)周期內(nèi)實(shí)現(xiàn)一次復(fù)數(shù)乘法需要四個(gè)實(shí)數(shù)乘法器和兩個(gè)實(shí)數(shù)加法器。本文對式(17)(18)加以改進(jìn)可節(jié)省運(yùn)算量。
yr=(xr+xi)cos θ-xi(sin θ+cos θ)(19)
yi=(xr+xi)cos θ+xr(sin θ-cos θ)(20)
可以把cos θ、sin θ+cos θ、sin θ-cos θ三個(gè)值預(yù)先計(jì)算出來存儲(chǔ)在ROM中，因而改進(jìn)之后，實(shí)現(xiàn)一次普通的復(fù)數(shù)乘法只需要三次實(shí)數(shù)乘法和三次實(shí)數(shù)加法。由于一個(gè)實(shí)數(shù)加法器占用的面積遠(yuǎn)小于一個(gè)實(shí)數(shù)乘法器，加法器的運(yùn)算速度也高于乘法器，面積和速度方面都得到了優(yōu)化。
4 性能比較
表1為幾種FFT硬件實(shí)現(xiàn)的性能比較，使用Stratix系列FPGA器件EP1S20F484C5設(shè)計(jì)的并行流水線蝶形運(yùn)算單元能夠穩(wěn)定地運(yùn)行于200 MHz時(shí)鐘下。對于N點(diǎn)FFT需要log2 N級，每級N/2個(gè)碟形單元。以1 024點(diǎn)實(shí)序列為例，由于本設(shè)計(jì)采用奇偶分離的改進(jìn)算法，在預(yù)處理階段將1 024點(diǎn)實(shí)序列FFT轉(zhuǎn)換為512點(diǎn)復(fù)序列FFT，共有九級，每級256次碟形運(yùn)算。
完成一次1 024點(diǎn)實(shí)序列FFT所需要的時(shí)間是
T=256×9×5ns=11.520 μs
5 結(jié)束語
根據(jù)DFT的理論定義，默認(rèn)輸入數(shù)據(jù)為復(fù)序列；而在工程實(shí)踐中，DFT的輸入數(shù)據(jù)一般為實(shí)序列。本文提出的加速實(shí)序列FFT的方法在FPGA中并行計(jì)算，節(jié)省了約一半的乘法運(yùn)算量，提高了運(yùn)算速度，具有很大的工程實(shí)踐應(yīng)用意義。
參考文獻(xiàn)：
［1］劉粵鉗，姚紅玉.分?jǐn)?shù)傅里葉變換在信號處理中的應(yīng)用研究[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，2007，28(3):198214.
[2]劉莉，周樸.基于分?jǐn)?shù)傅里葉變換的圖像隱藏技術(shù)[J].國防科技大學(xué)學(xué)報(bào)，2005，27(6):6771.
[3]CHEN He，ZHAO Zhongwu.ASIC design of floatingpoint FFT processor[J].北京理工大學(xué)學(xué)報(bào)：英文版，2004，13(4):389393.
[4]魯劍鋒.采用雙DSP的實(shí)時(shí)傅里葉變換系統(tǒng)設(shè)計(jì)[J].光電工程，2005，32(11): 9396.
[5]陳恒量，蔣勇.基于DSP的實(shí)數(shù)FFT算法研究與實(shí)現(xiàn)[J].動(dòng)力學(xué)與控制學(xué)報(bào)，2005，3(2):5053.
[6]GAN Liangcai，DING Yahui. CPLD implementation of 1024point FFT of shortwave wideband channel simulator[J].電子學(xué)報(bào)：英文版，2006，15(4):692696.
[7]劉國棟，陳伯孝，陳多芳.FFT處理器的FPGA設(shè)計(jì)[J].航空計(jì)算技術(shù)，2004，34(3): 101104.
[8]KUO J C，WEN C H，WU A Y.Implementation of a programmable 642048point FFT/IFFT processor for OFDM based communication systems[C]//Proc ofIEEE International Symposium on Circuits and Systems. 2003:121124.
[9]榮瑜，朱恩.一種高性能FFT蝶形運(yùn)算單元的設(shè)計(jì)[J].東南大學(xué)學(xué)報(bào)，2007，37(4):565568.
