一種高速2-D滑動FFT的設計實現

許丁鴻, 張多利, 陶相穎, 韓帥鵬, 宋宇鯤

(1.合肥工業大學 微電子學院,安徽 合肥 230601; 2.教育部IC設計網上合作研究中心,安徽 合肥 230601)

在波形分解[1]、醫學成像[2]、雷達[3]以及通信系統[4]等諸多應用中,都采用快速傅里葉變換(fast Fourier transform,FFT)算法完成信號的時域與頻域間變換。在遠程醫療、口語互譯等實時性應用領域,受傳統FFT算法對序列完整性依賴的約束,無法在序列采樣過程中開展FFT變換,變換的實時性難以達到目標應用要求。此外,在處理時域數據更新周期較短的長序列信號時,受運算復雜度和計算實時性約束,傳統FFT只能依靠大量計算器資源達到設計目標,導致整個硬件結構異常復雜[5]。

針對此類問題,文獻[6]使用滑動FFT算法。滑動FFT[7]的執行過程為迭代累加累乘計算,可在前一次滑動FFT結果基礎上遞推得到當前滑動FFT結果,復用前次運算結果,節約運算時間,提高信號分析實時性。利用滑動FFT運算遞推特征,可以根據信號分析需要選擇計算相應的頻譜區段,增加信號分析效果。

本文設計的滑動FFT變換器采用2-D FFT[8]硬件架構,即將1-D長序列FFT分解成由若干條短序列組成的矩陣,再分別對行/列向量進行1-D FFT運算,利用行/列向量間無關性實現多路并行。FFT處理器通常分為順序處理結構、級聯處理結構[9]、并行迭代結構[10]、陣列結構4類,本文根據應用需求選擇陣列結構。

1 算法原理

1.1 滑動FFT原理

設時域信號波形的采樣值為x(n),傅里葉變換截取的滑窗的長度為N,滑窗每次滑動的步長為p=2L。記i點起連續采樣N點數據的傅立葉變換為Xi(k),則i+p點至i+p+N點的N點傅立葉變換記為Xi+p(k)。

根據離散傅里葉變換(discrete Fourier transform,DFT)的公式假設,即

k=0,1,2,…,N-1

(1)

i+p點的FFT公式可以表示為:

(2)

(3)

令

(4)

式(4)中e(n)的前p=2L點為xi(n+N)-xi(n),后N-p點為0。

文獻[11]提出采用Pruning算法解決上述特征序列FFT計算問題。由式(3)可得,當窗長N=2M,步長p=2L時,僅需L級的FFT運算即可得出N點的FFT值。前M-L級的運算都是輸入的p=2L序列與0進行蝶形運算,也就是前p=2L點數據的復制,并不需要用到加法和乘法運算。N=16、L=2的Pruning算法圖如圖1所示。圖1中:虛線表示0值;實線表示輸入序列值。

從圖1可以看出:第1級與第2級是對輸入序列與0值做蝶形運算,可視為對輸入值的復制;第3級與第4級是2個輸入值為非零的復數蝶形運算。

圖1 輸入值為4點的16點FFT Pruning算法圖

假設已知Xi(k)的N點FFT值,滑動步長為p=2L,計算Xi+p(k)的FFT值,采用滑動FFT算法和使用傳統FFT算法節約的時間比值tτ計算公式為:

(5)

滑動FFT算法需要進行NL/2+N次復數乘法、NL+N次復數加法,而傳統的FFT算法需要進行(NlbN)/2次復數乘法、NlbN次復數加法。滑動窗長度為64～256點,滑動步長為4～16點,滑動窗每滑動一次產生的復數乘法計算量與復數加法計算量的關系見表1所列。由表1可知,滑動FFT所需復數乘法數量和復數加法數量都小于傳統FFT算法。

表1 滑動窗長、滑動步長與運算量的關系

1.2 2-D FFT算法原理

2-D FFT算法[12]是一種優化1-D長序列FFT的運算方法,該算法將N點長序列FFT運算分解為兩級短序列FFT級聯運算。

對N點數據進行FFT運算,令N=LC,即將輸入的N點分解為L行、C列的矩陣形式。根據DFT算法公式,有

X(k)=X(Lk1+k0)=X(k1,k0)=DFT[x(n)]=DFT[x(Cn1+n0)]=DFT[x(n1,n0)]

(6)

時間域變量n可以表示為:

n=Cn1+n0

(7)

其中,n1=0,1,…,L-1;n0=0,1,…,C-1。

頻率域變量k可以表示為:

k=Lk1+k0

(8)

其中:k1=0,1,…,C-1;k0=0,1,…,L-1。

(9)

令

(10)

(11)

得:

(12)

其中:G(k0,n0)為L點的列FFT變換;X1(k0,n0)為C點行FFT變換。因為輸出序列X(k)的序列順序是滿足先填寫行方向再填寫列方向,和輸入數據相逆,所以還需要一次整序才是正確的FFT序列。

對一個N點的大點數進行FFT運算可以分解為C組L點一維列變換與L組C點的一維行變換的級聯,其實現步驟如下:

1) 將N點表示為L行、C列的矩陣。

