基于改進的快速Fermat數(shù)變換的卷積算法及其FPGA實現(xiàn)

周勝文， 周云生， 詹 磊， 董 暉

(北京遙測技術(shù)研究所北京100076)

引 言

卷積的快速實現(xiàn)是信號和圖像處理中經(jīng)常遇到的問題，在實踐中，這些操作通常都采用快速傅里葉變換(FFT)算法實現(xiàn)，但在某些場合數(shù)論變換(NTT)要優(yōu)于FFT。1971年P(guān)ollard［1］在有限群上定義了NTT，并給出了NTT變換對存在的條件。Mcclellan［2］于1976年提出了一種針對模運算的編碼方式，從硬件角度分析了費爾馬數(shù)論變換(FNT)的實現(xiàn)結(jié)構(gòu)和流程，構(gòu)建了一個64點16bit的FNT原型，并在此基礎(chǔ)上探討了長數(shù)據(jù)FNT的實現(xiàn)。Leibowitz［3］改進了Mcclellan的模運算編碼方式，提出專門用于FNT的D1編碼。從上世紀(jì)80年代起，我國學(xué)者也針對數(shù)論變換的應(yīng)用展開研究。文獻［4］分析了基4快速數(shù)論變換的FPGA實現(xiàn)，利用Xilinx Virtex2芯片實現(xiàn)了64點費爾馬變換，但算法的運算時間太長。文獻［5］討論了基于蝶形運算的快速費爾馬數(shù)論變換的FPGA實現(xiàn)，但是并未涉及FNT變換長度和位寬等關(guān)鍵問題。文獻［6］研究了三個N點序列FNT合成一個2N點序列FNT的IP核實現(xiàn)方法，但該方法未能解決FNT變換長度增加與資源消耗急劇增大之間的矛盾。一直以來，人們不斷研究通過類似Good-Thomas映射的多維分解技術(shù)解決FNT變換長度增加的問題，但是作大點數(shù)FNT時需要剔除模運算中的壞因子，偽費爾馬變換(PFNT)能剔除壞因子卻不利于模運算的FPGA實現(xiàn)。

本文從分析快速費爾馬數(shù)論變換(FFNT)的基本蝶形算法出發(fā)，討論了FFNT算法的應(yīng)用局限及改進方法，之后提出了改進的FFNT(MFFNT)算法，將FFNT擴展到復(fù)數(shù)域來實現(xiàn)長復(fù)數(shù)序列的線性卷積。該MFFNT算法能與FPGA平臺緊密結(jié)合，易于實際工程應(yīng)用。

1 快速費爾馬數(shù)變換(FFNT)

1971年，Pollard在有限群上定義了NTT變換對

圖1 FNT時域抽取和變換域抽取算法的蝶形單元Fig.1 Butterfly unit of FNT algorithm with the decimation in time domain and transform domain

其中，gcd(α，M)是α和M的最大公因數(shù)。當(dāng)M為質(zhì)數(shù)時，上面的所有條件均自動滿足，模M=2b+1的NTT變換就稱為費爾馬數(shù)變換，特別是當(dāng)b=2t時FNT存在快速算法。觀察NTT變換對可以看出，不同于DFT變換，F(xiàn)NT變換為實變換，但FNT變換具有與DFT變換相似的周期性、循環(huán)卷積等特性。

因此，F(xiàn)NT算法可以通過基2時域抽取和變換域抽取的FFNT算法實現(xiàn)，如圖1所示。與FFT算法相似，F(xiàn)FNT算法的基本單元是蝶形運算單元。

蝶形運算用D1編碼系統(tǒng)實現(xiàn)只需要加法器和移位寄存器［8］，基于時域抽取的FFNT輸入數(shù)據(jù)為逆序，輸出為順序，每個蝶形運算需要2個加法器、2個存儲單元和1個運算周期?；谧儞Q域抽取的FFNT輸入數(shù)據(jù)為順序，輸出為逆序，每個蝶形運算需要2個加法器、3個存儲單元和2個運算周期?？梢?，F(xiàn)FNT算法不需要乘法器，其運算時間大大縮減。

2 改進的快速費爾馬數(shù)變換(MFFNT)算法

