低采樣率高分辨率壓縮功率譜估計方法的仿真研究

姚彥鑫

(北京信息科技大學,北京100010)

低采樣率高分辨率壓縮功率譜估計方法的仿真研究

姚彥鑫

(北京信息科技大學,北京100010)

低采樣率的寬帶功率譜估計在很多領域具有應用價值．采用壓縮多核采樣結構得到信號的壓縮測量值,然后建立測量值相關函數與信號相關函數之間的關系,用最小二乘法實現相關函數估計,最后實現功率譜的估計．該壓縮采樣方法的等效采樣率為M/N·fs,可在沒有任何對時域或頻域稀疏性的假設條件下降低采樣率．仿真分析表明,該方法的系統噪聲與加性噪聲性能比周期圖法略有降低,但只要系統設計合理,對于一定信噪比的信號,系統噪聲與加性噪聲基本可以忽略,并給出了對應的理論分析．估計分辨率與周期圖法相比,等效長度相同時略有提高;由于本文方法降低了測量值的數目,對于一定長度的數據來說,估計分辨率得到了極大的提高．本文方法適用于低信噪比信號的低采樣率高分辨率功率譜估計．

功率譜估計;降采樣率;自相關;頻率分辨率;信噪比

引 言

功率譜密度是隨機信號自相關序列的離散傅里葉變換,用來表征隨機信號的統計特征,可以用于諧波檢測、噪聲頻譜測量、弱信號檢測等．在雷達、聲納、通信、地址勘探、天文、生物醫學工程領域具有廣泛的應用[1-4]．在通信領域可以應用于多載波通信中的解調環節,或者在突發通信模式中進行頻率估計,從而保證可靠通信[2]．

在高頻雷達應用中,通過回波信號的功率譜密度分析可以得到有關目標的各種信息:根據功率譜密度峰值的寬度來判斷運動目標的位置;根據峰值的高度判斷運動目標的強度;根據峰值的位置來判斷目標的運動速度等[5]．在雷達信息處理中,多普勒頻率估計及多普勒頻率域處理是鑒別多目標的一種主要方法[6]．在電子對抗領域,雷達信號的頻率信息是信號分選、威脅識別、引導干擾的重要參數,如何對截獲雷達信號載頻進行高精度估計一直是電子戰接收機的設計重點．大規模集成電路的發展為數字測頻提供了硬件平臺,模擬數字混合接收機既可以進行實時處理,也可以將信號存儲起來,數字測頻的核心在于算法．

信號功率譜的估計方法可分為兩類:一為線性估計方法,有自相關估計、自協方差法及周期圖法等[7]．另一類為非線性估計方法,有最大似然法(MaximumLikelihood,ML)、最大熵法等．線性估計方法的譜分辨率隨數據長度的增加而提高,如周期圖方法等,然而這類方法在快速高分辨率估計中遇到了困難,難以實現多普勒的高精度估計[5]．快速傅里葉變換線性估計方法估計頻率也面臨頻率分辨率與處理數據量之間存在矛盾的問題[8]．非線性估計方法大多是無偏的譜估計方法,可以獲得高的譜分辨率．ML在所有方法中精度最高,信噪比性能最好,但需要解計算量超高的非線性問題,很多情況無法得到最優解,實時性很差,不適合實際應用．另外兩種非線性方法,多重信號分類算法(MultipleSignalClassification,MUSIC)可以達到高的頻率分辨率,但在對整個頻域內搜索峰值耗時較多,實時性差[3];子空間旋轉不變技術 (EstimatingSignalParametersviaRotationalInvarianceTechnique,ESPRIT)計算量比MUSIC小,但是精度沒有MUSIC高[4]．

實時信號處理希望處理的數據量越少越好,而高的頻率分辨率卻對數據量提出相反的要求．如果能夠降低所需采樣數據的數目,同時能夠以較高分辨率分辨信號,就成為理想的選擇．這方面的文章有基于壓縮感知理論構造冗余字典的方法等[9],但此類方法一般需要壓縮感知重構算法解算,計算量較大．文獻[1, 10]給出一種寬帶頻譜的功率譜估計方案,其性能分析側重于寬帶譜感知的估計和檢測性能方面[11-12]．本文則側重于頻點估計的應用,進行了頻率分辨率、噪聲性能等方面的分析,并給出了在最小稀疏尺度準則的可選方案中,可借鑒的參數方案的選取方法,得到了兼顧采樣數據數目、抗噪聲性能與頻率分辨率的參數方案,為該方法在頻率估計中應用提供重要參考．仿真發現該方法在降低采樣率的同時,頻率分辨率得到提高．至于噪聲性能,只要系統設計合理,該方法能夠處理一定程度低信噪比的信號．

