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

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

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

針對此類問題,文獻[6]使用滑動FFT算法?；瑒覨FT[7]的執(zhí)行過程為迭代累加累乘計算,可在前一次滑動FFT結(jié)果基礎(chǔ)上遞推得到當前滑動FFT結(jié)果,復(fù)用前次運算結(jié)果,節(jié)約運算時間,提高信號分析實時性。利用滑動FFT運算遞推特征,可以根據(jù)信號分析需要選擇計算相應(yīng)的頻譜區(qū)段,增加信號分析效果。

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

1 算法原理

1.1 滑動FFT原理

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

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

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點數(shù)據(jù)的復(fù)制,并不需要用到加法和乘法運算。N=16、L=2的Pruning算法圖如圖1所示。圖1中:虛線表示0值;實線表示輸入序列值。

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

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

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

(5)

滑動FFT算法需要進行NL/2+N次復(fù)數(shù)乘法、NL+N次復(fù)數(shù)加法,而傳統(tǒng)的FFT算法需要進行(NlbN)/2次復(fù)數(shù)乘法、NlbN次復(fù)數(shù)加法?；瑒哟伴L度為64～256點,滑動步長為4～16點,滑動窗每滑動一次產(chǎn)生的復(fù)數(shù)乘法計算量與復(fù)數(shù)加法計算量的關(guān)系見表1所列。由表1可知,滑動FFT所需復(fù)數(shù)乘法數(shù)量和復(fù)數(shù)加法數(shù)量都小于傳統(tǒng)FFT算法。

表1 滑動窗長、滑動步長與運算量的關(guān)系

1.2 2-D FFT算法原理

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

對N點數(shù)據(jù)進行FFT運算,令N=LC,即將輸入的N點分解為L行、C列的矩陣形式。根據(jù)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)的序列順序是滿足先填寫行方向再填寫列方向,和輸入數(shù)據(jù)相逆,所以還需要一次整序才是正確的FFT序列。

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

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 硬件架構(gòu)選擇

本文采用2-D滑動FFT硬件結(jié)構(gòu),從以下方面進行考慮。

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

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

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

為達到計算速度要求,本文設(shè)計行方向采用全展開FFT結(jié)構(gòu),對±1和±i的旋轉(zhuǎn)因子進行操作,采用交換虛實部的方法壓縮運算,對其他旋轉(zhuǎn)因子系數(shù)使用優(yōu)化編碼方法壓縮運算。

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

其次是壓縮浮點數(shù)乘法資源。滑動FFT運算中的矯正運算是影響設(shè)計速度的關(guān)鍵,以滑動窗長為256點、滑動步長為16為例,滑動FFT運算的矯正因子系數(shù)為±0.707 1、±0.382 7、±0.923 9。對上述因子使用高基Booth編碼操作將浮點數(shù)乘法轉(zhuǎn)變?yōu)槌?shù)乘法,優(yōu)化邏輯結(jié)構(gòu),減少資源。

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

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

2 系統(tǒng)結(jié)構(gòu)

本文設(shè)計的2-D FFT處理器結(jié)構(gòu)如圖3所示。

圖3 2-D FFT處理器結(jié)構(gòu)

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

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

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

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

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

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

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

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

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

2.1 基16蝶形運算器設(shè)計

設(shè)計的滑動FFT矯正模塊需要在每個周期內(nèi)完成16點的復(fù)數(shù)加法與復(fù)數(shù)乘法,相應(yīng)地,本文設(shè)計必須在每個周期輸出16點的FFT結(jié)果,因此采用基16核的全展開結(jié)構(gòu)。

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

圖5 16點蝶形單元結(jié)構(gòu)

2.2 使用高基Booth編碼優(yōu)化定點數(shù)乘法

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

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

本文設(shè)計使用高基Booth編碼進行系數(shù)優(yōu)化[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等),因為每個系數(shù)都可找到一個基獲得最小加法執(zhí)行次數(shù),所以本文設(shè)計中0.707 1、0.382 7使用基16編碼,0.923 9使用基32編碼,即上述3個系數(shù)的乘法操作分別轉(zhuǎn)換為8次加法操作,達到了資源最小化目標。

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

(a) 0.382 7

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

2.3 滑動FFT的矯正模塊旋轉(zhuǎn)因子結(jié)構(gòu)設(shè)計

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

每16個值一次循環(huán),對應(yīng)著2-D FFT運算結(jié)果,矯正因子的排列方式如圖7所示。

圖7 滑動FFT矯正因子排列

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

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

3 設(shè)計實現(xiàn)與性能分析

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

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

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

圖8 滑動窗FFT版圖

參考結(jié)果為MATLAB中FFT函數(shù)的輸出。

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的不同頻率的序列數(shù)據(jù)。滑動窗長度256,滑動步長16,一共產(chǎn)生256段數(shù)據(jù)。在滑動了不同長度之后進行數(shù)據(jù)復(fù)位,保證誤差精度后的誤差分析,結(jié)果見表3所列。

表3 滑動FFT復(fù)位前10段結(jié)果數(shù)據(jù)的誤差分析

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

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

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

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

4 結(jié) 論

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