2.1 FFNT、PFNT、多維 FNT算法

如第1節(jié)所述，F(xiàn)FNT算法在運算時間、資源消耗上具有很多優(yōu)點，但是實現(xiàn)長度N=2b=2t+1的FFNT算法的前提是模運算能夠在硬件上快速實現(xiàn)，而滿足如此特性的參數(shù)N和M并不多，更重要的是D1系統(tǒng)實現(xiàn)模運算的位寬W=b+1=2t+1，當(dāng)變換長度增加時位寬會成比例增大，此時FNT算法在FPGA平臺上實現(xiàn)時所需要的資源大大增加。

M非質(zhì)數(shù)、變換長度N不是費爾馬數(shù)的NTT稱作偽變換［8］。對于偽費爾馬數(shù)變換(PFNT)，首先根據(jù)變換長度N選取M，然后剔除M中的壞因子得到修正的M′，盡管PFNT算法的M′相比M小很多，但D1系統(tǒng)模運算的位寬仍很大，而且化簡后的模M′≠2+1，難以在FPGA平臺上實現(xiàn)。

因此，采用多維索引映射［9］，將長序列分解為二維或多維短序列，是增加FNT變換長度的一個有效方案。常用的多維映射變換有Cooley-Tukey變換、Good-Thomas變換及Agarwal-Burrus變換等，其基本原理是將總長度N=N1N2的變換分解為長度分別為N1和N2的兩個短序列的變換。輸入域的索引變換和輸出域的索引變換分別如下

Cooley和Tukey提出的索引變換是最簡單的映射，即令Cooley-Tukey變換的分解因子N1、N2可以是不互質(zhì)的。Cooley-Tukey索引變換算法首先對輸入序列進行索引變換，然后計算N2個N1點變換，變換結(jié)果乘以旋轉(zhuǎn)因子后進行N1個N2點變換，最后對輸出序列進行索引變換。

Good和Thomas提出的索引變換的最大特點是不需要乘以旋轉(zhuǎn)因子，其代價是分解因子N1、N2必須是互質(zhì)的，其中Good-Thomas索引變換算法與Cooley-Tukey索引變換算法類似，只是兩次短序列變換運算之間不需要乘以旋轉(zhuǎn)因子。

綜上，采用多維索引映射雖然可以大大增加變換長度，但是必須要找到易于FPGA平臺實現(xiàn)的FNT算法參數(shù)，由此可見將實數(shù)FFNT變換擴展到復(fù)數(shù)域是必然之舉。

2.2 MFFNT算法

Dimitrov.V［10］在研究基3-FFT算法時提出了Eisenstein Residue Number System(ERNS)。ERNS是在閉集Z［μ］上定義的加法和乘法運算，其中

通過選取不同的FNT復(fù)數(shù)變換基α，可以大大增加定義在艾森斯坦余數(shù)系統(tǒng)上的FNT變換的有效長度。例如，當(dāng)α=-2μ、M=2b+1時，式(1)拓展為ERNS域上的復(fù)數(shù)變換，變換長度擴展為N=6b，相比FFNT這一長度增加了3倍。因此，MFFNT算法將FFNT多維映射算法拓展到ERNS域，變換長度大大增加而實現(xiàn)變換所需要的位寬卻不變，這非常有利于算法在FPGA平臺上實現(xiàn)。另外，IMFFNT為MFFNT的逆變換，其變換過程與MFFNT算法相似，變換基為α的逆元即β=α-1modM，此處不再贅述。

MFFNT是在ERNS域的復(fù)變換，下面以變換基α=-2μ的96點復(fù)序列為例來闡述MFFNT算法的實現(xiàn)流程，具體的FPGA實現(xiàn)流程如圖2所示。

①采用Good-Thomas變換將96點MFFNT分解為N1=3和N2=32的二維MFFNT，3個32點變換的變換基為-8，32個3點變換的變換基32點MFFNT變換為實變換，3點MFFNT變換為復(fù)變換。

②32點MFFNT變換按照Cooley-Tukey映射分解為L1=8、L2=4的二維變換，即4個8點MFFNT變換后乘以旋轉(zhuǎn)因子再作8個4點MFFNT變換，乘以旋轉(zhuǎn)因子矩陣通過移位寄存器實現(xiàn)，8點MFFNT變換可類似FFT由奇4點FFNT和偶4點FFNT合成。

