劉明曉，王旭光

（1.西安鐵路職業(yè)技術(shù)學(xué)院 陜西 西安 710014；2.西安航天自動(dòng)化股份有限公司 陜西 西安 710014）

基于MATLAB實(shí)現(xiàn)的AR模型功率譜估計(jì)

劉明曉1，王旭光2

（1.西安鐵路職業(yè)技術(shù)學(xué)院 陜西 西安 710014；2.西安航天自動(dòng)化股份有限公司 陜西 西安 710014）

現(xiàn)代功率譜估計(jì)和經(jīng)典功率譜估計(jì)是兩種常用的功率譜估計(jì)，同時(shí)也是分析隨機(jī)信號的常用方法。本文詳細(xì)介紹了現(xiàn)代功率譜估計(jì)中有關(guān)AR模型參數(shù)的功率譜估計(jì)，具體包括自相關(guān)算法、Burg算法、協(xié)方差算法以及改進(jìn)的協(xié)方差算法，在此基礎(chǔ)上又對四種不同功率譜估計(jì)方法的性能指標(biāo)進(jìn)行了比較分析。在MATLAB仿真軟件平臺上對AR模型參數(shù)的四種不同功率譜估計(jì)算法進(jìn)行了仿真，同時(shí)對功率譜估計(jì)結(jié)果進(jìn)行了分析比較并得到了預(yù)期的譜估計(jì)效果。最后從實(shí)際應(yīng)用角度出發(fā)討論了AR模型參數(shù)不同功率譜估計(jì)算法的特點(diǎn)，以便在應(yīng)用中能夠選擇合適的功率譜估計(jì)算法。

功率譜估計(jì)；AR模型；算法；MATLAB

在頻譜分析中功率譜估計(jì)是常被采用的一種重要方法，功率譜估計(jì)廣泛應(yīng)用于人們的日常生活、聲納、生物醫(yī)學(xué)、雷達(dá)以及農(nóng)業(yè)活動(dòng)等不同領(lǐng)域的信號處理過程中[1-2]。如何從頻帶較寬的信號中檢測出頻帶較窄的信號也是功率譜估計(jì)在信號處理中的重要應(yīng)用[3]。此外功率譜估計(jì)是用有限長的數(shù)據(jù)來對信號進(jìn)行功率譜估計(jì)，是分析隨機(jī)信號時(shí)常被采用的一種重要方法。

為了能夠良好地解決經(jīng)典功率譜估計(jì) （包括自相關(guān)法和周期圖法）中方差性能不好以及分辨率較低等問題提出了現(xiàn)代功率譜估計(jì)[4]，其涉及的領(lǐng)域以及學(xué)科相當(dāng)廣泛且內(nèi)容比較豐富[5]?，F(xiàn)代功率譜估計(jì)包括非參數(shù)模型估計(jì)和參數(shù)模型估計(jì)兩種，前者有最小方差法，多分量的MUSIC方法等；后者有AR模型、MA模型、ARMA模型、PRONY指數(shù)模型等[6-7]。文中對現(xiàn)代功率譜估計(jì)中有關(guān)AR模型參數(shù)的幾種典型求解算法進(jìn)行分析比較，同時(shí)在MATLAB仿真軟件平臺上對AR模型參數(shù)中不同算法的功率譜估計(jì)進(jìn)行了仿真，最后從實(shí)際應(yīng)用的角度出發(fā)討論了不同功率估計(jì)算法的特點(diǎn)，以便在實(shí)際應(yīng)用中能夠選擇最合適的功率譜估計(jì)算法[8]。

1 現(xiàn)代功率譜估計(jì)

現(xiàn)代功率譜估計(jì)又稱參數(shù)譜估計(jì)[9]，是先對數(shù)據(jù)進(jìn)行觀察后再進(jìn)一步對參數(shù)模型進(jìn)行估計(jì)，然后根據(jù)得到參數(shù)模型輸出的功率來對信號進(jìn)行譜估計(jì)。采用現(xiàn)代功率譜估計(jì)能夠較好地改善經(jīng)典譜估計(jì)中方差性能差以及分辨率低等問題[10]?，F(xiàn)代功率譜估計(jì)中AR模型是應(yīng)用較多的一種模型。

