一種稀疏度自適應(yīng)的稀疏傅里葉變換算法

劉 仲,李立春,李慧啟

(信息工程大學(xué) 信息系統(tǒng)工程學(xué)院,鄭州 450001)

0 概述

在信號(hào)處理領(lǐng)域,離散傅里葉變換(Discrete Fourier Transform,DFT)是最基本的算法之一。而快速傅里葉變換(Fast Fourier Transform,FFT)是最快的DFT算法之一,其算法復(fù)雜度只有O(nlogan)。目前,計(jì)算DFT的算法復(fù)雜度為Θ(n),均與信號(hào)長(zhǎng)度成正比。在圖像處理[1]、寬帶頻譜感知[2-4]、機(jī)器學(xué)習(xí)理論[5]、多尺度分析[6]、衛(wèi)星通信[7]等眾多科學(xué)領(lǐng)域中,許多實(shí)際信號(hào)都有一定程度的稀疏性或可壓縮性。這種結(jié)構(gòu)特殊、廣泛存在的信號(hào)成為了信號(hào)處理領(lǐng)域熱門(mén)的研究方向。針對(duì)離散傅里葉變換,幾種頻域稀疏的亞線(xiàn)性算法[8-11]相繼被提出,但這些算法由于結(jié)構(gòu)的限制只適用于非常稀疏的信號(hào)。

近幾年,稀疏快速傅里葉變換(sparse Fast Fourier Transform,sFFT)算法成為研究稀疏信號(hào)的一大熱點(diǎn),并相繼提出了一些實(shí)用的sFFT算法。該類(lèi)算法利用信號(hào)頻域的稀疏性對(duì)頻率進(jìn)行“分桶”,使DFT運(yùn)算的點(diǎn)數(shù)成倍減少,然后通過(guò)其他方法實(shí)現(xiàn)對(duì)原信號(hào)完整頻譜的重構(gòu),從而使算法的復(fù)雜度與信號(hào)長(zhǎng)度成亞線(xiàn)性關(guān)系,進(jìn)而提高信號(hào)處理速度。近年來(lái)稀疏傅里葉變換的算法碩果累累,研究人員提出了許多算法。這些算法降低了算法復(fù)雜度[12-13],拓寬了維度適用范圍[14],并實(shí)現(xiàn)了算法的并行加速[15],并且推廣到了其他變換域[16]。但稀疏傅里葉變換都需要信號(hào)以傅氏域的稀疏度為先驗(yàn)信息,通常稀疏度K是未知的,在一定程度上限制了上述算法的應(yīng)用。為此,針對(duì)未知稀疏度的信號(hào),本文提出一種稀疏度自適應(yīng)的稀疏傅里葉變換算法(Sparsity Adaptive Algorithm for Sparse Fast Fourier Transform,SAsFFT)。

1 稀疏傅里葉變換算法

稀疏傅里葉變換算法通過(guò)時(shí)域降維對(duì)頻域稀疏的信號(hào)進(jìn)行處理,將信號(hào)的頻域信息壓縮到低維空間中,在快速傅里葉變換后,通過(guò)定位和估值2個(gè)環(huán)節(jié)重構(gòu)得到頻域信號(hào)。稀疏傅里葉變換算法的流程如圖1所示。

圖1 稀疏傅里葉變換算法流程

1.1 信號(hào)的稀疏表示

x=Fs

(1)

其中,s是信號(hào)x在正交基F上的投影系數(shù)向量。若s中只有K個(gè)系數(shù)不為0,其他N-K個(gè)系數(shù)全為0,那么信號(hào)x在頻域上是K-稀疏的,K表示信號(hào)x的稀疏度。若s中K個(gè)系數(shù)值較大,其他N-K個(gè)系數(shù)很小,且?guī)缀鯙?,則稱(chēng)信號(hào)x在頻域上是近似K-稀疏的。若一個(gè)信號(hào)在頻域上能夠被稀疏表示,則該信號(hào)可進(jìn)行稀疏傅里葉變換,也可以降維到低維空間中。

1.2 時(shí)域降維

