黃 超，王 偉，單 涼

(電子工程學(xué)院，合肥 230037)

?

基于Burg算法功率譜估計(jì)的參數(shù)選擇

黃 超，王 偉，單 涼

(電子工程學(xué)院，合肥 230037)

介紹了現(xiàn)代功率譜估計(jì)中常用的基于自回歸(AR)模型的Burg算法，分析了Burg算法中參數(shù)選擇對(duì)功率譜估計(jì)的影響，并給出綜合最優(yōu)的Burg算法參數(shù)選擇方法，最后將本文參數(shù)選擇方法應(yīng)用于列車測(cè)速系統(tǒng)進(jìn)行實(shí)驗(yàn)驗(yàn)證，準(zhǔn)確地估計(jì)出了列車的多普勒頻率。

功率譜估計(jì)；自回歸模型；Burg算法；參數(shù)選擇

0 引 言

在現(xiàn)代信號(hào)處理中，對(duì)于具有各態(tài)歷經(jīng)性的平穩(wěn)隨機(jī)信號(hào)很難用數(shù)學(xué)關(guān)系式來清楚地描述。為了分析平穩(wěn)隨機(jī)信號(hào)，可以利用采集的N個(gè)樣本數(shù)據(jù)來估計(jì)該信號(hào)的功率譜密度(PSD)，這種方法就被稱為功率譜估計(jì)[1-3]。功率譜估計(jì)常被應(yīng)用于許多實(shí)際問題中，例如雷達(dá)信號(hào)處理和生物醫(yī)學(xué)工程等。一般來說, 功率譜估計(jì)方法可分為兩大類：經(jīng)典譜估計(jì)法(非參數(shù)估計(jì)法)和現(xiàn)代譜估計(jì)法(參數(shù)估計(jì)法)。經(jīng)典譜估計(jì)法分為直接法和間接法，均是以數(shù)字傅里葉變換(DFT)為基礎(chǔ)，具有分辨率低、頻譜混疊等固有缺點(diǎn)?，F(xiàn)代譜估計(jì)法先建立功率譜估計(jì)模型，然后利用觀測(cè)數(shù)據(jù)對(duì)建立的模型進(jìn)行參數(shù)估計(jì)，最后完成功率譜估計(jì)。該方法分辨率很高，擁有比經(jīng)典譜估計(jì)法更好的參數(shù)估計(jì)性能。

目前常用的功率譜估計(jì)模型有自回歸(AR)模型、滑動(dòng)平均(MA)模型和自回歸滑動(dòng)平均(ARMA)模型[4]。其中，AR模型是現(xiàn)代譜估計(jì)中最常用的一種功率譜估計(jì)模型，這是因?yàn)锳R模型的實(shí)際物理系統(tǒng)是全極點(diǎn)系統(tǒng)；而且該模型的參數(shù)估計(jì)算法是線性方程組，運(yùn)算量較小，計(jì)算比較簡(jiǎn)便。在基于AR模型進(jìn)行功率譜估計(jì)時(shí)，必須先提取出AR模型的參數(shù)。目前這些參數(shù)的提取算法主要有Levinson-Durbin算法[5]、Burg算法[6-7]和Marple算法[8]3種。其中，Levinson-Durbin算法的復(fù)雜度最低，而且能夠保證預(yù)測(cè)誤差濾波器最小相位，但是分辨率較低；Marple算法分辨率很高，性能最好，但是對(duì)硬件要求較高，計(jì)算量較大；Burg算法分辨率較高，計(jì)算也不太復(fù)雜，綜合性能最好，常常被用來進(jìn)行功率譜估計(jì)。