1.1 AR模型的建立

AR模型是全極點(diǎn)模型同時(shí)又稱為自回歸模型[11]?？捎霉剑?）所示的差分方程來描述：

其中u（n）是一個(gè)白噪聲序列，其均值為0、方差為 σ2，p 是 AR 模型的階數(shù)，ap（i），i=1，2，…，p。AR 模型系統(tǒng)的轉(zhuǎn)移函數(shù)H（z）為：

由公式（1）和公式（2）可得到AR模型參數(shù)的功率譜估計(jì)計(jì)算公式：

AR模型的參數(shù)和x（n）自相關(guān)函數(shù)有如下的關(guān)系：

將上式寫成矩陣的形式：

式（5）為AR模型的正則方程，又稱尤拉-沃克（Yule-Walker）方程[12]。

1.2 AR模型參數(shù)估計(jì)的典型算法

1.2.1 自相關(guān)法

自相關(guān)法是AR模型參數(shù)估計(jì)中較為簡單的一種功率譜估計(jì)方法。按照模型階數(shù)由小到大的順序進(jìn)行計(jì)算，先計(jì)算階次 m=1 時(shí)的系數(shù) am（k）=a1（1）和，再計(jì)算 m=2 時(shí)的 a2（1），a2（2）和，按此方法依次計(jì)算到階次 m=p 時(shí)的 ap（1），ap（2），…，ap（p），及，當(dāng)?shù)木葷M足要求時(shí)即可停止。遞推公式為：

1.2.2 Burg算法

Burg算法進(jìn)行功率譜估計(jì)時(shí)其前后項(xiàng)序列誤差的前后都不加窗，使用Levinson-Durbin遞推可快速求解AR模型參數(shù)的系數(shù)。

Burg算法與自相關(guān)算法不同之處在于使序列x（n）的前后項(xiàng)功率預(yù)測誤差之和：

再根據(jù)Levinson-Durbin進(jìn)行遞推得出AR模型系數(shù)如下：

Burg算法不需先計(jì)算自相關(guān)函數(shù)進(jìn)而計(jì)算速度有所提高同時(shí)計(jì)算也相對簡單。與自相關(guān)算法相比分辨率也有所提高。但在處理白噪聲加正弦信號時(shí)功率譜線會(huì)出現(xiàn)分裂現(xiàn)象。

1.2.3 協(xié)方差算法

協(xié)方差算法與其他算法的不同之處在于預(yù)測誤差功率求和的上下限取值范圍不同。在協(xié)方差算法區(qū)間[0，N-1]外 x（n）樣本不為 0，故表達(dá)式中的 x（n-k）總是落在區(qū)間[0，N-1]中，因此預(yù)測誤差功率求和的上下限必然會(huì)落在[p，N-1]之間。得到的矩陣是半正定矩陣，且不具有Toeplitz性質(zhì)，故不能采用Levinson-Durbin遞歸算法進(jìn)行參數(shù)求解，同時(shí)得到的AR模型也不穩(wěn)定。

1.2.4 改進(jìn)協(xié)方差算法

采用改進(jìn)的協(xié)方差算法進(jìn)行功率譜估計(jì)時(shí)，為了使前后項(xiàng)功率誤差之和最小其前后項(xiàng)都不加窗，得到的協(xié)方差矩陣不是Toeplitz矩陣，故不能用Levinson遞推算法對模型參數(shù)進(jìn)行求解。Marple在1980年提出了改進(jìn)協(xié)方差方程求解的快速算法，該算法能夠大大提高功率譜估計(jì)的性能。

1.3 用于功率譜估計(jì)的數(shù)據(jù)

