對稀疏傅里葉變換并行算法研究并在FPGA上設(shè)計(jì)實(shí)現(xiàn)

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

0引言

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

2010年以來，一種并行結(jié)構(gòu)的SFT算法受到了廣泛的關(guān)注［5］。并行SFT首先通過并行下采樣，采集計(jì)算所需的所有時(shí)域數(shù)據(jù)，然后再通過FFT，通過亞線性頻譜估計(jì)方法獲得信號的稀疏頻率及其幅值。由于該類方法以并行取代迭代獲得頻譜估計(jì)所需的所有信息，因此可以實(shí)時(shí)處理各種頻域稀疏信號，使得經(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)方案都存在主存儲區(qū)與GPU存儲區(qū)的連接交互問題，因此數(shù)據(jù)間的正常流動不能得到更好的促進(jìn)。

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

1SFT并行算法

SFT并行算法主要由下采樣、頻率估計(jì)、幅度估計(jì)三個(gè)部分組成。在下采樣過程中，將輸信號劃分為N個(gè)組，每個(gè)組的采樣因子分別為σ1,σ2,…,σN。利用中國余數(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為采樣樣本。之后，對X~進(jìn)行DFT變換，得到X^=DFT(X~)。按幅度從大到小對X^向量中的值進(jìn)行排序，計(jì)算hk：

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

其中K為指定的重構(gòu)信號的參數(shù)。得到hk之后則可通過求余運(yùn)算獲取余數(shù)信息r1,khk的位置modq1。通過并行查詢的方式搜索余數(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通過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)，信號在經(jīng)過下采樣之后，通過FFT運(yùn)算得到復(fù)數(shù)的輸出信號，為了對該復(fù)數(shù)信號進(jìn)行排序，將該復(fù)數(shù)信號取模后送入排序網(wǎng)絡(luò)，由于每個(gè)通道送入排序網(wǎng)絡(luò)的點(diǎn)數(shù)不同，排序網(wǎng)絡(luò)的結(jié)構(gòu)會稍有差異。在利用CRT估計(jì)信號的幅度和頻率之前，需要對信號進(jìn)行求余、求最優(yōu)解等運(yùn)算。其中，最優(yōu)解運(yùn)算的核心是排序網(wǎng)絡(luò)，利用排序網(wǎng)絡(luò)的思想求取輸入信號的最大值以及獲取排序后的信號在原輸入信號中的位置；CRT模塊由一些加法器和乘法器組成。