1 壓縮采樣

設x(t)是一帶寬為1/T的復值廣義平穩信號．圖1是一個被稱為模擬信息轉換(AnalogInformationConversion,AIC)的采樣裝置[1],能夠對模擬信號進行壓縮采樣．該裝置有M個分支,第i個分支包含一周期為NT的隨機波形調制器和一周期也為NT的積分鎖存器,則第i支的第k個輸出結果可表示為

(1)

式中,ci(t)是偽隨機調制波形一個周期的信號pi(t)乘以1/N．圖1實現了對原始信號的壓縮,使得觀測值含有了原始信號的所有信息．

圖1 信號壓縮采樣框圖

在數字仿真模擬數字轉換采樣過程的平臺上,式(1)可表示為

(2)

取調制波形為分段常數函數,滿足ci(t)=ci[-n],nT

2 功率譜估計的時域重構方法

下面基于所獲得的{yi[k]}i,k樣本,進行x[n]的功率譜重構．x[n]的功率譜或功率譜密度表達式如下:

(3)

定義rx(n)=E(x[m]x*[m-n])為x(n)的自相關函數．我們可利用M2個yi[k]與yj[k]的互相關值估計自相關函數rx(n),實現壓縮估計而不需要任何關于x(t)的時域或頻域稀疏約束條件．

yi[k]與yj[k]的互相關表示為

(4)

yi是一個對zi的N折疊抽取,即yi[k]=zi[kN]．

rzi,zj[n] =rci,cj[n]?rx[n]

(5)

式中,?表示卷積運算．由式(5)可得

(6)

式中:

rci,cj[0]= [rci,cj[0],rci,cj[-1],…,

rci,cj[-N+1]]T;

(7)

rci,cj[1]= [rci,cj[N],rci,cj[N-1],…,

rci,cj[1]]T;

(8)

rx[k]= [rx[kN],rx[kN+1],…,

rx[(k+1)N-1]]T．

(9)

rci,cj[l]為ci[·]信號的自相關;

通過級聯M2個不同的互相關函數ryi,yj[k],i,j=0,1,…,M-1,得到了M2×1的矢量ry[k]=[…,ryi,yj[k],…]T．根據式(7)、(8)、(9)可推出

(10)

式中,Rc[0]=[…,rci,cj[0],…]T,Rc[1]=[…,rci,cj[1],…]T是M2×N的矩陣．雖然x(n)是帶限的,ry會有無限的支撐,但在許多實際應用中,ry[k]只有在-L≤k≤L范圍內具較大的值,在此范圍外都是可忽略的．又由式(10)知,ry[k]取決于rx[k]與rx[k-1],則可認為rx[k]也需要被限制在-L≤k≤L．那么所有的有效觀測量為:

(11)

(12)

注意以下兩點:

1)據上討論,rx[k]被限制于-L≤k≤L,且rx[n]有復共軛對稱性,所以rx[n]被限制在-NL≤n≤NL,rx[L]的最后N-1個值都是0．

2)根據rci,cj[l]的定義,rci,cj[n]被限制在1-N≤n≤N-1,Rc[1]的第一列是0．

這兩點使得可以將式(10)中的線性卷積表達為在-L≤k≤L上沒有任何額外補零的循環卷積.因此,我們最終將ry以及rx的關系表示為

ry=Rcrx．

(13)

式中,Rc是[(2L+1)M2×1]×[(2L+1)N×1]的矩陣,

(14)

如果Rc是滿列秩的,即要求M2≥N,則式(14)可以利用最小二乘法求解．

最后利用式(15)求出(2L+1)N×1階的能量譜矢量Sx:

=F(2L+1)Nrx.

(15)

式中,F(2L+1)N是(2L+1)N×1的離散傅里葉矩陣．

3 最小稀疏尺度抽樣方案

