AR模型譜估計實用教案

1、平穩(wěn)隨機信號的參數(cnsh)模型qkkpkkknubknxanx01)()()(0)()()(kknukhnx)()()(zAzBzHpkkkzazA11)(qkkkzbzB11)(0)()(kkzkhzH)(jwxeP222)()(jwjweAeB212221)(pkjwkkjweaeA0kb為白噪聲(zoshng) 的方差2)(kuAR模型全極點：該模型現(xiàn)在的輸出(shch)是現(xiàn)在的輸入和過去p個輸出(shch)的加權和。第2頁/共21頁第1頁/共21頁第一頁，共22頁。 經常使用的參數(cnsh)模型是線性模型, 其中以有理分式 型應用最為普遍, 有理分式模型一般表現(xiàn)為自回歸動均 模型

2、, 即ARMA模型, 自回歸模型( AR模型) 和動均模 型( MA模型) 都是YuleWalker 模型的特殊情況。事 實上, 白噪聲序列通過全極點型、全零點型濾波器會分 別產生AR, MA和ARMA過程。在這三種參數(cnsh)模型 中, AR模型得到了普遍應用, 因為AR模型的參數(cnsh)計 算是線性方程, 比較簡單, 與建立在外推自相關函數時 保持原概率空間的最大熵法是等價的, 同時很適合表示 很窄的頻譜, 在作譜估計時, 由于具有遞推特性所以所 需的數據較短; 而MA模型表示窄譜時一般需要數量很多的參數(cnsh); ARMA模型雖然所需的參數(cnsh)數量最少, 但 參數(c

3、nsh)估計的算法是非線性方程組, 其運算遠比AR模型 復雜, 故AR模型參數(cnsh)估計是重點。第3頁/共21頁第2頁/共21頁第二頁，共22頁。 AR模型法 AR模型法的基本理論 任何具有功率譜密度的隨機信號都可以看成由白 噪聲 激勵一物理網絡所形成, 可寫成: 該形式稱為(chn wi)p階自回歸模型, 簡稱AR模型。將其進行z 變換可得AR模型的傳遞函數為: 自回歸模型的H(z) 只有極點, 沒有除原點以外的零點, 因此又稱為(chn wi)全極點型。當用自回歸模型時, 功率譜密度的表達式寫成:第4頁/共21頁第3頁/共21頁第三頁，共22頁。 式中: 為 白噪聲的功率譜密度。因此

4、只要求解出 及所有 的值, 就可以得到隨機(su j)信號x(n) 的功率譜。目前估計AR模型參數的方法有: ( 1) 相關函數類算法: ( 2) 反射系數類算法: ( 3) 最小二乘類算法: 第5頁/共21頁第4頁/共21頁第四頁，共22頁。 ( 1) 相關函數類算法: 先估計自相關序列, 然后解YuleWalker方程算出AR系數, 計算出功率譜。經 常使用有偏自相關估計, 以保證自相關矩陣正定性, 此 法適于較長的序列。 ( 2) 反射系數類算法: 它不需要(xyo)估計過程的自相關函數, 而是按照Levinson遞推公式直接從序列值遞推計算預測誤差和反射系數( 遞推估計反射系數時,是使

5、各階的平均預測誤差功率最小) , 最后得到模型系數, 計算出功率譜。優(yōu)點是分辨率高, 適于短序列。具體的算法有Burg 算法和Itakura算法, 兩者的主要區(qū)別在于計算平均預測誤差功率的方法不同, 前者采用向前、向后預測誤差功率的算術平均而后者取幾何平均。 ( 3) 最小二乘類算法: 可有多種算法, 例如直接用LS方法擬合出模型參數或是采用FTF算法求出模 型參數。第6頁/共21頁第5頁/共21頁第五頁，共22頁。 AR模型法的仿真 YuleWalker 方程 采樣頻率為200Hz, 采樣點數為50和200時, 采 用YuleWalker 方程的仿真圖。 當采樣點數為50時, 在頻率為40H

6、z處的譜峰明顯向左偏移, 且在頻率為60 70Hz之間出現(xiàn)了虛假譜峰。當采樣點數為200時, 譜峰回到40Hz處, 且虛假譜峰的幅度變小。因此(ync)YuleWalker 方程的性能可通過增加采樣點數來解決。第7頁/共21頁第6頁/共21頁第六頁，共22頁。 AR模型的參數和x(n)自相關函數有如下的關系: 將上式寫成矩陣的形式(xngsh): 即是AR模型的正則方程, 又稱尤拉沃克 (YuleWalker)方程。第8頁/共21頁第7頁/共21頁第七頁，共22頁。AR模型(mxng)參數估計的典型算法 1.自相關法 自相關法是AR模型參數求解中最簡單的一種方法。L-D遞推算法是在滿足前向預測

7、均方誤差最小的前提下,先求得觀測數據的自相關函數,然后利用YuleWalker 方程的遞推性質求得模型參數,進而求得功率譜的估值。它是模型階次逐次(zh c)加大的一種算法,即先計算階次m=1時的預測系數 ，再計算m=2時的a2(1) , a2(2)和2,按此依次計算到階次m=p時的ap(1) , ap(2) , , ap(p)及2p ,當2p 滿足精度要求時即可停止遞推。第9頁/共21頁第8頁/共21頁第八頁，共22頁。 遞推公式為: Burg算法 用Burg算法進行功率譜估計時令前后向預測誤差功率之和最小,即對前向序列誤差和后向序列誤差前后都不加(b ji)窗,使用LevinsonDurb

8、in遞推可快速的求解AR系數。Burg算法與自相關法不同,它是使序列x(n)的前后向預測誤差功率之和:第10頁/共21頁第9頁/共21頁第九頁，共22頁。 最小。在上式中,當階次由1至p時, 和 （n）有以下的遞推關系(gun x): 可知 僅為km 的函數。令 ,可得到:第11頁/共21頁第10頁/共21頁第十頁，共22頁。 再利用Levinson2Durbin遞推算法可得AR模型(mxng)系數: Burg算法是建立在數據基礎之上的,避免了先計算自相關函數從而提高計算速度;是較為通用的方法,計算不太復雜,且分辨率優(yōu)于自相關法,但對于白噪聲加正弦信號有時會出現(xiàn)譜線分裂現(xiàn)象。第12頁/共21頁

9、第11頁/共21頁第十一頁，共22頁。改進(gijn)協(xié)方差算法 同Burg算法一樣,改進協(xié)方差（修正協(xié)方差）算法進行功率譜估計時令前后向預測誤差功率之和最小,即對前后向預測誤差都不加(b ji)窗,但得到的協(xié)方差矩陣不是Toeplitz矩陣,因此正則方程不能用Levinson遞推算法求解。Marple于1980年提出了實現(xiàn)協(xié)方差方程求解的快速算法,大大提高了譜估計的性能。第13頁/共21頁第12頁/共21頁第十二頁，共22頁。 參數提取時,利用Matlab工具箱中的信號處理中的Levinson函數和arburg函數來分別進行自相關算法(sun f)和Burg算法(sun f)的AR模型參數估

10、計,以AR模型功率譜估計及Matlab實現(xiàn) 平穩(wěn)隨機信號第14頁/共21頁第13頁/共21頁第十三頁，共22頁。 clc; clear;close all; fs=1000;% 采樣率 1000 N=1;% 改變數據長度 p=50;%AR 模型階數 nfft=512;%fft長度 t=0:1/ fs: N; wn=sqrt(1)+randn(1,N*fs+1); %白噪聲,均值0, 方差1 s1=sqrt(20)*sin(2*pi*100*t); % 正弦信號1, 信噪比10db s2=sqrt(2000)*sin(2*pi*110*t); % 正弦信號 2, 信噪比30db x=s1+s2+

11、wn; % 觀測(gunc)數據 %figure,plot(t,x); x1=xcorr(x,biased); Pxx,f=pyulear(x1,p,nfft,fs); %Yule -Walker 方程 figure,plot(f,10*log10(Pxx); grid on; title(自相關法); 第15頁/共21頁第14頁/共21頁第十四頁，共22頁。 Pxx1,f1=pcov(x,p,nfft,fs); figure,plot(f1,10*log10(Pxx1); grid on; title(協(xié)方差法); Pxx2,f2=pmcov(x,p,nfft,fs); figure,plo

12、t(f2,10*log10(Pxx2); grid on; title(修正(xizhng)協(xié)方差法); Pxx3,f3=pburg(x,p,nff t,fs); figure,plot(f3,10*log10(Pxx3); grid on; title(Burg 法); 第16頁/共21頁第15頁/共21頁第十五頁，共22頁。從上圖可以看出，采用參數建模的譜估計方法得到的功率譜曲線平滑（方差?。?，分辨率高，可以明顯地觀察(gunch)到兩個譜峰。第17頁/共21頁第16頁/共21頁第十六頁，共22頁。第18頁/共21頁第17頁/共21頁第十七頁，共22頁。 降低( jingd)模型階次后，可以發(fā)現(xiàn)，譜的分辨率降低( jingd)（兩個譜峰幾乎變成一個譜峰），但是曲線平滑性變好（估計誤差降低( jingd)第19頁/共21頁第18頁/共21頁第十八頁，共22頁。第20頁/共21頁第19頁/共21頁第十九頁，共22頁。thanks 通過對功率譜的仿真比較可以得到:自相關法的計算簡單( jindn),但譜估計的分辨率較差,而Burg算法和改進協(xié)方差算法是較為通用的方法,計算不太復雜,且具有較好的譜估計質量。第21頁/共21頁第20頁/共21頁第二十頁，共