2) 對矩陣的所有列進行L點FFT列變換得到G(k0,n0)。

4) 對H(k0,n0)矩陣進行C點的FFT行變換得到X1(k0,k1)。

5) 將X1(k0,k1)整序為X(k1,k0),然后輸出。

1.3 硬件架構選擇

本文采用2-D滑動FFT硬件結構,從以下方面進行考慮。

1) 算法性能。由式(3)可知,滑動窗由p個有效值與N-p個0值構成,以本文設計的參數滑動窗長為16、滑動步長為4舉例,使用2-D FFT結構可以構成第1行為滑動步長的有效值,其余3行都為0的二維矩陣,如圖2所示。

圖2 2-D滑動FFT的2-D運算圖

根據二維矩陣的運算順序,先執行列方向矯正計算,因為列構成為1個有效值與3個0值,該FFT結果等于有效值的列方向復制,所以可以忽略列方向FFT,直接對輸入數據進行列矯正。

為達到計算速度要求,本文設計行方向采用全展開FFT結構,對±1和±i的旋轉因子進行操作,采用交換虛實部的方法壓縮運算,對其他旋轉因子系數使用優化編碼方法壓縮運算。

2) 壓縮資源。首先計算過程中間數據存儲資源的消耗,1-D架構為全序列長度,2-D架構為短序列中最長的序列長度。在設計中短序列采用全展開結構,無需緩沖中間結果。

其次是壓縮浮點數乘法資源。滑動FFT運算中的矯正運算是影響設計速度的關鍵,以滑動窗長為256點、滑動步長為16為例,滑動FFT運算的矯正因子系數為±0.707 1、±0.382 7、±0.923 9。對上述因子使用高基Booth編碼操作將浮點數乘法轉變為常數乘法,優化邏輯結構,減少資源。

最后是公用因子優化復用。式(3)中滑動FFT運算的矯正因子系數為±0.707 1、 ±0.382 7、±0.923 9,也是16點FFT運算的旋轉因子系數。對滑動FFT矯正運算的優化方法可以在16點FFT運算中得到復用。

綜上所述,本文設計采用2-D滑動FFT結構使其可多路并行并減少電路復雜度及FFT運算時間,通過對基16核的優化和一部分定點數乘法的優化來減少資源消耗。

2 系統結構

本文設計的2-D FFT處理器結構如圖3所示。

圖3 2-D FFT處理器結構

處理器結構包括以下部分:① 控制單元模塊(control unit),用來調節各個模塊的工作內容以及狀態,對模塊置零的控制;② 數據預處理模塊(data preproccessing unit),完成式(3)中xi(n+N)-xi(n)的步驟以及對數據的倒序處理;③ FFT運算模塊(FFT calculate module),其中包含了輸入數據與2-D FFT的列方向矯正模塊(revolve unit)和基16的FFT運算單元(FFT computing unit);④ 滑動FFT的運算模塊(sliding FFT calculate),完成FFT運算模塊的結果與上一次滑動FFT運算結果進行相加的復數加法模塊(complex add),以及對該結果進行復乘調整操作的復數乘法模塊(complex mul);⑤ 存儲控制單元模塊(memory control unit),實現地址無沖突規則,管理存儲單元的訪存;⑥ 存儲單元組(memory unit),由多塊memory組成,負責存儲當前時刻滑動FFT的結果。

滑動FFT處理器的工作流程如圖4所示。

1) 滑動窗FFT接收到啟動信號后開始工作,單次滑動窗計算結束后判斷當前序列處理是否完成,否則讀取源數據并寫入先入先出(first input first output,FIFO)和預處理模塊中。

2) 當累計輸入數據的長度等于滑動窗長度時,在INPUT端口讀取xi(n+N),從FIFO讀取xi(n),在數據預處理模塊中完成xi(n+N)-xi(n)操作以及倒序操作,每16個周期從預處理模塊輸出16個浮點數。

3) 數據預處理模塊輸出的16個結果送入FFT列方向矯正模塊(revolve unit),每個周期完成16次與矯正因子復乘運算,得到16個結果;再將其結果分輸入FFT單元(FFT computing unit)。

4) FFT單元完成16點FFT計算。

5) FFT單元結果輸入滑動FFT運算模塊(sliding FFT calculate)進行矯正運算,即與前一次滑動FFT結果進行復數加法與復數乘法。

6) 上述結果分別寫入結果memory并輸出,同時作為前一次的FFT結果暫存在滑動FFT運算模塊,參與下一次的運算。

圖4 2-D滑動FFT處理器工作流程

2.1 基16蝶形運算器設計

設計的滑動FFT矯正模塊需要在每個周期內完成16點的復數加法與復數乘法,相應地,本文設計必須在每個周期輸出16點的FFT結果,因此采用基16核的全展開結構。

當旋轉因子系數為±1和±i時,可以通過交換符號位以及實部或虛部的方法規避復數乘法。剩余部分運算使用基2來構建基16核需要10次復乘與64次復加;使用基4來構成基16核需要8次復乘與96次復加,因此本文設計使用基2來構建基16的核,如圖5所示[13]。