為了保證最小二乘(Least Squares,LS)解的唯一性,隨機調制波形可以有很多不同的實現方法,如高斯抽樣、貝努力抽樣等．研究如何選取隨機調制波形,才能在保證估計性能的條件下,使支路的數目M盡量?。?/p>

本文采用最小稀疏尺度抽樣多核實現方法．該方法基于最小稀疏尺度問題,可以通過簡單地設置ci[n]來實現,對于不同的i,設置不同的ci[n]．ci[n]=1如果-n=ni;ci[n]=0,如果-n≠ni,其中i≠j時,ni≠nj．

以上ci[n]的設置方法也可以由從單位陣IN中選擇M個不同的行實現．然而,注意到這種行的選擇不會是隨機的,因為我們要保證式(14)中Rc的列滿秩性,其實只需保證Rc中每列至少有一個1．

而由自相關函數的對稱性和式(7)、式(8)知,如果Rc[0]對應于前[N/2]+1列至少有一個1,那么對應地,Rc[1]于后[N/2]列至少有一個1,[x]表示不大于x的最大整數．因此,問題成為選擇IN的行組合,保證Rc[0]在其前[N/2]+1列至少有一個1．為了最小化壓縮比M/N,希望M盡量小．

(16)

式中,|S|表示集合S的勢．

盡管式(16)可通過窮舉或貪婪搜索法求解,但運算量太大．一個可能找到式(16)次優解的方法是通過解所謂的最小長度[N/2]稀疏尺度問題,該問題已經在文獻[13-14]中被很好地研究．表1中是一些符合最小稀疏準則的可選M、N對的取值．

本文方法的優點是將N個采樣點壓縮為M個支路的測量值,壓縮比為M/N,達到壓縮功率譜估計的目標．并且直接在模擬域實現信號變換,而后進

表1 最小稀疏準則的例子

行壓縮采樣,大大降低了數字信號的采樣率和處理復雜度,求解的時候,除了必要的傅里葉變換矩陣轉換,即時域與頻域的轉換外,只用到了LS方法,避開了壓縮感知恢復算法等復雜的運算環節,是功率譜估計一個可選方法．因此,本方法中只需要滿足M2≥N和Rc矩陣的列滿秩,即成為過定方程,可用最小二乘法求解．

4 壓縮功率譜估計方法的性能

4.1 分辨率性能分析

功率譜估計的分辨率是分辨兩個頻率相近的正弦(或窄帶)信號的能力,常用功率譜的主瓣寬度來度量,即主瓣在半功率點處的半邊寬度．

仿真方案:信號為單頻信號x=10cos(2πf0t),最小稀疏尺度抽樣多核采樣方案,參數選取為(M=6、N/2=13),隨機調制波形的碼速率是1/T=800 Hz．周期圖法中,參數為fs=800 Hz．為了方便比較,兩種方法處理信號的時間長度均為N(L+1)T．對于同樣時間的采樣數據,采樣數目不同,周期圖法是N(L+1),本文方案是M(L+1)．

兩種方法分別對參數f0=100 Hz、200 Hz、300 Hz,L=5、10、20、30的每對組合進行100次試驗．統計本文方法和周期圖法功率譜估計分辨率的平均值,如表2所示．圖2為某次實驗中周期圖法的頻率分辨率的計算示意圖．

表2 兩種估計方法的頻率分辨率對比(M=6、N/2=13)

圖2 周期圖法的頻率分辨率L=30

由表2可知:在同樣測量數目的情況下,本文方法的分辨率明顯高于周期圖法;而對于同樣的等效數據長度,本文方法的分辨率也略高于周期圖法．

為了比較最小稀疏尺度抽樣多核采樣方案在不同參數設置下分辨率的情況,選取L=20,對于f0=100 Hz、200 Hz、300 Hz,不同的采樣方案(M、N/2)分別進行100次試驗,統計分辨率如表3所示．觀察得到當L一定時,最小稀疏尺度抽樣多核采樣方案使用實際數據的長度越長,分辨率越高．

周期圖的頻率分辨率嚴格與所使用數據的長度成反比[15]．而本文方法當(M,N/2)取值一定的時候,分辨率隨著L長度的增加而近似成反比例規律提高,如表2所示．實驗結果不成嚴格反比關系可能是因為仿真次數較為有限,結果存在一些偏差．本文方法當L取值一定,參數(M,N/2)變化時,分辨率與等效數據長度表現出類似于反比例的關系,如表3所示．

