基于改進的快速Fermat數變換的卷積算法及其FPGA實現

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

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

引 言

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

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

1 快速費爾馬數變換(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的最大公因數。當M為質數時，上面的所有條件均自動滿足，模M=2b+1的NTT變換就稱為費爾馬數變換，特別是當b=2t時FNT存在快速算法。觀察NTT變換對可以看出，不同于DFT變換，FNT變換為實變換，但FNT變換具有與DFT變換相似的周期性、循環卷積等特性。

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

蝶形運算用D1編碼系統實現只需要加法器和移位寄存器［8］，基于時域抽取的FFNT輸入數據為逆序，輸出為順序，每個蝶形運算需要2個加法器、2個存儲單元和1個運算周期。基于變換域抽取的FFNT輸入數據為順序，輸出為逆序，每個蝶形運算需要2個加法器、3個存儲單元和2個運算周期。可見，FFNT算法不需要乘法器，其運算時間大大縮減。

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

2.1 FFNT、PFNT、多維 FNT算法

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

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

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

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

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

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

2.2 MFFNT算法

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

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

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

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

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

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

2.3 ERNS復乘

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

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

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

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

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

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

3 基于MFFNT的卷積算法

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

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

MFFNT算法繼承了FNT算法的優點，但是與FNT不同，MFFNT是ERNS域上的復數變換，因此基于MFFNT的ERNS卷積與圓周卷積是有區別的。

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

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

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

圖4是在Virtex6-LX240T平臺上實現卷積時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平臺上實現卷積時MFFNT算法和FFT算法的資源消耗對比，圖中左邊是基于MFFNT的卷積算法的資源消耗，右邊是基于FFT的卷積算法的資源消耗。從圖中可以看出，實現長復數序列卷積時，基于MFFNT的算法會消耗更多的邏輯資源LUT，但是比基于FFT的算法節約了約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 結束語

本文從FFNT的周期性、對稱性出發，對FFNT的蝶形運算結構進行分析，通過研究FFNT、PFNT及多維映射類FNT的優缺點，提出了一種復數域的MFFNT算法，并利用MATALB、ISE驗證了Viretex6 FPGA平臺上采用MFFNT算法實現線性卷積的可行性。仿真及應用結果表明，基于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］陶 濤，初建朋，賴宗聲，等.一種可參數化快速FNT的FPGA實現［J］.微電子學與計算機，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的數論變換算法及應用的研究［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變換完成大點數循環卷積IP核的VLSI實現［J］.微電子學與計算機，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 著.數論概述［M］.孫智偉，等，譯.北京:機械工業出版社，2008，5.

［8］U.Meyer-Baese著.數字信號處理的FPGA實現 ［M］.劉 凌，譯.北京:清華大學出版社，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.