摘要：提出了一種基于最優(yōu)搜索的稀疏傅里葉變換(SFT)的并行實(shí)現(xiàn)設(shè)計(jì)。首先將輸入信號(hào)分為并行N組，分別進(jìn)行快速傅里葉變換(FFT)，實(shí)現(xiàn)信號(hào)頻率分量的取模處理，然后通過(guò)排序搜索獲得。經(jīng)驗(yàn)證，相較于FFTW，當(dāng)信號(hào)長(zhǎng)度大于524 288時(shí)，執(zhí)行時(shí)間會(huì)有更好的表現(xiàn)；相較于正交匹配算法及其他SFT的FPGA實(shí)現(xiàn)，其系統(tǒng)的復(fù)雜度降低了。
?

0引言

稀疏傅里葉變換（Sparse Fourier Transform, SFT）是一種新的算法框架，也是快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）在處理稀疏頻譜信號(hào)上的延伸。2003年AYDINER A A等人提出了針對(duì)頻域稀疏信號(hào)的傅里葉變換基本思想［1］。對(duì)于頻域稀疏信號(hào)來(lái)說(shuō)，其頻譜可以通過(guò)其多級(jí)子集頻譜獲得。之后，IWEN M A等人從壓縮感知得到啟發(fā)，將采樣和頻率估計(jì)整合到快速傅里葉變換并提出了經(jīng)典SFT框架［2］。之后SFT廣泛運(yùn)用于稀疏頻譜信號(hào)（諸如音頻信號(hào)、醫(yī)學(xué)圖像信號(hào)）的處理以及頻譜感知領(lǐng)域［3］。大量的SFT算法被提出，它們多利用經(jīng)典的頻率估計(jì)算法通過(guò)亞奈奎斯特采樣點(diǎn)的子集的傅立葉變換重構(gòu)稀疏頻點(diǎn)［4］。但由于經(jīng)典SFT的亞奈奎斯特率樣本是通過(guò)多次采樣獲得的，因此，經(jīng)典SFT 不可能代替 FFT來(lái)處理實(shí)時(shí)信號(hào)，比如雷達(dá)信號(hào)等。

2010年以來(lái)，一種并行結(jié)構(gòu)的SFT算法受到了廣泛的關(guān)注［5］。并行SFT首先通過(guò)并行下采樣，采集計(jì)算所需的所有時(shí)域數(shù)據(jù)，然后再通過(guò)FFT，通過(guò)亞線性頻譜估計(jì)方法獲得信號(hào)的稀疏頻率及其幅值。由于該類方法以并行取代迭代獲得頻譜估計(jì)所需的所有信息，因此可以實(shí)時(shí)處理各種頻域稀疏信號(hào)，使得經(jīng)典SFT得到了改善。基于此，參考文獻(xiàn)［6］~［8］探討了稀疏傅里葉變換在GPU以及多核CPU上的實(shí)現(xiàn)方式。這些研究顯示，基于GPU加速的實(shí)現(xiàn)方案運(yùn)行速度要顯著高于基于CPU的實(shí)現(xiàn)方案。然而，基于GPU的實(shí)現(xiàn)方案都存在主存儲(chǔ)區(qū)與GPU存儲(chǔ)區(qū)的連接交互問(wèn)題，因此數(shù)據(jù)間的正常流動(dòng)不能得到更好的促進(jìn)。

為解決GPU的數(shù)據(jù)并行處理的局限性，本文研究SFT的并行算法并在FPGA上對(duì)其進(jìn)行實(shí)現(xiàn)，同時(shí)應(yīng)用中國(guó)余數(shù)定理（CRT）的基本原理對(duì)信號(hào)進(jìn)行重構(gòu)。相較于傳統(tǒng)的SFT，本文的方法可以極大地降低系統(tǒng)的復(fù)雜度，減少了硬件的開(kāi)銷。本文，首先介紹SFT的并行框架，然后討論SFT的FPGA實(shí)現(xiàn)架構(gòu)，最后從仿真結(jié)果以及硬件實(shí)現(xiàn)兩方面對(duì)系統(tǒng)進(jìn)行評(píng)估。

1SFT并行算法