文中采用由復(fù)數(shù)噪聲加上4個(gè)復(fù)正弦組成的數(shù)據(jù)對功率譜進(jìn)行估計(jì)。其中歸一化頻率分別是f1=0．15，f2=0．16，f3=0．252，f4=-0．16。 該數(shù)據(jù)的 128 點(diǎn)復(fù)序列如圖1所示。

圖1 功率譜估計(jì)數(shù)據(jù)

2 AR模型功率譜估計(jì)的MATLAB仿真

借助MATLAB對AR模型功率譜進(jìn)行估計(jì)時(shí)，上述4種算法分別采用pyulear函數(shù)、pburg函數(shù)、pcov函數(shù)和pmcov函數(shù)來完成功率譜估計(jì)。分別用這4個(gè)函數(shù)對圖1中的隨機(jī)序列進(jìn)行功率譜估計(jì)便可得到相應(yīng)的功率譜線。圖2～4是自相關(guān)算法得出的階次p分別為10，20和30時(shí)的AR模型功率譜曲線。

圖2 階次p=10時(shí)自相關(guān)算法功率譜估計(jì)

從圖2～4的功率譜估計(jì)圖形可以得出，在階次p較低時(shí) （如圖2所示）功率譜估計(jì)的分辨率不理想。當(dāng)階次p提高到30時(shí)，在頻率f1和f2處兩個(gè)正弦信號剛好被分離開，在頻率f3和f4處也能夠分離出兩個(gè)正弦信號。在自相關(guān)算法中由于對序列的前后項(xiàng)加窗降低了分離效果，數(shù)據(jù)越短分離效果越不好。

圖3 階次p=20時(shí)自相關(guān)算法功率譜估計(jì)

圖4 階次p=30時(shí)自相關(guān)算法功率譜估計(jì)

圖5～圖6是Burg算法對隨機(jī)數(shù)列信號進(jìn)行功率譜估計(jì)時(shí)的功率譜曲線。

圖5 階次p=10時(shí)Burg算法功率譜估計(jì)

圖6 階次p=15時(shí)Burg算法功率譜估計(jì)

比較分析圖5和圖6可以得出在階次p=15時(shí)得到了較好的譜估計(jì)。另外對比圖4和圖5可以看到Burg算法的分辨率明顯優(yōu)于自相關(guān)算法，在階次較低時(shí)（如圖5所示）也能較好地分辨出間隔較小的頻率點(diǎn)。

圖7為階數(shù)p=10時(shí)協(xié)方差算法對隨機(jī)數(shù)列信號進(jìn)行估計(jì)得到的功率譜曲線。

圖7 階次p=10時(shí)協(xié)方差算法功率譜估計(jì)

從圖7可得出功率譜估計(jì)在信號源頻率f1=0．15，f2=0．16，f3=0．252，f4=-0．16 處功率譜線相對較窄，在其他頻率處功率譜線起伏劇烈。采用協(xié)方差算法對功率譜進(jìn)行估計(jì)能夠較真實(shí)地反映信號實(shí)際模型。

圖8～圖9是階數(shù)p=10和p=15時(shí)改進(jìn)的協(xié)方差算法得出的AR模型功率譜估計(jì)曲線。

圖8 階次p=10時(shí)改進(jìn)協(xié)方差算法功率譜估計(jì)

圖9 階次p=15時(shí)改進(jìn)協(xié)方差算法功率譜估計(jì)

比較圖7～圖9可以得出兩種協(xié)方差算法的功率譜估計(jì)大致相同，但改進(jìn)的協(xié)方差算法在信號頻率處的峰值更加尖銳、更加突出且更有利于辨識。分析圖8和圖9可以得出在階次較高時(shí)能得到非常滿意的功率譜估計(jì)曲線。

3 結(jié) 論