圖5 16點蝶形單元結構

2.2 使用高基Booth編碼優化定點數乘法

本文設計的滑動FFT運算矯正模塊需要在16個周期內完成256點的復加與復乘,因此必須壓縮浮點數乘法器的關鍵路徑長度。

IEEE-754浮點數由1位符號位、8位階碼、23位尾碼構成,2個浮點數相乘是階碼相加,尾碼相乘。在本文設計中,16點FFT運算的旋轉因子系數與滑動FFT運算的矯正因子系數相同,僅有±0.707 1、±0.382 7、±0.923 9。因此可對系數0.707 1、 0.382 7、 0.923 9進行Booth編碼,使用常數與變量乘法替代通用浮點乘法,壓縮計算時間[14]。

本文設計使用高基Booth編碼進行系數優化[15],編碼的公式為:

xi+k+1=-2k+1,xi+k=2k,x(i-1)=20,k=0,1,…,a-2

(13)

其中,2a為Booth的基。0.707 1、0.382 9、0.923 9的高基編碼情況如圖6所示。

圖6中2的冪次方乘法可以通過移位來獲得,其余通過相加獲得(如3、6等),因為每個系數都可找到一個基獲得最小加法執行次數,所以本文設計中0.707 1、0.382 7使用基16編碼,0.923 9使用基32編碼,即上述3個系數的乘法操作分別轉換為8次加法操作,達到了資源最小化目標。

在TSMC 28 nm工藝下,高基Booth運算器模塊的DC綜合頻率在700 MHz以上,資源的消耗情況見表2所列。

(a) 0.382 7

表2 Booth編碼乘法器與DW乘法器的資源對比

2.3 滑動FFT的矯正模塊旋轉因子結構設計

本文設計的參數為窗長度256,滑動窗步長為16。式由(3)可知:

每16個值一次循環,對應著2-D FFT運算結果,矯正因子的排列方式如圖7所示。

圖7 滑動FFT矯正因子排列

由圖7可知: 1、5、9、13行的矯正因子的實部與虛部由±1構成;2、6、10、14行的矯正因子的實部與虛部由±0.382 7、±0.923 9構成;3、7、11、15行的矯正因子的實部與虛部由±0.707 1構成;4、8、12、16行的矯正因子的實部與虛部由±0.382 7、±0.923 9 構成。

設輸入的復數為a+bi,使a和b分別與定點數的常數0.707 1、0.382 7、0.923 9進行乘法計算,然后根據狀態機來進行帶有正負的實部與虛部的組合排列。

3 設計實現與性能分析

SMIC 28 nm工藝下,本文設計使用Synopsys公司Design Compiler和IC Compiler完成后端設計工作,設計的滑動窗FFT變換器工作主頻高于700 MHz,核心面積為1 980 μm×2 060 μm,版圖如圖8所示。

因為本文設計采用的滑動窗算法是一個乘法累加的情況,所以計算誤差隨滑動窗滑動次數增多而累加。運算周期性地清零[16]可以有效防止誤差累積,使設計的輸出穩定在一個誤差范圍內。

本文設計在Xilinx Virtex-7 FPGA芯片上進行功能和性能測試,輸入序列生成采用不同頻率的單頻正弦信號疊加隨機噪聲。

圖8 滑動窗FFT版圖

參考結果為MATLAB中FFT函數的輸出。

y1(t)=sin(2π×50×t×1/n)+5randn(size(t)),

y2(t)=sin(2π×100×t×1/n)+5randn(size(t)),

y3(t)=sin(2π×150×t×1/n)+5randn(size(t)),

y4(t)=sin(2π×200×t×1/n)+5randn(size(t))。

對以上4段信號進行n=1 024的采樣,采樣間隔為1/n。由以上4段正弦隨機信號拼接成一個4 096的不同頻率的序列數據。滑動窗長度256,滑動步長16,一共產生256段數據。在滑動了不同長度之后進行數據復位,保證誤差精度后的誤差分析,結果見表3所列。

表3 滑動FFT復位前10段結果數據的誤差分析

設計要求在600 MHz以上,使用Booth編碼對旋轉/矯正因子系數的常數化操作不僅減少了消耗資源,也提高了工作速度。600 MHz的工作頻率,Design Ware乘法器需要5級流水,而使用Booth編碼設計在單周期即可達到頻率要求。Booth編碼的優化系數也用于FFT計算模塊。2種實現工程的運算周期數對比見表4所列。

表4 2-D滑動FFT不同乘法器類型運算周期數

在參數滑動窗長為256、滑動步長為16的案例下,2-D 滑動FFT處理器與1-D滑動FFT處理器的單次滑動后,計算滑動窗長所需周期見表5所列。

表5 不同滑動FFT處理器單次滑動后滑動窗計算所需周期

4 結 論

本文給出了一種高速高精度的2-D滑動FFT處理器設計,采用IEEE-754單精度浮點數數據格式,支持256點步長16點的FFT計算,并完成了TSMC 28 nm工藝下后端實現。實驗結果表明,本文所設計的FFT處理器的工作頻率大于等于600 MHz,計算精度達到10-4,SNR大于130 dB。