表3 最小稀疏尺度抽樣多核采樣方案在不同參數方案的頻率分辨率的對比

4.2 噪聲條件下性能分析

衡量抗系統噪聲性能指標的實驗方案:取有用信號為x=10·cos(2πf0t)．為考察最小稀疏尺度抽樣多核采樣方案的噪聲性能,以L=30,(M=8,N/2=23)為例,1/T=800 Hz,那么測量數據長度為(L+1)M,等效數據長度為(L+1)N．周期圖法采樣率也為800 Hz,采樣數據長度為(L+1)N．未在有用信號x上添加噪聲時,計算有用頻率處的譜峰值是p,代表有用信號的功率;計算其他位置處譜峰的最大值是q,代表該頻率處的系統噪聲功率,衡量抗系統噪聲性能用10log2(p/q)dB．如圖3(a)一次試驗的局部圖所示,抗系統噪聲性能為13.05 dB．

為了驗證最小稀疏尺度抽樣多核采樣方案取不同最小稀疏尺度時的有效性和抗噪聲性能,隨機選用0～400 Hz之間的單頻點作為f0,用不同隨機波形的選取方法,反復進行實驗,得到抗系統噪聲指標的平均值．為了公平,每組參數的第二次試驗中,調整L使得等效的數據長度即(L+1)N接近相同．

由表4可見,除了(M=5,N/2=9)組合外,基本上抗自噪聲性能都在13 dB左右．那么,綜合壓縮比的因素,得到M=8,N/2=23時,信號的壓縮比與系統自噪聲性能都相對較好．觀察得到新方法的抗自噪聲性能隨著信號數據增加,沒有明顯改善．這與傳統周期圖法的規律是不同的．考慮主要是由于方法在L=30時數據自相關信息已經基本包括,L不需要取得更長．

(a) 局部圖

(b) 全局圖圖3 壓縮功率譜估計系統噪聲性能 (L=30,M=8,N=46)

周期圖法的抗系統噪聲性能:f0從0～400 Hz隨機選取,采樣率800 Hz,統計100次試驗中最大相鄰旁瓣與峰值之差的平均值得到抗系統噪聲性能,在18 dB左右,較為理想．f0=100 Hz的示例如圖4所示．

衡量抗加性噪聲性能的實驗方案:有用信號x=a·cos(2πf0t),在有用信號x基礎上添加了加性噪聲,該加性噪聲是均值為0、方差為σ2的限帶高斯白噪聲,噪聲覆蓋帶寬為0～800 Hz．周期圖法采用fs=800 Hz,使用了(31·N)的采樣數據長度．多核采樣方法參數采用M=6、N/2=13、L=30．

圖4 周期圖法的抗系統噪聲性能

M[N/2]L選擇的尺子抗系統干擾性能/dB3330116[0,1,3]14.9613.15593043[0,1,2,5,9]8.909.006133053[0,1,4,5,11,13]13.1013.2182330[0,1,2,11,15,18,21,23]13.00

兩種方法都隨機選用0～400 Hz之間的單頻作為f0,當正弦信號幅度a和噪聲方差σ2取不同值時,反復進行Monte Carlo實驗,得到的抗加性噪聲性能如表5所示．將兩種方法得到的功率譜有用頻率處譜峰值,與代表加性噪聲最大譜峰值的比值作為抗加性噪聲性能指標．

表5 兩種方法的抗加性噪聲性能對比

從表5可以看出,新方法相比于周期圖法,其抗系統噪聲性能略有降低．這是因為經過理論推導,發現周期圖法等效于使用自相關函數加上三角Bartlett窗來計算功率譜[15],而新方法相當于使用自相關函數加上近似為矩形的窗來估計功率譜．由于抗系統噪聲性能主要取決于所使用窗函數第一旁瓣與主瓣的高度比值,而兩種方法分別使用的窗函數頻譜的主瓣與第一旁瓣比值為:矩形窗約為9dB,Bartlett窗為17dB．這個理論分析的結果與表5中仿真得到的兩種方法的平均抗系統噪聲11.80dB和18.48dB基本吻合．

