面向FPGA的稀疏傅里葉并行算法實(shí)現(xiàn)

楊 超，錢 慧

(福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350108)

面向FPGA的稀疏傅里葉并行算法實(shí)現(xiàn)

楊 超，錢 慧

(福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350108)

提出了一種基于最優(yōu)搜索的稀疏傅里葉變換(SFT)的并行實(shí)現(xiàn)設(shè)計(jì)。首先將輸入信號(hào)分為并行N組，分別進(jìn)行快速傅里葉變換(FFT)，實(shí)現(xiàn)信號(hào)頻率分量的取模處理，然后通過排序搜索獲得。經(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ù)雜度降低了。

稀疏傅里葉變換；并行框架；現(xiàn)場(chǎng)可編程門陣列

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)來說，其頻譜可以通過其多級(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ì)算法通過亞奈奎斯特采樣點(diǎn)的子集的傅立葉變換重構(gòu)稀疏頻點(diǎn)[4]。但由于經(jīng)典SFT的亞奈奎斯特率樣本是通過多次采樣獲得的，因此，經(jīng)典SFT 不可能代替 FFT來處理實(shí)時(shí)信號(hào)，比如雷達(dá)信號(hào)等。

2010年以來，一種并行結(jié)構(gòu)的SFT算法受到了廣泛的關(guān)注[5]。并行SFT首先通過并行下采樣，采集計(jì)算所需的所有時(shí)域數(shù)據(jù)，然后再通過FFT，通過亞線性頻譜估計(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ū)的連接交互問題，因此數(shù)據(jù)間的正常流動(dòng)不能得到更好的促進(jìn)。

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

1 SFT并行算法

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

(1)

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

(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)得出：

(4)

ak=|x+iy|

(5)

2 SFT主要部分的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è)部分。

圖1 SFT的并行架構(gòu)

如圖1所示，頻率估計(jì)與幅度估計(jì)共用部分相同的硬件結(jié)構(gòu)，信號(hào)在經(jīng)過下采樣之后，通過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)過多路選擇器獲得下采樣信號(hào)，所以該部分主要研究下采樣信號(hào)的頻率估計(jì)以及幅度估計(jì)，頻率估計(jì)包括最優(yōu)解模塊以及CRT重構(gòu)模塊。另外，硬件構(gòu)成部門還包含了存儲(chǔ)和控制單元，各通道采樣因子數(shù)ql、參數(shù)t、排序位置信息等都在存儲(chǔ)單元中保存，控制單元產(chǎn)生地址值來執(zhí)行讀寫存儲(chǔ)器的操作，并輸出必要的控制信號(hào)來初始化運(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，因此，通過中國余數(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)

圖2 最優(yōu)解模塊架構(gòu)

根據(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)過排序網(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ǔ)器中。

圖3 排序網(wǎng)絡(luò)結(jié)構(gòu)(3輸入結(jié)構(gòu))

2.1.2 CRT模塊架構(gòu)

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

(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ù)，則有：

(7)

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

圖4 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)。

圖5 幅度估計(jì)計(jì)算框架

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

為評(píng)估該算法框架的有效性，將其與FFTW做對(duì)比，F(xiàn)FTW是一個(gè)快速計(jì)算離散傅里葉變換的庫,這個(gè)庫可以在多核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)來與本文算法架構(gòu)進(jìn)行對(duì)比，其結(jié)果如表1所示。

表1 本文提出架構(gòu)與已有架構(gòu)實(shí)施效果以及FPGA資源利用率對(duì)比 (%)

相較于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開發(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ù)雜度。

[1] AYDINER A A,WENG C C,SONG J,et al.A sparse data fast Fourier transform (SDFFT)[J].IEEE Transactions on Antennas & Propagation,2003,51(11):3161-3170.

[2] IWEN M A.A deterministic sub-linear time sparse Fourier algorithm via non-adaptive compressed sensing methods[C].Proceedings of the nineteenth annual ACM-SIAM symposium on Discrete Algorithms,2008:20-29.

[3] 那美麗,周志剛,李霈霈.基于稀疏傅里葉變換的低采樣率寬帶頻譜感知[J].電子技術(shù)應(yīng)用,2015,41(11):85-88.

[4] GILBERT A C,STRAUSS M J,TROPP J A.A tutorial on fast fourier sampling[J].IEEE Signal Processing Magazine,2008,25(2):57-66.

[5] Wang Cheng,ARAYA-POLO M,CHANDRASEKARAN S,et al.Parallel sparse FFT[C].Proceedings of the 3rd Workshop on Irregular Applications:Architectures and Algorithms.ACM,2013:10:1-10:8.

[6] Hu Jiaxi,Wang Zhaosen,Qiu Qiyuan,et al.Sparse fast Fourier transform on GPUs and multi-core CPUs[C].2012 IEEE 24th International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD),IEEE,2012:83-91.

[7] BRAUN T R.An evaluation of GPU acceleration for sparse reconstruction[J].Proceedings of SPIE-The International Society for Optical Engineering,2010,7697.

[8] Wang Cheng,CHANDRASEKARAN S,CHAPMAN B.cusFFT:a high-performance sparse fast fourier transform algorithm on GPUs[C].Proceedings of 30th IEEE International Parallel & Distributed Processing Symposium (IPDPS),2016:936-972.

[9] RABAH H,AMIRA A,MOHANTY B K,et al.FPGA implementation of orthogonal matching pursuit for compressive sensing reconstruction[J].IEEE Transactions on Very Large Scale Integration Systems,2015,23(10):2209-2220.

[10] AGARWAL A,HASSANIEH H,ABARI O,et al.High-throughput implementation of a million-point sparse Fourier Transform[C].IEEE Conference on Field Programmable Logic and Applications (FPL),2014:1-6.

The parallel framework of sparse Fourier transform on FPGA

Yang Chao,Qian Hui

(College of Physics and Information Engineering,Fuzhou University,Fuzhou 350108,China)

This paper proposed the parallel framework of sparse Fourier transform and implemented all the design on FPGA .The input signal is divided intoNgroups to search for the optimal solution.It is proved that the design compared to the fast Fourier transform in the West(FFTW),the execution time has a better performance when the signal length is greater than a specific value 524 288,and compared to the FPGA implementation of orthogonal matching pursuit algorithm and other implementation of sparse Fourier transform algorithm,this design reduces the complexity of the system.

sparse Fourier transform；parallel framework；FPGA

TN911.7

A

10.19358/j.issn.1674- 7720.2017.10.020

楊超，錢慧.面向FPGA的稀疏傅里葉并行算法實(shí)現(xiàn)[J].微型機(jī)與應(yīng)用，2017,36(10)：70-73.

2016-11-21)

楊超(1992-)，男，碩士研究生，主要研究方向：信號(hào)處理、圖像處理。

錢慧(1977-)，通信作者，女，博士，講師，主要研究方向：通信信號(hào)處理、壓縮采樣等。E-mail:qianhui@fzu.edu.cn。