如果信號(hào)在頻域能夠被稀疏表示,當(dāng)設(shè)置若干個(gè)桶且桶數(shù)比頻率數(shù)更多時(shí),那么按照一定規(guī)則將各個(gè)有效頻率分到這些桶中,每個(gè)桶中只存在唯一的頻率。從頻域角度考慮,將每個(gè)桶視為一個(gè)序列單元,那么長(zhǎng)序列可被映射為短序列。對(duì)短序列進(jìn)行快速傅里葉變換,則算法復(fù)雜度將大幅度降低。sFFT的時(shí)域降維過(guò)程可以表示為:

y=ΨB×Nx

(2)

將每一桶作為一個(gè)序列單元進(jìn)行處理,其方法是取桶中的一個(gè)有效頻點(diǎn)代表該桶的值,該步驟稱(chēng)為下采樣。下采樣對(duì)x中的元素按索引進(jìn)行等間隔的累加,其處理過(guò)程如式(3)所示。

(3)

為保證取到每一個(gè)有效頻點(diǎn),利用頻域卷積定理,將時(shí)域信號(hào)與窗函數(shù)2個(gè)矢量相乘,則在頻域上各桶內(nèi)唯一存在的有效頻點(diǎn)展寬成多個(gè)頻點(diǎn),這樣在下采樣時(shí)就能夠取到有效頻點(diǎn)。本文所使用的窗函數(shù)取自文獻(xiàn)[11]的平滑窗函數(shù)。

1.3 頻域重構(gòu)

在正交頻分復(fù)用(OFDM)技術(shù)的頻偏估計(jì)方法中,2個(gè)采樣點(diǎn)之間的相位差與頻點(diǎn)坐標(biāo)呈線(xiàn)性關(guān)系。漸近法從相位差角度出發(fā),運(yùn)用該原理中不同頻率所對(duì)應(yīng)相位不同的特點(diǎn),通過(guò)對(duì)相位區(qū)間進(jìn)行能量檢測(cè)判斷有效頻率的存在與否。該類(lèi)方法可對(duì)相位二分查找,或劃分成多區(qū)段查找[17]逐步迭代得到頻率。二分漸近法保證了算法的魯棒性,但復(fù)雜度相對(duì)較高;相比之下多區(qū)段漸進(jìn)法速度較快,復(fù)雜度較低,但魯棒性較差。

統(tǒng)計(jì)法主要從重復(fù)統(tǒng)計(jì)角度出發(fā),利用哈希映射、孫子定理[11]等方法先確定下采樣所對(duì)應(yīng)頻域中存在頻點(diǎn)的桶,反推出原始信號(hào)中可能包含的頻率,通過(guò)改變重排參數(shù)對(duì)哈希逆映射的多次統(tǒng)計(jì),或?qū)Σ煌瑓?shù)下得到的線(xiàn)性同余方程求解,得到原始信號(hào)中的有效頻率。由于統(tǒng)計(jì)法的重排隨機(jī)進(jìn)行,因此會(huì)產(chǎn)生一定的頻率沖突。在每次循環(huán)中存在的虛警也會(huì)產(chǎn)生極少量的重構(gòu)信息誤差[17]。但統(tǒng)計(jì)法均衡了算法復(fù)雜度和魯棒性,該方法性能較為穩(wěn)定,具有較好的研究前景。

2 稀疏度自適應(yīng)的稀疏傅里葉變換算法

為解決未知稀疏度的問(wèn)題,基于sFFT-1算法的優(yōu)勢(shì),本文提出了SAsFFT算法。本節(jié)從算法總體流程、稀疏度估計(jì)的影響因素和稀疏度的確定3個(gè)部分介紹SAsFFT算法。

2.1 算法總體流程

SAsFFT利用了信號(hào)的頻域稀疏性,在信號(hào)稀疏度未知的前提下,通過(guò)迭代逐步估計(jì)得到信號(hào)的稀疏度,從而使算法適用于未知稀疏度的信號(hào)。SAsFFT與現(xiàn)有的sFFT算法相比,其改進(jìn)在于增加了對(duì)稀疏度的估計(jì)和檢測(cè)過(guò)程。該算法可以分為3個(gè)步驟,其總體流程如圖2所示。

步驟3有效分量篩選。對(duì)稀疏傅里葉變換的結(jié)果進(jìn)行能量檢測(cè),剔除無(wú)效的虛警分量,得到最終的k個(gè)傅氏域信息。

圖2SAsFFT算法流程

2.2 稀疏度估計(jì)的影響因素