SFT并行算法主要由下采樣、頻率估計(jì)、幅度估計(jì)三個(gè)部分組成。在下采樣過(guò)程中，將輸信號(hào)劃分為N個(gè)組，每個(gè)組的采樣因子分別為σ1,σ2,…,σN。利用中國(guó)余數(shù)定理（Chinese Residue Theorem，CRT）進(jìn)行頻率以及幅度的估計(jì)，設(shè)定各組的采樣因子兩兩互質(zhì)。

并行SFT算法從L個(gè)采樣樣點(diǎn)中重構(gòu)2 K算個(gè)參數(shù)，其中包括K個(gè)時(shí)延參數(shù)以及K個(gè)幅度參數(shù)。重構(gòu)法的具體計(jì)算步驟如下：首先，設(shè)采樣通道數(shù)為N，每個(gè)通道的采樣點(diǎn)數(shù)分別為qi,i=1,2,…,N。計(jì)算X~=M*X，其中M為下采樣矩陣，X為采樣樣本。之后，對(duì)X~進(jìn)行DFT變換，得到X^=DFT(X~)。按幅度從大到小對(duì)X^向量中的值進(jìn)行排序，計(jì)算hk：

hk=kthMAX|X^j|,k∈［0，K］（1）

其中K為指定的重構(gòu)信號(hào)的參數(shù)。得到hk之后則可通過(guò)求余運(yùn)算獲取余數(shù)信息r1,khk的位置modq1。通過(guò)并行查詢的方式搜索余數(shù)的最優(yōu)解：

tminmint∈［0，qj］|hk-X^t*q1+r1,k|,j∈［2，N］（2）

rj,k=r1,k+tmin*q1modqj,j∈［2，N］（3）

利用CRT通過(guò)r1,b,…，rN，b重構(gòu)時(shí)延參數(shù)τk，幅度估計(jì)參數(shù)可由公式（4）和 （5）得出：

x=real(hk)|τk

y=imag(hk)|τk(4)

ak=|x+iy|（5）

2SFT主要部分的FPGA實(shí)現(xiàn)

本文考慮使用MATLAB Simulink工具構(gòu)建SFT采樣算法的FPGA實(shí)現(xiàn)架構(gòu)。圖1展示了當(dāng)采樣通道數(shù)N=3時(shí)的SFT并行結(jié)構(gòu)，其主要包括下采樣、頻率估計(jì)、幅度估計(jì)三個(gè)部分?！?/p>

如圖1所示，頻率估計(jì)與幅度估計(jì)共用部分相同的硬件結(jié)構(gòu)，信號(hào)在經(jīng)過(guò)下采樣之后，通過(guò)FFT運(yùn)算得到復(fù)數(shù)的輸出信號(hào)，為了對(duì)該復(fù)數(shù)信號(hào)進(jìn)行排序，將該復(fù)數(shù)信號(hào)取模后送入排序網(wǎng)絡(luò)，由于每個(gè)通道送入排序網(wǎng)絡(luò)的點(diǎn)數(shù)不同，排序網(wǎng)絡(luò)的結(jié)構(gòu)會(huì)稍有差異。在利用CRT估計(jì)信號(hào)的幅度和頻率之前，需要對(duì)信號(hào)進(jìn)行求余、求最優(yōu)解等運(yùn)算。其中，最優(yōu)解運(yùn)算的核心是排序網(wǎng)絡(luò)，利用排序網(wǎng)絡(luò)的思想求取輸入信號(hào)的最大值以及獲取排序后的信號(hào)在原輸入信號(hào)中的位置；CRT模塊由一些加法器和乘法器組成。

輸入信號(hào)經(jīng)過(guò)多路選擇器獲得下采樣信號(hào)，所以該部分主要研究下采樣信號(hào)的頻率估計(jì)以及幅度估計(jì)，頻率估計(jì)包括最優(yōu)解模塊以及CRT重構(gòu)模塊。另外，硬件構(gòu)成部門還包含了存儲(chǔ)和控制單元，各通道采樣因子數(shù)ql、參數(shù)t、排序位置信息等都在存儲(chǔ)單元中保存，控制單元產(chǎn)生地址值來(lái)執(zhí)行讀寫存儲(chǔ)器的操作，并輸出必要的控制信號(hào)來(lái)初始化運(yùn)算模塊。