③4點FFNT的計算分兩個階段，每個階段按時域抽取算法進行兩個蝶形單元的運算，每個蝶形單元的運算時間為一個時鐘周期，資源消耗為4個加法器和4個存儲單元。

2.3 ERNS復(fù)乘

MFFNT與FFNT的根本區(qū)別在于ERNS域復(fù)數(shù)的實部和虛部相比普通復(fù)數(shù)存在交織，ERNS域的復(fù)數(shù)乘法滿足式(9)，如何在FPGA平臺上用最少的乘法器資源實現(xiàn)ERNS域復(fù)數(shù)乘法是優(yōu)化MFFNT卷積算法資源消耗的重要一環(huán)。如果采用實、虛部分別計算的方法，則需要4個加法器和4個乘法器。

Virtex6平臺中的專用乘法器DSP48E1提供了強大的運算輔助功能，利用DSP48E1的預(yù)加器及級聯(lián)端口，將式(9)改寫為式(10)可知，任意兩個復(fù)數(shù)的ERNS域乘法只需要3個DSP48E1就可以實現(xiàn)。

圖2 96點MFFNT的FPGA實現(xiàn)流程Fig.2 Flow chart of ninety-six point MFFNT implemented on FPGA

而對于3點MFFNT中的復(fù)數(shù)乘法，由于變換基αpoint3組成的變換矩陣為

因此，3點MFFNT中的乘法同32點MFFNT中的蝶形單元相比運算時間更短，資源消耗更少，只需要兩個寄存器和一個加法器。

另外，DSP48E1專用乘法器實現(xiàn)的是二進制補碼的乘法運算，最終的運算結(jié)果還必須通過模化簡轉(zhuǎn)換成D1碼，以方便進行IMFFNT運算。兩個16位寬的數(shù)A、B相乘的結(jié)果可表示為高16位L和低16位S的組合，如式(12)所示，MFFNT的模運算在FPGA平臺上用一個加法器就可以實現(xiàn)。

3 基于MFFNT的卷積算法

若x和y是模M定義的長度為N的序列，則x和y的循環(huán)卷積可以利用NTT來實現(xiàn)，令X、Y分別是x、y的NTT，則有

其中，Θ為ERNS域卷積。從式(13)可知，序列x和y的循環(huán)卷積可以通過兩次NTT變換、一次N點復(fù)乘和一次NTT逆變換來實現(xiàn)，該特性與DFT類似，稱為NTT的循環(huán)卷積特性。利用NTT實現(xiàn)循環(huán)卷積相比FFT實現(xiàn)卷積，其優(yōu)勢在于不需要消耗任何乘法器資源，并且運算時間可以大大減少。

MFFNT算法繼承了FNT算法的優(yōu)點，但是與FNT不同，MFFNT是ERNS域上的復(fù)數(shù)變換，因此基于MFFNT的ERNS卷積與圓周卷積是有區(qū)別的。

其中，?為圓周卷積，?為ERNS域復(fù)乘。由式(14)可知，從ERNS域卷積變換到圓周卷積需要進行單位向量的變換，即對MFFNT的輸入進行ERNS序列變換，對IMFFNT的輸出進行IERNS序列變換，ERNS變換對如式(15)所示。

圖3是基于MFFNT的卷積算法仿真結(jié)果。首先由Matlab分別生成32點和352點復(fù)隨機序列，然后用ISE和Modelsim仿真軟件讀入相應(yīng)的文本文檔，之后將卷積算法的FPGA結(jié)果寫到輸出文件，最后通過Matlab軟件讀出卷積結(jié)果，并與理論值進行對比分析。從仿真對比圖可以看出，基于MFFNT的卷積算法正確可行。

圖3 基于MFFNT的卷積算法仿真結(jié)果Fig.3 Simulation results of the convolution algorithm based on MFFNT

圖4是在Virtex6-LX240T平臺上實現(xiàn)卷積時MFFNT算法和FFT算法的運算時間對比，圖中左邊基于MFFNT的卷積算法的運算時間為34個時鐘周期，右邊基于FFT的卷積算法的運算時間為40個時鐘周期?？梢?，基于MFFNT的卷積算法相比基于FFT的卷積算法運算時間更少。