利用第2.1節(jié)中信號(hào)的稀疏特性,以及稀疏傅里葉變換中時(shí)域降維的性質(zhì)對(duì)頻域稀疏的離散信號(hào)進(jìn)行稀疏度的估計(jì)得到預(yù)測(cè)稀疏度。稀疏度估計(jì)的流程如下:

子步驟1初始化預(yù)測(cè)稀疏度;

子步驟2時(shí)域降維;

子步驟3快速傅里葉變換;

對(duì)預(yù)測(cè)稀疏度準(zhǔn)確性的影響因素主要有2個(gè):頻譜碰撞和窗函數(shù)卷積。本節(jié)主要從以下的2個(gè)方面進(jìn)行分析。

2.2.1 頻譜碰撞的影響

稀疏度決定了一個(gè)信號(hào)變換到下采樣域后的矢量維度,進(jìn)而影響頻譜碰撞概率。首先可推導(dǎo)在未知稀疏度的條件下,設(shè)定稀疏度與碰撞概率的關(guān)系。文獻(xiàn)[11]給出了定理1,該定理限定了非零系數(shù)落入同一桶中的概率。

定理1若j≠0,n是2的冪次方,σ∈[1 ,n]是一個(gè)隨機(jī)奇數(shù),則Pr[σj∈[-C,C]mod n]≤4C/n。

證明:令bi為下采樣域的第i個(gè)元素。考慮任意i,j∈S,由定理1有:

Prσ,τ[bi=bj]≤ Prσ,τ[|(σi+τ)mod n-(σj+τ)mod n|<

n/B]=Prσ,τ[|σ(i-j)|mod n

那么,從原時(shí)域的角度考慮,對(duì)于信號(hào)矢量x有:

得證。

2.2.2 窗函數(shù)卷積的影響

本節(jié)闡明窗函數(shù)的卷積會(huì)對(duì)信號(hào)中有效頻率產(chǎn)生虛警,使得估計(jì)的頻率數(shù)增加,導(dǎo)致預(yù)測(cè)稀疏度大于真實(shí)稀疏度。

根據(jù)文獻(xiàn)[17]中的定理3.3,以及對(duì)本文算法中信號(hào)和窗函數(shù)的假設(shè),對(duì)降維過(guò)程結(jié)果可得如下的推論。

2.3 信號(hào)稀疏度的確定

2.4 算法復(fù)雜度分析

(4)

3 仿真結(jié)果與分析

為了驗(yàn)證SAsFFT的可行性和魯棒性,本節(jié)對(duì)SAsFFT、FFTW和sFFT-1 3種算法進(jìn)行了仿真實(shí)驗(yàn)。仿真實(shí)驗(yàn)分別是在不同信號(hào)長(zhǎng)度、不同稀疏度下的算法對(duì)比,以及SAsFFT算法的魯棒性仿真。本文算法中所采用測(cè)試信號(hào)的產(chǎn)生方式參考了文獻(xiàn)[9]。時(shí)域降維過(guò)程使用的窗函數(shù)為2.2節(jié)中所定義的窗函數(shù)。

信號(hào)參數(shù)設(shè)置如下:信號(hào)x的矢量長(zhǎng)度為n,在信號(hào)的傅氏域上從[1,n]中隨機(jī)選擇k個(gè)分量,各分量的幅度為1,相位隨機(jī),其他分量設(shè)為0。在仿真結(jié)果中,每個(gè)點(diǎn)代表不同參數(shù)的信號(hào),每個(gè)點(diǎn)是經(jīng)過(guò)100次獨(dú)立的重復(fù)實(shí)驗(yàn)的平均結(jié)果。在實(shí)驗(yàn)中,與SAsFFT算法進(jìn)行對(duì)比的是FFTW(the Faster Fourier Transform in the West)和sFFT-1[11]。

實(shí)驗(yàn)平臺(tái)采用Ubuntu Linux 16.04版本,處理器為Intel(R) Core(TM) i5-4210M CPU @ 2.60 GHz,內(nèi)存為8 GB。

3.1 信號(hào)長(zhǎng)度對(duì)不同算法的影響比較

對(duì)于本文實(shí)驗(yàn),固定信號(hào)稀疏度參數(shù)k=50,選擇了8個(gè)不同的信號(hào)長(zhǎng)度分別為n=217,218,…,224。在不同的信號(hào)長(zhǎng)度下,對(duì)同一信號(hào)分別運(yùn)行SAsFFT、 FFTW和sFFT-1,得到3種算法的平均運(yùn)行時(shí)間對(duì)比如圖3所示。

圖3 不同信號(hào)長(zhǎng)度的算法平均運(yùn)行時(shí)間對(duì)比

從圖3中可以發(fā)現(xiàn),在橫坐標(biāo)為對(duì)數(shù)分度時(shí),SAsFFT、FFTW和sFFT-1的運(yùn)行時(shí)間均呈線(xiàn)性關(guān)系,SAsFFT和sFFT-1的斜率基本相同,且比FFTW的斜率小。SAsFFT與sFFT-1相比,在不同信號(hào)長(zhǎng)度下SAsFFT的平均運(yùn)算時(shí)間均相對(duì)較長(zhǎng),但平均運(yùn)算時(shí)間差比較穩(wěn)定。同時(shí),由圖3可知,若信號(hào)稀疏度k=50,當(dāng)信號(hào)長(zhǎng)度n>219時(shí),則SAsFFT運(yùn)算速度比FFTW快。信號(hào)長(zhǎng)度越長(zhǎng),SAsFFT的平均運(yùn)算時(shí)間越短,其運(yùn)算速度優(yōu)勢(shì)體現(xiàn)就越明顯。

3.2 信號(hào)稀疏度對(duì)不同算法的影響比較

對(duì)于本文實(shí)驗(yàn),固定信號(hào)長(zhǎng)度參數(shù)n=222,選擇了11個(gè)不同的信號(hào)稀疏度k=50,100,200,…,1 000。在不同的信號(hào)稀疏度下,對(duì)同一信號(hào)分別運(yùn)行SAsFFT、FFTW和sFFT-1,得到3種算法的平均運(yùn)行時(shí)間對(duì)比如圖4所示。

圖4 不同信號(hào)稀疏度的算法平均運(yùn)行時(shí)間對(duì)比

從圖4可以觀(guān)察到,FFTW在不同信號(hào)稀疏度下運(yùn)算時(shí)間基本保持不變,而SAsFFT和sFFT-1的平均運(yùn)算時(shí)間均呈線(xiàn)性關(guān)系。在不同信號(hào)稀疏度下SAsFFT的平均運(yùn)算時(shí)間比sFFT-1長(zhǎng),但兩者的平均運(yùn)算時(shí)間差也比較穩(wěn)定。當(dāng)信號(hào)長(zhǎng)度固定不變時(shí),信號(hào)稀疏度為900,SAsFFT與FFTW的平均運(yùn)算時(shí)間基本相等。因此,綜合來(lái)看,對(duì)于小于900的信號(hào)稀疏度,SAsFFT的運(yùn)算時(shí)間更快。

除此之外,圖4給出了不同稀疏度下SAsFFT算法的迭代次數(shù)。當(dāng)稀疏度k≤900時(shí),迭代次數(shù)維持在2次,可見(jiàn)SAsFFT的稀疏度估計(jì)過(guò)程比較穩(wěn)定,即SAsFFT算法的運(yùn)算時(shí)間比較穩(wěn)定,不會(huì)出現(xiàn)較大波動(dòng)。

3.3 信號(hào)稀疏度對(duì)算法魯棒性的影響

(5)

在不同的信號(hào)稀疏度下,對(duì)同一信號(hào)分別運(yùn)行SAsFFT,得到該算法的每符號(hào)平均l1誤差對(duì)比,如圖5所示。

圖5 不同信號(hào)稀疏度的魯棒性對(duì)比

由圖5可以得到,對(duì)于固定長(zhǎng)度的信號(hào),算法不同信號(hào)稀疏度的每符號(hào)平均l1誤差都小于10-6,算法保持了較高的穩(wěn)定性,完全可以滿(mǎn)足絕大部分信號(hào)的實(shí)際應(yīng)用需求。

4 結(jié)束語(yǔ)

本文提出一種基于稀疏度自適應(yīng)的稀疏傅里葉變換算法。該算法根據(jù)頻譜碰撞和窗函數(shù)卷積對(duì)頻率檢測(cè)的影響估計(jì),得到一個(gè)稍大于真實(shí)稀疏度的預(yù)測(cè)稀疏度,通過(guò)已有的稀疏傅里葉變換后,在算法的輸出中剔除冗余數(shù)據(jù),從而得到較好的結(jié)果。根據(jù)數(shù)學(xué)理論分析和算法運(yùn)行的實(shí)驗(yàn)結(jié)果表明,該算法實(shí)現(xiàn)了未知稀疏度條件下的稀疏傅里葉變換,擴(kuò)展了稀疏傅里葉變換算法的應(yīng)用范圍;算法運(yùn)算量小,復(fù)雜度低,在運(yùn)算時(shí)間上快于傳統(tǒng)FFT算法,且完成傅里葉變換后的效果理想,驗(yàn)證了算法的可靠性。

[1] CHANDRAKASAN A,GUTNIK V,XANTHOPOULOS T.Data Driven Signal Processing:An Approach for Energy Efficient Computing[C]//Proceedings of International Symposium on Low Power Electronics and Design.Washington D.C.,USA:IEEE Press,1996:347-352.

[2] DONOHO D L.Compressed Sensing[J].IEEE Transac-tions on Information Theory,2006,52(4):1289-1306.

[3] CANDES E J,WAKIN M B.An Introduction to Com-pressive Sampling[J].IEEE Signal Processing Magazine,2008,25(2):21-30.

[4] 盧 迅,李立春,李 鷗,等.非重構(gòu)壓縮樣值跳頻信號(hào)時(shí)頻聯(lián)合估計(jì)算法[J].信息工程大學(xué)學(xué)報(bào),2016,17(4):425-430.

[5] LINIAL N,MANSOUR Y,NISAN N.Constant Depth Circuits,Fourier Transform,and Learnability[J].Journal of the ACM,1989,40(3):607-620.

[6] DAUBECHIES I,RUNBORG O,ZOU A J.A Sparse Spectral Method for Homogenization Multiscale Problems[J].Siam Journal on Multiscale Modeling & Simulation,2007,6(3):711-740.

[7] 范星宇.低軌衛(wèi)星星載通信信號(hào)處理關(guān)鍵技術(shù)研究[D].北京:北京理工大學(xué),2014.

[8] GILBERT A C,GUHA S,INDYK P,et al.Near-optimal Sparse Fourier Representations via Sampling[C]//Proceed-ings of ACM Symposium on Theory of Computing.New York,USA:ACM Press,2002:152-161.

[9] GILBERT A C,MUTHUKRISHNAN S,STRAUSS M.Improved Time Bounds for Near-optimal Sparse Fourier Representations[C]//Proceedings of SPIE’05.Washington D.C.,USA:IEEE Press,2005:398-412.

[10] IWEN M A,GILBERT A,STRAUSS M.Empirical Evaluation of a Sub-linear Time Sparse DFT Algorithm[J].Communica-tions in Mathematical Sciences,2007,5(4):981-998.

[11] IWEN M A.Combinatorial Sublinear-time Fourier Algori-thms[J].Foundations of Computational Mathematics,2010,10(3):303-338.

[12] HASSANIEH H,INDYK P,KATABI D,et al.Simple and Practical Algorithm for Sparse Fourier Transform[C]// Proceedings of ACM-SIAM Symposium on Discrete Algorithms,Society for Industrial and Applied Mathematics.New York,USA:ACM Press,2012:1183-1194.

[13] 孟祥玉,李 靜,彭 華.一種基于迭代更新的稀疏傅里葉變換改進(jìn)算法[J].信息工程大學(xué)學(xué)報(bào),2016,17(6):669-674.

[14] INDYK P,KAPRALOV M.Sample-optimal Fourier Sampling in Any Constant Dimension[C]//Proceedings of the 55th IEEE Annual Symposium on Foundations of Computer Science.Philadelphia,USA:IEEE Press, 2014:514-523.

[15] WANG C,MAURICIO A P,CHANDRASEKARAN S,et al.Parallel Sparse FFT[C]//Proceedings of Workshop on Irregular Applications:Architectures and Algorithms.New York,USA:ACM Press,2013:10-23.

[16] 周廣福,文成林,高敬禮.基于小波變換與稀疏傅里葉變換相結(jié)合的光場(chǎng)重構(gòu)方法[J].電子學(xué)報(bào),2017,45(4):782-790.

[17] HAITHAM H,INDYK P,KATABI D,et al.Nearly Optimal Sparse Fourier Transform[C]//Proceedings of the 44th ACM Symposium on Theory of Computing.New York,USA:ACM Press,2012:563-578.