在本設(shè)計(jì)中，設(shè)定信號(hào)長(zhǎng)度N=223，參數(shù)個(gè)數(shù)K，采樣通道數(shù)M=3，其中，各個(gè)通道的采樣點(diǎn)分別為q1，q2，q3；q1,q2,q3兩兩互質(zhì)且乘積大于信號(hào)長(zhǎng)度N，因此，通過(guò)中國(guó)余數(shù)定理可由q1+q2+q3個(gè)采樣點(diǎn)數(shù)獲取原始信號(hào)所有的信息，降低了幅度以及頻率估計(jì)時(shí)所需的采樣點(diǎn)數(shù)。下面介紹各個(gè)主要功能模塊的設(shè)計(jì)。

2.1頻率估計(jì)

2.1.1最優(yōu)解模塊實(shí)現(xiàn)架構(gòu)

頻率估計(jì)模塊的核心部分在于最優(yōu)解的獲得，最優(yōu)解實(shí)現(xiàn)架構(gòu)如圖2所示。圖中X^t*q1+r1,k可由采樣樣本經(jīng)過(guò)多路選擇器MUX獲得，根據(jù)t的取值不同，可以得到多個(gè)定點(diǎn)采樣的樣點(diǎn)值。當(dāng)排序網(wǎng)絡(luò)為三輸入結(jié)構(gòu)時(shí)，r1,k代表排序網(wǎng)圖2最優(yōu)解模塊架構(gòu)絡(luò)模塊中獲得的k1,1,k1,2,k1,3對(duì)應(yīng)位置信息數(shù)據(jù)對(duì)q1取余所得的值。t的取值范圍為0~qj,j=2,3，其中當(dāng)j=2時(shí)，q2=4，同樣地，有q3=5；兩種情況下只有t值的變化導(dǎo)致定點(diǎn)采樣的樣點(diǎn)值的變化，而其對(duì)應(yīng)的模塊架構(gòu)是相同的，所以這里僅對(duì)j=2時(shí)的情況加以分析。

根據(jù)排序網(wǎng)絡(luò)結(jié)構(gòu)，需要的輸入數(shù)據(jù)有兩組，一組為需要排序的數(shù)據(jù)，以便求得最小值，另外一組則為數(shù)據(jù)對(duì)應(yīng)的位置信息t。這樣在排序網(wǎng)絡(luò)求取完最小值后可以直接獲取相應(yīng)的t值而不需要進(jìn)行其他的運(yùn)算處理。為此，將需要排序的數(shù)據(jù)并行導(dǎo)入排序網(wǎng)絡(luò)的數(shù)據(jù)輸入接口，將對(duì)應(yīng)的位置信息t值也并行導(dǎo)入排序網(wǎng)絡(luò)的位置信息接口。

如圖3所示，原有輸入的3路信號(hào)序號(hào)為1,2,3。該模塊實(shí)現(xiàn)對(duì)這3路信號(hào)進(jìn)行從大到小的排序，并獲得排序后的信號(hào)在原序列中的序號(hào)，即取位。圖3顯示了3輸入結(jié)構(gòu)的排序圖，4輸入乃至更多輸入結(jié)構(gòu)圖原理相同，圖中比較器的輸出作為多路選擇器的sel選擇端輸入，利用比較器以及多路選擇器的硬件電路連接實(shí)現(xiàn)邏輯上的比較選擇排序。k1,1,k1,2,k1,3為3輸入信號(hào)經(jīng)過(guò)排序網(wǎng)絡(luò)的輸出信號(hào)，有k1,1>k1,2>k1,3。k1,1_loc,k1,2_loc,k1,3_loc分別記錄了k1,1,k1,2,k1,3在原序列中的位置。同時(shí)將位置信息存儲(chǔ)到位置信息存儲(chǔ)器中。

2.1.2CRT模塊架構(gòu)

最優(yōu)解模塊輸出一組余數(shù)信息的集合，利用中國(guó)余數(shù)定理可以輕易地通過(guò)一組累加求和運(yùn)算獲取頻率集合，進(jìn)一步便可獲取時(shí)延參數(shù)τk。由中國(guó)余數(shù)定理可以得到如下方程組：