基于Burg算法對(duì)實(shí)際工程問題進(jìn)行功率譜估計(jì)時(shí)，如何構(gòu)建適合實(shí)際問題的Burg算法參數(shù)模型至關(guān)重要。本文首先通過MATLAB分析采樣頻率、采樣點(diǎn)數(shù)、AR模型階數(shù)和快速傅里葉(FFT)長(zhǎng)度等Burg算法參數(shù)對(duì)功率譜估計(jì)的影響，并給出綜合最優(yōu)的Burg算法參數(shù)選擇方法，最后采用列車測(cè)速系統(tǒng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行仿真驗(yàn)證。

1 基于Burg算法的功率譜估計(jì)

1.1 基于AR模型的功率譜估計(jì)

AR模型又稱為自回歸模型,是一個(gè)全極點(diǎn)的模型,可用如下差分方程來表示:

(1)

式中：p為AR模型的階數(shù)；{a(k)|k=1,2,…,p}，為p階AR模型的參數(shù)。

將該模型記為AR(p)，它的功率譜為：

(2)

(3)

將式(1)乘以x(n-m) 求平均(數(shù)學(xué)期望)，可以求得觀測(cè)數(shù)據(jù)的AR(p)模型參數(shù)與自相關(guān)函數(shù)的關(guān)系式為:

(4)

由式(4)可知自相關(guān)函數(shù)Rxx具有遞推的性質(zhì)，將式(4)寫成矩陣形式得:

最后，大力發(fā)展科技服務(wù)中介機(jī)構(gòu)。加快落實(shí)《市政府辦公室關(guān)于加快科技服務(wù)業(yè)發(fā)展的實(shí)施意見》要求，重點(diǎn)培育市場(chǎng)化、專業(yè)化的人才中介、研究開發(fā)、技術(shù)轉(zhuǎn)移、檢驗(yàn)檢測(cè)認(rèn)證、科技金融等科技中介服務(wù)機(jī)構(gòu)，加強(qiáng)對(duì)責(zé)任部門實(shí)施情況績(jī)效考核。著力打造“科技領(lǐng)軍型人才創(chuàng)新驅(qū)動(dòng)中心”服務(wù)品牌，完善領(lǐng)軍型人才驅(qū)動(dòng)中心服務(wù)職能建設(shè)，提升政府招才引才育才工作社會(huì)影響力。

(5)

式(5)就是著名的Yule-Walker(Y-W)方程。式(5)表明，只要獲取觀測(cè)數(shù)據(jù){x(0),x(1),…,x(N-1)}的自相關(guān)函數(shù)，就能估計(jì)出AR模型的參數(shù){a(k)|k=1,2,…,p}和σ2，進(jìn)而按式(3)求得信號(hào)功率譜的估值。目前AR模型參數(shù)的提取算法主要有Levinson-Durbin算法、Burg算法和Marple算法3種，本文主要研究基于Burg算法的功率譜估計(jì)。

1.2 Burg算法

Burg算法的基本思想是基于線性預(yù)測(cè)器的前、后向預(yù)測(cè)的總均方誤差之和最小的準(zhǔn)則直接從觀測(cè)數(shù)據(jù)來估計(jì)反射系數(shù)，然后通過Lenvinson-Durbin算法的遞推公式求出AR模型參數(shù)。這種方法的優(yōu)點(diǎn)就是對(duì)未知數(shù)據(jù)不需要做任何假設(shè)，估計(jì)精度較高。假設(shè)觀察到的N個(gè)數(shù)據(jù)為X(0),X(1),X(2)，…，X(N-1)，則具體算法如下：

(1) 初始化前、后向預(yù)測(cè)誤差以及預(yù)測(cè)誤差功率：

(6)

(7)

(2) 計(jì)算反射系數(shù):

m=1,2，…，p

(8)

(3) 計(jì)算濾波器系數(shù)及預(yù)測(cè)誤差功率:

am(m)=ρm

(9)

am(k)=am-1(k)+ρmam-1(m-k)，

k=1,2,…，m-1

(10)

(11)

(4) 遞推高一階前、后向預(yù)測(cè)誤差:

(12)

(13)

把m更新為m+1，重復(fù)步驟(2)至步驟(4),直到預(yù)測(cè)誤差功率Em滿足要求。

2 Burg算法的參數(shù)選擇分析

為了分析Burg算法的參數(shù)選擇對(duì)功率譜估計(jì)的影響，本文首先構(gòu)建仿真信號(hào)為：

x(t)=5sin(2πf1t)+3sin(2πf2t)+n(t)

(14)

2.1 采樣頻率的選擇對(duì)功率譜估計(jì)的影響

為了分析Burg算法中采樣頻率的選擇對(duì)功率譜估計(jì)的影響，必須保證采樣點(diǎn)數(shù)、AR模型階數(shù)和FFT長(zhǎng)度恒定?，F(xiàn)令采樣點(diǎn)數(shù)N=512，AR模型階數(shù)p=100，F(xiàn)FT長(zhǎng)度為2 048，采樣頻率分別選擇500 Hz,1 kHz,5 kHz和10 kHz。對(duì)信號(hào)x(t)進(jìn)行功率譜估計(jì)，仿真結(jié)果如圖1所示。

圖1 采樣率不同時(shí)的功率譜估計(jì)

由仿真圖形可知，當(dāng)采樣頻率過低，fs=500 Hz<710 Hz，無法滿足采樣定理時(shí)，不能完成功率譜估計(jì)；采樣頻率逐漸增大，fs=1 000 Hz時(shí)，能夠清晰分辨出信號(hào)的2個(gè)譜峰；當(dāng)采樣頻率選擇過大，fs=5 kHz和10 kHz時(shí)，信號(hào)的2個(gè)頻率出現(xiàn)頻譜混疊，僅能發(fā)現(xiàn)一個(gè)譜峰，無法完成功率譜估計(jì)。因此，采樣頻率的選擇必須滿足采樣定理，不能過??；同時(shí)也要避免采樣頻率選擇過高導(dǎo)致的頻譜混疊，一般選擇fs∈[3fmax,4fmax]。

2.2 采樣點(diǎn)數(shù)的選擇對(duì)功率譜估計(jì)的影響

設(shè)定采樣頻率為1 kHz，AR模型階數(shù)p=100，F(xiàn)FT長(zhǎng)度為2 048，采樣點(diǎn)數(shù)分別為N=128，512，2 048和8 192，仿真結(jié)果如圖2所示。

圖2 采樣點(diǎn)數(shù)不同時(shí)的功率譜估計(jì)

由圖2可知，當(dāng)采樣點(diǎn)數(shù)很小，N=128時(shí)，譜估計(jì)的分辨率很低，而且出現(xiàn)了譜峰分裂的現(xiàn)象。隨著采樣點(diǎn)數(shù)的增加，獲得的數(shù)據(jù)會(huì)越來越多，譜估計(jì)的分辨率也越來越高，但同時(shí)也會(huì)增大整個(gè)譜估計(jì)的運(yùn)算量。對(duì)此必須綜合考慮，選擇合適的采樣點(diǎn)數(shù)，在保證分辨率的同時(shí)，盡量減小運(yùn)算量。

2.3 FFT長(zhǎng)度的選擇對(duì)功率譜估計(jì)的影響

設(shè)定采樣頻率為1 kHz，AR模型階數(shù)p=100，采樣點(diǎn)數(shù)為N=512。FFT長(zhǎng)度分別為512,1 024,2 048和4 096。仿真結(jié)果如圖3所示。

圖3 FFT長(zhǎng)度不同時(shí)的功率譜估計(jì)

由圖3可知，隨著FFT長(zhǎng)度的增大，2個(gè)譜峰變得更清晰陡峭，功率譜估計(jì)的分辨率也更高，但同時(shí)也必將導(dǎo)致算法運(yùn)算量的增大。對(duì)此，選擇FFT長(zhǎng)度時(shí)必須綜合考慮分辨率和算法運(yùn)算量，通常選擇FFT長(zhǎng)度為2 048或4 096。

2.4 AR模型階數(shù)的選擇對(duì)功率譜估計(jì)的影響

設(shè)定采樣頻率為1 kHz，F(xiàn)FT長(zhǎng)度=2 048，采樣點(diǎn)數(shù)為N=512，AR模型階數(shù)p=30，100，200和400。仿真結(jié)果如圖4所示。

圖4 AR模型階數(shù)不同時(shí)的功率譜估計(jì)

由圖4可知，當(dāng)AR模型階數(shù)較小，p=30時(shí)，分辨率很低，無法清晰分辨出2個(gè)譜峰。隨著階數(shù)的增大，分辨率逐漸增大。當(dāng)p=100和200時(shí)，能夠清晰分辨出2個(gè)譜峰；但當(dāng)階數(shù)過大，p=400時(shí)，2個(gè)譜峰處出現(xiàn)了明顯的譜峰分裂，整個(gè)譜估計(jì)出現(xiàn)了大量的虛假峰，無法正確估計(jì)出頻率所在位置。因此，在選擇AR模型階數(shù)時(shí)，必須在保證分辨率的同時(shí)避免譜峰分裂情況的出現(xiàn)。綜合比較，本文選擇p=100較為合適。

3 實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證

列車測(cè)速系統(tǒng)采用多普勒雷達(dá)體制，工作原理基于多普勒效應(yīng)。測(cè)速天線置于列車車廂下方，以角度θ向地面發(fā)射雷達(dá)波。當(dāng)列車以速度v(t)運(yùn)動(dòng)時(shí)，根據(jù)多普勒效應(yīng)，接收到的雷達(dá)波會(huì)產(chǎn)生多普勒頻移fd，如圖5所示。多普勒頻移fd與列車運(yùn)行速度v(t)的關(guān)系如下：

圖5 列車測(cè)速原理圖

(15)

式中：λ為雷達(dá)信號(hào)波長(zhǎng)；θ為波束主瓣方向與列車行駛方向的夾角；v(t)為列車相對(duì)地面的行駛速度。

為了更好地探測(cè)列車速度，避免因顛簸震動(dòng)導(dǎo)致的雷達(dá)天線與地面夾角的變化，本文采用的列車測(cè)速系統(tǒng)使用角度補(bǔ)償?shù)姆椒ǎ虻孛姘l(fā)射2個(gè)存在固定夾角的雷達(dá)波束，保證了多普勒雷達(dá)列車測(cè)速系統(tǒng)的穩(wěn)定性。