兩種方法的抗加性噪聲的性能隨信噪比的降低而減弱．新方法的抗加性噪聲的性能并沒有隨著噪聲強度變化成正比變化,這是因為根據式(10),當x(t)包括噪聲時,估計自相關中會包括無噪信號的自相關、噪聲自相關與無噪信號和噪聲的互相關三項．仿真發現,由于所使用估計信號長度有限,限帶“白噪聲”的互相關也并不是完美地為“0”,加上噪聲與信號的互相關項的影響,使得估計出來的自相關函數隨著噪聲能量的增加,并不呈現明顯的線性關系．自相關函數再經過傅里葉變換得到功率譜(見式(15)),則新方法的抗加性噪聲性能,雖然具有單調性,但呈現出較為復雜的非線性關系,如表5所示．

雖然新方法的抗噪聲性能不如周期圖法,但能夠降低采樣數據的數目,且分辨率佳．新方法的測量值數目壓縮比可高達0.27．另外,新方法的頻率分辨率不隨測量數據長度的減少成正比例關系下降,降低較少,抗噪聲性能在大多數應用場合能夠適應．因此,對于很多應用,最小稀疏尺度抽樣多核采樣方案是具有實用性的,可以達到對一定程度低信噪比的信號在總數據處理量增加不多的前提下,提高頻率分辨率的效果．

5 結 論

1) 本文方法能夠在對時域和頻域都沒有稀疏性假設條件的前提下,以較低的采樣率,實現低信噪比信號的功率譜估計．

2) 系統的噪聲性能,包括系統噪聲與加性噪聲,主要取決于系統參數．系統噪聲雖然較大,但通過對基于最小稀疏尺度準則的采樣結構的研究,發現只要參數設計合理,就能夠完全克服這個問題;當選取參數合理時,抗加性噪聲性能也可以接受．

3) 本文方法的估計分辨率與相同等效長度的周期法相比,分辨率略有提高;由于本文方法降低了測量值的數據長度,那么對于一定長度的測量數據來說,估計分辨率得到了極大的提高．

4) 對于基于最小稀疏準則的多核采樣方案,如何設計參數能夠取得更好的性能還有待進一步研究．

[1] ARIANANDA D D, LEUS G. Compressive wideband power spectrum estimation[J]. IEEE transactions on signal processing, 2012, 60(9): 4775-4789.

[2] 石盛超, 索龍龍, 張衛同. 基于自相關函數的前向頻率估計方法綜述[J].電子質量, 2012(12):38-42.

SHI S C, SUO L L, ZHANG W T. Autoeorrelation-based feedforward carrier frequency estimator[J]. Electronics quality, 2012(12):38-42. (in Chinese)

[3] 蔣瓊. 基于MUSIC算法的電力諧波、間諧波頻譜估計方法研究[D]. 長沙:長沙理工大學, 2012.

JIANG Q. The research of spectrum estimation for the electrical harmonic and interharmonic based on MUSIC algorithm[D]. Changsha: Changsha University of Science & Technology, 2012. (in Chinese)

[4] 司偉建, 藍曉宇, 劉學. 提高二維DOA估計分辨率的改進MUSIC算法[J]. 西安電子科技大學學報(自然科學版), 2012, 39(2):138-144.

SI W J, LAN X Y, LIU X. Modified MUSIC algorithm for improving the resolution of 2-D direction of arrival estimation[J]. Journal of Xidian University(natural science), 2012, 39(2): 138-144. (in Chinese)

[5] 梁宏光. 基于壓縮感知的頻率估計技術的研究[D]. 哈爾濱:哈爾濱工業大學, 2013.

LIANG H G. Research on frequency estimation on techniques based on compressed sensing[D]. Harbin: Harbin Institute of Technology, 2013. (in Chinese)

[6] 劉銀恩. 高精度頻率估計算法研究[D]. 南京: 南京理工大學, 2007.

LIU Y E. Study on high precision frequency estimation algorithm[D]. Najing: Nanjing University of Science and Technology, 2007. (in Chinese)

[7] 單涼, 王偉, 黃超. 列車測速中多普勒頻率估計方法研究[J].電子技術, 2016, 45(2): 1-5.

SHAN L,WANG W, HUANG C. Train speed estimation method of Doppler frequency [J]. Electronic technology, 2016, 45(2):1-5. (in Chinese)