x≡r1(modq1)

x≡r2(modq2)



x≡rn(modqn) （6）

其中ri（i=1,2,…,n）為頻率點(diǎn)的集合，qi（i=1,2,…,n）為采樣點(diǎn)數(shù)的集合。假設(shè)Q為q1到qn的乘積，并設(shè)Qi=Q/qi,i∈{1,2,…,n}，ti為Qi模qi的數(shù)論倒數(shù)，則有：

x=∑ni=1ritiQi（7）

圖4顯示了一個(gè)頻率點(diǎn)的CRT重構(gòu)模塊架構(gòu)。

2.2幅度估計(jì)

幅度估計(jì)中，利用CRT重構(gòu)模塊中獲取的頻率集合w1,w2分別與L1,L2作求余運(yùn)算，以此為基礎(chǔ)求得hk,利用前面式（4）和式（5）可求得原始信號(hào)的幅度估計(jì)。其中頻率集合w1,w2由CRT模塊獲得，圖5中求余的作用為頻率集合w1,w2分別對(duì)采樣點(diǎn)數(shù)L1,L2作求余運(yùn)算。輸入序列xl、稀疏度值、采樣通道數(shù)、每個(gè)通道的采樣點(diǎn)數(shù)存儲(chǔ)在寄存器中供乘法器調(diào)用。利用排序網(wǎng)絡(luò)分別求得輸入信號(hào)實(shí)部與虛部的最大值，再對(duì)其進(jìn)行取模則可得到幅度值的估計(jì)。幅度估計(jì)的模型如圖5所示，其中，排序網(wǎng)絡(luò)為4輸入結(jié)構(gòu)。

3結(jié)果分析以及性能評(píng)估

為評(píng)估該算法框架的有效性，將其與FFTW做對(duì)比，F(xiàn)FTW是一個(gè)快速計(jì)算離散傅里葉變換的庫(kù),這個(gè)庫(kù)可以在多核CPU以及GPU上運(yùn)行。分別考慮稀疏度k恒定為1 000時(shí)信號(hào)長(zhǎng)度的變化對(duì)執(zhí)行時(shí)間的影響，以及信號(hào)長(zhǎng)度N恒定為223時(shí)稀疏度的變化對(duì)執(zhí)行時(shí)間的影響，比較結(jié)果如圖6所示。

將本文討論的稀疏傅里葉變換采樣框架與已知的OMP算法框架作性能上的對(duì)比，實(shí)現(xiàn)了信號(hào)長(zhǎng)度N=32，參數(shù)個(gè)數(shù)K=2以及采樣點(diǎn)數(shù)的采樣框架。其中，使用RAM塊實(shí)現(xiàn)所有所需的向量、常數(shù)或矩陣的存儲(chǔ)。將OMP架構(gòu)［9］以及SFT架構(gòu)［10］在同樣的平臺(tái)下做了實(shí)現(xiàn)來(lái)與本文算法架構(gòu)進(jìn)行對(duì)比，其結(jié)果如表1所示。

相較于OMP架構(gòu)，本文提出架構(gòu)大大減少了DSP48E以及所需寄存器的數(shù)量。相較于文獻(xiàn)［10］提出的SFT架構(gòu)，本文架構(gòu)依舊能夠有良好的表現(xiàn)。

4結(jié)論

本文提出了SFT的FPGA并行實(shí)現(xiàn)方案，使用Simulink中的XSG開(kāi)發(fā)工具構(gòu)建FPGA實(shí)現(xiàn)框架。對(duì)獨(dú)立功能塊的并行化處理可以大大減少執(zhí)行時(shí)間。之后對(duì)FPGA上的硬件實(shí)現(xiàn)進(jìn)行了評(píng)估，相對(duì)于FFTW的實(shí)現(xiàn)方案，在采樣點(diǎn)數(shù)的量級(jí)足夠大時(shí)，提高了系統(tǒng)運(yùn)行速度，降低了計(jì)算所需的時(shí)間；相對(duì)于其他OMP等算法的FPGA實(shí)現(xiàn)方案，減少了資源的消耗，降低了系統(tǒng)的復(fù)雜度。