為了驗(yàn)證第2節(jié)中Burg算法參數(shù)選擇分析正確與否，本文采用多普勒雷達(dá)列車測(cè)速系統(tǒng)進(jìn)行仿真實(shí)驗(yàn)。其中雷達(dá)天線分別以與地面成40°和50°的方向向地面發(fā)射雷達(dá)波，2個(gè)雷達(dá)波束的夾角固定為10°，測(cè)速雷達(dá)的工作頻率為f0=24.125GHz。

列車測(cè)速系統(tǒng)分別采集了列車速度為10～80km/h的數(shù)據(jù)，其對(duì)應(yīng)的多普勒頻移范圍約為0.45～3.6GHz。由上節(jié)分析可知，選取采樣率fs=4fmax，取整后fs=15 GHz；選取采樣點(diǎn)數(shù)N=512；選取FFT長(zhǎng)度為2 048；AR模型階數(shù)選為p=100。對(duì)雷達(dá)天線與地面成40°角、列車速度為40km/h的數(shù)據(jù)進(jìn)行仿真，實(shí)驗(yàn)結(jié)果如圖6、圖7所示。由多普勒頻率公式可求出：

圖6 列車速度為40 km/h時(shí)數(shù)據(jù)時(shí)域圖

圖7 列車速度為40 km/h時(shí)功率譜估計(jì)

