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

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

引 言

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

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

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

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

圖1 FNT時(shí)域抽取和變換域抽取算法的蝶形單元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ù)時(shí)，上面的所有條件均自動滿足，模M=2b+1的NTT變換就稱為費(fèi)爾馬數(shù)變換，特別是當(dāng)b=2t時(shí)FNT存在快速算法。觀察NTT變換對可以看出，不同于DFT變換，F(xiàn)NT變換為實(shí)變換，但FNT變換具有與DFT變換相似的周期性、循環(huán)卷積等特性。

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

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

2 改進(jìn)的快速費(fèi)爾馬數(shù)變換(MFFNT)算法

2.1 FFNT、PFNT、多維 FNT算法

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

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

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

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

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

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

2.2 MFFNT算法

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

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

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

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

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

③4點(diǎn)FFNT的計(jì)算分兩個(gè)階段，每個(gè)階段按時(shí)域抽取算法進(jìn)行兩個(gè)蝶形單元的運(yùn)算，每個(gè)蝶形單元的運(yùn)算時(shí)間為一個(gè)時(shí)鐘周期，資源消耗為4個(gè)加法器和4個(gè)存儲單元。

2.3 ERNS復(fù)乘

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

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

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

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

因此，3點(diǎn)MFFNT中的乘法同32點(diǎn)MFFNT中的蝶形單元相比運(yùn)算時(shí)間更短，資源消耗更少，只需要兩個(gè)寄存器和一個(gè)加法器。

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

3 基于MFFNT的卷積算法

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

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

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

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

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

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

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

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

圖5是在Virtex6-LX240T平臺上實(shí)現(xiàn)卷積時(shí)MFFNT算法和FFT算法的資源消耗對比，圖中左邊是基于MFFNT的卷積算法的資源消耗，右邊是基于FFT的卷積算法的資源消耗。從圖中可以看出，實(shí)現(xiàn)長復(fù)數(shù)序列卷積時(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的蝶形運(yùn)算結(jié)構(gòu)進(jìn)行分析，通過研究FFNT、PFNT及多維映射類FNT的優(yōu)缺點(diǎn)，提出了一種復(fù)數(shù)域的MFFNT算法，并利用MATALB、ISE驗(yàn)證了Viretex6 FPGA平臺上采用MFFNT算法實(shí)現(xiàn)線性卷積的可行性。仿真及應(yīng)用結(jié)果表明，基于MFFNT的卷積算法相比基于FFT的卷積算法需要更少的乘法器資源，運(yùn)算速度更快。

［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實(shí)現(xiàn)［J］.微電子學(xué)與計(jì)算機(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的數(shù)論變換算法及應(yīng)用的研究［J］.微計(jì)算機(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變換完成大點(diǎn)數(shù)循環(huán)卷積IP核的VLSI實(shí)現(xiàn)［J］.微電子學(xué)與計(jì)算機(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 著.數(shù)論概述［M］.孫智偉，等，譯.北京:機(jī)械工業(yè)出版社，2008，5.

［8］U.Meyer-Baese著.數(shù)字信號處理的FPGA實(shí)現(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.