圖4 基于MFFNT和FFT的卷積算法在Virtex6-LX240T平臺上的運算時間對比Fig.4 Comparison of calculation time between the convolution algorithms based on MFFNT and FFT on the Virtex6-LX240T platform

圖5是在Virtex6-LX240T平臺上實現(xiàn)卷積時MFFNT算法和FFT算法的資源消耗對比，圖中左邊是基于MFFNT的卷積算法的資源消耗，右邊是基于FFT的卷積算法的資源消耗。從圖中可以看出，實現(xiàn)長復(fù)數(shù)序列卷積時，基于MFFNT的算法會消耗更多的邏輯資源LUT，但是比基于FFT的算法節(jié)約了約13%的專用乘法器DSP48E1，因此功耗相對較低。

圖5 基于MFFNT和FFT的卷積算法在Virtex6-LX240T平臺上的資源消耗對比Fig.5 Comparison of resource consumption between the convolution algorithms based on MFFNT and FFT on the Virtex6-LX240T platform

4 結(jié)束語

本文從FFNT的周期性、對稱性出發(fā)，對FFNT的蝶形運算結(jié)構(gòu)進行分析，通過研究FFNT、PFNT及多維映射類FNT的優(yōu)缺點，提出了一種復(fù)數(shù)域的MFFNT算法，并利用MATALB、ISE驗證了Viretex6 FPGA平臺上采用MFFNT算法實現(xiàn)線性卷積的可行性。仿真及應(yīng)用結(jié)果表明，基于MFFNT的卷積算法相比基于FFT的卷積算法需要更少的乘法器資源，運算速度更快。

［1］Pollard J.The Fast Fourier Tranform in a Finite Field［J］.Mathmatics of Computation，1971，25:365 ～374.

［2］Mcclellan JH.Harware Realization of a Fermat Number Transform［J］.IEEE Trans.on ASSP，1976，24(3):216 ～226.

［3］Leibowitz LM.A Binary Arithmetic for the Fermat Number Transform［R］.Naval Research Laboratory Report，1976.

［4］陶 濤，初建朋，賴宗聲，等.一種可參數(shù)化快速FNT的FPGA實現(xiàn)［J］.微電子學(xué)與計算機，2004，21(10):165～168.Tao Tao，Chu Jianpeng，Lai Zongsheng，etc.A Parameterized FPGA Realization of High-Speed FNT Processor［J］.Microelectronics & Computer，2004，21(10):165 ～168.

［5］余漢成，王成華，邵 杰，夏永君.基于FPGA的數(shù)論變換算法及應(yīng)用的研究［J］.微計算機信息，2006，22(11-2):212～214.Yu Hancheng，Wang Chenghua，Shao Jie，Xia Yongjun.The Research into FPGA-based NTTAlgorithm and Application［J］.Microcomputer Information，2006，22(11-2):212 ～214.

［6］李新兵，初建朋，賴宗聲，等.一種用FNT變換完成大點數(shù)循環(huán)卷積IP核的VLSI實現(xiàn)［J］.微電子學(xué)與計算機，2004，21(11):158～160.Li Xinbing，Chu Jianpeng，Lai Zongsheng，etc.A Parameterized FPGA Realization of High-Speed FNT Processor［J］.Microelectronics & Computer， 2004，21(11):158 ～160.

［7］Joseph H.Silverman 著.數(shù)論概述［M］.孫智偉，等，譯.北京:機械工業(yè)出版社，2008，5.

［8］U.Meyer-Baese著.數(shù)字信號處理的FPGA實現(xiàn) ［M］.劉 凌，譯.北京:清華大學(xué)出版社，2012，4.

［9］Burrus C.Index Mappings for Multidimensional Formulation of the DFT and Convolution［J］.IEEE Transactions on Acoustic，Speech and Signal Processing， 1977，25，239 ～242.

［10］Dimitrov V，Jullien G A，MillerW C.Eisenstein Residue Number System with Applications to DSP［C］//Proceedings of the 40th Midwest Symposium on Circuits and Systems，1997，675 ～678.