對(duì)比圖6可知，按照本文第2節(jié)中分析結(jié)果選擇Burg算法仿真參數(shù)，成功地完成了對(duì)列車測(cè)速數(shù)據(jù)的功率譜估計(jì)，并且估計(jì)結(jié)果較為準(zhǔn)確。

4 結(jié)束語

本文首先通過MATLAB仿真信號(hào)分析了采樣頻率、采樣點(diǎn)數(shù)、AR模型階數(shù)和FFT長(zhǎng)度等Burg算法參數(shù)對(duì)功率譜估計(jì)的影響。由分析可知，采樣頻率的選擇必須滿足采樣定理，同時(shí)也要避免采樣頻率選擇過高導(dǎo)致的頻譜混疊，一般選擇fs∈[3fmax,4fmax]；選擇合適的采樣點(diǎn)數(shù)，在保證分辨率的同時(shí)，盡量減小運(yùn)算量；選擇FFT長(zhǎng)度時(shí)必須綜合考慮分辨率和算法運(yùn)算量；在選擇AR模型階數(shù)時(shí)，必須在保證分辨率的同時(shí)避免譜峰分裂情況的出現(xiàn)。

為了驗(yàn)證分析結(jié)果，本文采用多普勒雷達(dá)列車測(cè)速系統(tǒng)進(jìn)行仿真實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，按照本文分析選擇Burg算法仿真參數(shù)，能夠成功地完成對(duì)列車測(cè)速數(shù)據(jù)的功率譜估計(jì)，并得到良好的功率譜估計(jì)效果。