文中介紹了現(xiàn)代功率譜估計(jì)中AR模型參數(shù)的幾種功率譜估計(jì)算法。其中自相關(guān)算法需要對前項(xiàng)預(yù)測的誤差序列前后加窗，加窗降低了自相關(guān)算法的分辨率，數(shù)據(jù)越短分辨率越低。與自相關(guān)算法相比Burg算法有較好的分辨率，但有時(shí)會(huì)出現(xiàn)譜估計(jì)曲線分裂現(xiàn)象。改進(jìn)的協(xié)方差算法性能最好，但計(jì)算與編程都相對繁瑣。在對功率譜進(jìn)行估計(jì)時(shí)，需要依據(jù)實(shí)際數(shù)據(jù)序列選擇更合適的算法進(jìn)行功率譜估計(jì)。

[1]鄧澤懷，劉波波，李彥良．常見的功率譜估計(jì)方法及其Matlab仿真[J]．電子科技，2014（2）:50-52．

[2]李軍，王凱．三種現(xiàn)代功率譜估計(jì)方法性能研究[J]．科技信息，2010（30）:554-555．

[3]李文超．現(xiàn)代功率譜估計(jì)在設(shè)備故障診斷應(yīng)用[J]．傳感器世界，2012（2）:25-29．

[4]王春興．基于Matlab實(shí)現(xiàn)現(xiàn)代功率譜估計(jì)[J]．現(xiàn)代電子技術(shù)，2011，34（16）:65-67．

[5]范澤華，白鐵成．基于改進(jìn)差分演化的高效圖像分割算法[J]．現(xiàn)代電子技術(shù)，2016，39（14）:87-91．

[6]楊蓉．多普勒雷達(dá)測速系統(tǒng)設(shè)計(jì)及信號處理方法研究[D]．武漢：華中科技大學(xué)，2012．

[7]汪勝波．功率譜估計(jì)在寬帶ADCP信號檢測中的研究與應(yīng)用[D]．青島：中國海洋大學(xué)，2011．

[8]閆慶華，程兆剛，段云龍．AR模型功率譜估計(jì)及Matlab實(shí)現(xiàn) [J]．計(jì)算機(jī)與數(shù)字工程，2010，38（4）:154-156．

[9]郭會(huì)娟．基于AR參數(shù)模型和小波分析的電機(jī)故障診斷[D]．鄭州：華北水利水電大學(xué)，2013．

[10]崔曉榮．電力系統(tǒng)間諧波的檢測算法研究與實(shí)現(xiàn)[D]．成都：西華大學(xué)，2013．

[11]張峰，伊麗，石現(xiàn)峰．振動(dòng)信號的最小方差譜估計(jì)算法[J]．計(jì)算機(jī)技術(shù)與發(fā)展，2015（10）:76-79．

[12]袁華，趙越，藺永政，等．基于最小方差的功率譜估計(jì)方法的研究[J]．?dāng)?shù)字通信世界，2015（9）:254-255．

The realization of AR model power spectrum estimation based on MATLAB

LIU Ming-xiao1，WANG Xu-guang2
（1.Xi'an Railway Vocational&Technical Institute，Xi'an 710014，China； 2.Xi'an Aerospace Automation CO，.LTD.，Xi'an 710014，China）

Power spectrum estimation can be divided into modern spectral estimation and classical spectral estimation，it is an method for analyzing random signal.It describes the autocorrelation algorithm，Burg algorithm，covariance algorithm and improved covariance algorithm of AR model parameters in Modern Spectral Estimation，and their perform indicators are analyzed in this paper.These methods for AR model parameter algorithm are implemented by MATLAB，the results of simulation are analyzed and get a better spectrum estimation.Moreover，the advantages and disadvantages of these methods are discussed from the experimental view so as to make reasonable selection in practical work.

PSD estimation； AR model； algorithm；MATLAB

TN015

A

1674－6236（2017）17-0129-04

2016-07-23稿件編號：201607167

劉明曉（1987—），女，河南焦作人，碩士，助教。研究方向：電力系統(tǒng)及其自動(dòng)化。