[8] 黃翔東, 王越冬, 靳旭康, 等. 無窗全相位FFT/FFT相位差頻移補償頻率估計器[J]. 電子與信息學報, 2016, 38 (5): 907-912.

HUANG X D, WANG Y D, JIN X K, et al. N o-windowed apFFT/FFT phase difference frequency estimator based on frequency-shift & compensation[J]. Journal of electronics & information technology, 2016, 38 (5): 907-912.(in Chinese)

[9] 沈志博, 董春曦, 黃龍, 等. 一種基于稀疏分解的窄帶信號頻率估計算法[J]. 電子與信息學報, 2015, 37(4): 907-912.

SHEN Z B, DONG C X, HUANG L, et al. A frequency estimation algorithm of narrow-band signal based on sparse decomposition[J]. Journal of electronics & information technology, 2015, (4):907-912. (in Chinese)

[10]LEUS G, ARIANANDA D D. Power spectrum blind sampling[J]. IEEE signal processing letters, 2011, 18(8): 443-446.

[11]盧光躍, 彌寅, 包志強, 等. 基于特征結構的頻譜感知算法[J]. 西安郵電大學學報, 2014, 9(2): 1-10.

LU G Y, MI Y, BAO Z Q, et al. Cooperative spectrum sensing algorithms based on eigenvalue structure of the received signal[J]. Journal of Xi’an University of Posts and Telecommunications, 2014, 19(2): 1-10. (in Chinese)

[12]劉學謙. 支持向量機在頻率估計算法中的應用研究[D]. 鄭州: 解放軍信息工程大學, 2014.

LIU X Q. Research on application of support vector machoine in frequency estimation algorithm[D]. Zhengzhou: The PLA Information Engineering University, 2014. (in Chinese)

[13]Sparse ruler [EB/OL]. [2016-08-20]. https://en.wikipedia.org/wiki/Sparse-ruler, 2013. 01.

[14]LEECH J. On the representation of 1, 2… n by differences[J]. Journal of the London Mathematical Society, 1956, 31: 160-169.

[15]姚天任, 江太輝. 數字信號處理[M]. 2版. 武漢: 華中科技大學, 2000.

Simulation on low sampling rate high resolution compressed power spectrum estimation method

YAO Yanxin

(BeijingInformationandScienceTechnologyUniversity,Beijing100010,China)

Low sampling rate power spectrum estimation could be applied in many domains. In the paper, firstly the compressed multi-coset sampling structure is adopted to obtain the compressed measuring values, and the relationship between correlation of measuring values and autocorrelation is built. Secondly, the estimation for signal autocorrelation is realized using least squares. At last, the power spectrum estimation is realized through frequency domain transformation. To reduce the compression rate, the realization structure based on minimal sparse rule is studied. The equivalent sampling rate for the method isM/N·fs, which enables low sampling rate spectrum estimation without any sparse assumptions about the frequency or time domain signals. Through simulations, it proves that the system noise and additive noise performance is not as good as periodogram method. But if system design parameters are properly designed, the noise could be ignored. The corresponding theoretical analysis is given as well. The frequency resolution performance is improved compared to periodogram method, however, the method reduces the number of measured data, so for certain measured data, the frequency resolution performance is elevated greatly. Thus, the method is applicable to the low sampling power spectrum estimation of low SNR signals.

power spectrum estimation; low sampling rate; autocorrelation;frequency resolution;signal-to-noise ratio

10.13443/j.cjors.2016082001

2016-08-20

國家自然科學基金(61302073)；北京市自然科學基金(4172021)；北京市教委面上項目；北京市自然科學基金(Z160002)

TP391

A

1005-0388(2016)06-1172-08

姚彥鑫 (1982-),女,河北人,副教授,碩士生導師,研究方向為衛星導航、認知無線電．

姚彥鑫. 低采樣率高分辨率壓縮功率譜估計方法的仿真研究[J]. 電波科學學報,2016,31(6):1172-1179.

YAO Y X. Simulation on low sampling rate high resolution compressed power spectrum estimation method [J]. Chinese journal of radio science,2016,31(6):1172-1179.(in Chinese). DOI: 10.13443/j.cjors.2016082001

聯系人： 姚彥鑫 E-mail: yanxin_buaa@126.com

DOI 10.13443/j.cjors.2016082001