[1] 王鳳瑛,張麗麗.功率譜估計(jì)及其MATLAB仿真[J].仿真技術(shù),2006,22(11):287-289.

[2] 羅豐,段沛沛,吳順君.基于Burg算法的短序列譜估計(jì)研究[J].西安電子科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2005,32(5):724-728.

[3] 段沛沛,羅豐,吳順君.一種有效的短序列功率譜估計(jì)算法及其應(yīng)用[J].雷達(dá)科學(xué)與技術(shù),2004,2(6):272- 275.

[4] 蔡方凱.已調(diào)波功率譜分析研究[D].成都:電子科技大學(xué),2005.

[5] 沈慧芳,呂容川,蔡立娜,等.基于修正Levinson算法的地雜波仿真及其應(yīng)用[J].火控雷達(dá)技術(shù),2007,36(2):17-19.

[6] 姚文俊.自相關(guān)法和Burg法在AR模型功率譜估計(jì)中的仿真研究[J].計(jì)算機(jī)與數(shù)字工程,2006,35(10):32- 34.

[7] 田坦,李延.伯格譜估計(jì)算法的一種改進(jìn)[J].數(shù)據(jù)采集與處理,2002,17(3):276-278.

[8]MARPLESL.Anewautoregressivespectrumanalysisalgorithm[J].IEEETransactionsonAcoustics,Speech,andSignalProcession,1980,28(4):441-454.

ParametersSelectionofPowerSpectrumEstimationBasedonBurgAlgorithm

HUANGChao,WANGWei,SHANLiang

(ElectronicEngineeringInstitute,Hefei230037,China)

ThispaperintroducesBurgalgorithmbasedonautoregressive(AR)modelinmodernpowerspectrumestimation,analyzestheinfluenceofparameterselectioninBurgalgorithmonpowerspectrumestimation,andgivesthecomprehensiveoptimalparameterselectionmethodbasedonBurgalgorithm,finallyappliestheparameterselectionmethodproposedinthispapertotrainspeeddetectionsystemforexperimentvalidation,estimatestheDopplerfrequencyoftrainaccurately.

powerspectrumestimation;autoregressivemodel;Burgalgorithm;parameterselection

2015-08-31

TN

A

CN32-1413(2016)05-0085-05

10.16426/j.cnki.jcdzdk.2016.05.022