AR模型和ARMA模型譜估計(jì)仿真

1、AR模型和ARMA模型譜估計(jì)仿真一、問題重述有兩個(gè)ARMA過程，其中信號(hào)1是寬帶信號(hào)，信號(hào)2是窄帶信號(hào)，分別用AR譜估計(jì)算法、ARMA譜估計(jì)算法和周期圖法估計(jì)其功率譜。產(chǎn)生信號(hào)1的系統(tǒng)函數(shù)為：Hz=1+0.3544z-1+0.3508z-2+0.1736z-3+0.2401z-41-1.3817z-1+1.5632z-2-0.8843z-3+0.4906z-4激勵(lì)白噪聲的方差為1.產(chǎn)生信號(hào)2的系統(tǒng)函數(shù)為：Hz=1+1.5857z-1+0.9604z-21-1.6408z-1+2.2044z-2-1.4808z-3+0.8145z-4激勵(lì)白噪聲的方差為1.每次實(shí)驗(yàn)使用的數(shù)據(jù)長度為256.二、模型

2、分析+很多隨機(jī)過程可以由或近似由均值為零、方差為2的白噪聲序列un經(jīng)過具有有理想傳輸函數(shù)H(z)的ARMA線性系統(tǒng)來得到。稱該隨機(jī)過程為ARMA過程。U(z)X(z)B(z)A(z)-1Hz=i=0qbiz-ii=0paiz-i=B(z)A(z)Pxxw=2H(w)2由上式可知只要估計(jì)出模型的參數(shù)（ai和bi），即可求出功率譜。1.AR模型的建立：AR模型是一種特殊的ARMA模型，利用AR（p）模型，即：xn=-i=1paixn-i+u(n)逼近采樣樣本,此時(shí)功率譜表達(dá)式為：Px=21+i=1Paie-jwi2需要求解得未知量為參數(shù)ai，當(dāng)階數(shù)p已知時(shí)，利用x(n)的自相關(guān)函數(shù)與AR模型參數(shù)的

3、關(guān)系，可建立Y-W方程，解該方程，即可得到AR參數(shù)。Rxm+i=1paiRxm-i=0 m=1,2pRxm+i=1paiRxm-i=2 m=0R(0)R(-1)R(1)R(0)R(-p)R(1-p)R(p)R(p-1)R(0)1a1ap=200對于(p+1)元線性方程，若采用matlab中的函數(shù)，則使用的是高斯消去法運(yùn)算量為p3數(shù)量級。采用下面的Levinson-Durbin算法。它可以將運(yùn)算量減少到p2數(shù)量級。遞推公式的獲取方法如下：令p=1，得一階AR模型對應(yīng)的尤勒沃克方程為Rx0+a11Rx1=1Rx1+a11Rx0=0可解得a11=-Rx1Rx01=Rx0-Rx21Rx0=01-a12

4、(1)令p=2,得二階AR模型對應(yīng)的尤勒沃克方程為Rx0+a21Rx1+a22Rx2=2Rx1+a21Rx0+a22Rx1=0Rx2+a21Rx1+a22Rx0=0解此方程組，得a22=-Rx2+a1(1)Rx(1)/1a21=a11+a22a11 2=11-a222 令p=3,4,由遞推規(guī)律可得到下式：kp=-Rxp+i=1p-1ap-1(i)Rx(p-i)/p-1 api=ap-1i+appap-1p-i, i=1,2,p-1p=p-11-kp2 式中，kp=app, m=1,2,p 0=Rx(0).按照此遞推式即可方便解出Ar模型參數(shù)a。2.ARMA模型ARMA(p,q)的差分方程如下：

5、xn=-i=1paixn-i+i=0qbiu(n-i)由此可得出自相關(guān)的關(guān)系：i=0paiRxm-i=2k=0h*kbk+mRm=-i=1paiRm-i+2i=0q-mh*ibk+m m=0,1,2,q-i=1paiRm-i+2i=-mq-mh*ibi+m m<0 -i=1paiRm-i m>q 當(dāng)m>q時(shí)的自相關(guān)函數(shù)關(guān)系可得,建立超定方程（M>p+q）：R(q)R(q-1)R(q+1)R(q)R(q-p-1)R(q-p+2)R(M-1)R(M-p)R(q)a1ap=-R(q+1)-R(M)Ra=r利用最小二乘解得： a=(RHR)-1RHr利用Kaveh譜估計(jì)方法，在

6、計(jì)算出a后直接計(jì)算參數(shù)ck，則此時(shí)功率譜估計(jì)式變?yōu)椋篜xz=k=-qqckz-kA(z)A(z-1)而參數(shù)ck可以由以下關(guān)系式求出：ck=i=0pj=0paiai*Rxk-i+j k=0,1,q三、實(shí)驗(yàn)內(nèi)容信號(hào)一：按照題目要求要對信號(hào)一分別使用AR(4),AR(8),ARMA(4,4),ARMA(8,8)模型和周期圖法進(jìn)行估計(jì)，做20次獨(dú)立實(shí)驗(yàn)結(jié)果繪制在一張圖上，并繪制20次的平均譜和真實(shí)譜。信號(hào)一是均值為0，方差為1的白噪聲通過如下系統(tǒng)函數(shù)產(chǎn)生的。得到如下實(shí)驗(yàn)結(jié)果：信號(hào)二：按照題目要求要對信號(hào)二分別使用AR(4),AR(8),AR(12),AR(16),ARMA(4,2),ARMA(8,4)

7、, ARMA(12,6)模型和周期圖法進(jìn)行估計(jì)，做20次獨(dú)立實(shí)驗(yàn)結(jié)果繪制在一張圖上，并繪制20次的平均譜和真實(shí)譜。信號(hào)一是均值為0，方差為1的白噪聲通過如下系統(tǒng)函數(shù)產(chǎn)生的。實(shí)驗(yàn)結(jié)果如下圖所示：四、結(jié)果分析為比較算法的不同，現(xiàn)對兩種不同的信道做以分析，信道零極點(diǎn)圖如下：由此可知信道二比信道一更不穩(wěn)定。實(shí)驗(yàn)結(jié)果展示出采用不同方法有著不同的特點(diǎn)，周期圖法做出的譜估計(jì)圖最差，毛刺較多，對穩(wěn)定性較差的系統(tǒng)估計(jì)不是很好。AR模型有較好的性能，估計(jì)譜比較平滑，但對穩(wěn)定性較差的系統(tǒng)的估計(jì)也不是很好。ARMA模型估計(jì)譜也比較平滑，其性能的最主要優(yōu)點(diǎn)體現(xiàn)在其對穩(wěn)定性較差的系統(tǒng)的估計(jì)能力比較強(qiáng)，但是在每次估計(jì)時(shí)會(huì)出

8、現(xiàn)某些極差的數(shù)據(jù)點(diǎn)，其算法穩(wěn)定性較差。同時(shí)通過實(shí)驗(yàn)可以看出，AR和ARMA的性能并不是階數(shù)越大越好，所以應(yīng)該采用一定的方法確定其階數(shù)，能得到較好效果。附 件AR模型法clear;clc;echo off;close all;h1=1 0.3544 0.3508 0.1736 0.2401 1 -1.3817 1.5632 -0.8843 0.4906;h2=1 1.5857 0.9604 1 -1.6408 2.2044 -1.4808 0.8145;sigma=1;h=h2;% 繪制信道功率譜H,w=freqz(cell2mat(h(1),cell2mat(h(2),400);Htol=ze

9、ros(400,20);% 產(chǎn)生500點(diǎn)數(shù)據(jù)for k=1:20, %做20次譜估計(jì)實(shí)驗(yàn)N=500;ha=cell2mat(h(1);hb=cell2mat(h(2);Na=length(ha);%噪聲項(xiàng)系數(shù)長度Nb=length(hb);%信號(hào)項(xiàng)系數(shù)長度x=zeros(1,N+Nb-1);e=wgn(1,N+Na-1,sigma,'linear','real');for i=1:1:N, x(i+Nb-1)=fliplr(ha)*e(i:i+Na-1)'-fliplr(hb(2:Nb)*x(i:i+Nb-2)'endx=x(5:end);% 奇

10、異值分解法確定模型階數(shù)% for k=1:20,temp=xcorr(x,'biased');cor_l=round(length(temp)/2);Cxx=zeros(cor_l,cor_l);for i=1:cor_lCxx(i,:)=temp(cor_l+1-i:length(temp)+1-i);endX,B=eig(Cxx);V=diag(B)./max(diag(B);%歸一化特征值% figure;% plot(1:N,V);V(V<0.5)=0;P=0;for i=1:N, if V(i)>0 P=P+1; endend% 計(jì)算AR模型參數(shù)P=16;

11、Rxx=temp(cor_l:cor_l+P);Perr=zeros(1,P+1);Perr(1)=Rxx(1);G=ones(P+1,P+1);G(2,2)=-(Rxx(2)/Rxx(1);Perr(2)=Perr(1)*(1-G(2,2)2);for i=2:P,% for j=1:i, G(i+1,i+1)=-(fliplr(Rxx(2:i+1)*G(i,1:i)')/Perr(i); G(i+1,2:i)=G(i,2:i)+G(i+1,i+1)*fliplr(G(i,2:i); Perr(i+1)=Perr(i)*(1-G(i+1,i+1)2);endHa,wa=freqz(1

12、,G(P+1,:),400); Hf=(abs(Ha).*sigma).2; Htol(:,k)=10*log10(Hf); figure(1); plot(wa,10*log10(Hf),'g');hold on;endsum=0;for j=1:20; sum=sum+Htol(:,j);endfigure(1);plot(wa,sum/20,'r');hold on;H,w=freqz(cell2mat(h(1),cell2mat(h(2),400); Hf=(abs(H).*sigma).2; figure(1); plot(w,10*log10(Hf)

13、; Je=num2str(P); str1=strcat('離散系統(tǒng)功率譜 AR(',Je,')');title(str1);grid on;ARMA模型法clear;clc;echo off;close all;h1=1 0.3544 0.3508 0.1736 0.2401 1 -1.3817 1.5632 -0.8843 0.4906;h2=1 1.5857 0.9604 1 -1.6408 2.2044 -1.4808 0.8145;sigma=1;h=h2;% 繪制信道功率譜H,w=freqz(cell2mat(h(1),cell2mat(h(2),4

14、00);Htol=zeros(400,20);% 產(chǎn)生500點(diǎn)數(shù)據(jù) for k=1:20, %做20次譜估計(jì)實(shí)驗(yàn)N=500;ha=cell2mat(h(1);hb=cell2mat(h(2);Na=length(ha);%噪聲項(xiàng)系數(shù)長度Nb=length(hb);%信號(hào)項(xiàng)系數(shù)長度x=zeros(1,N+Nb-1);e=wgn(1,N+Na-1,sigma,'linear','real');for i=1:1:N, x(i+Nb-1)=fliplr(ha)*e(i:i+Na-1)'-fliplr(hb(2:Nb)*x(i:i+Nb-2)'endx=x(5:end);% P=12;Q=6;arma1=armax(x',P Q);H,w=freqz(arma1.c,arma1.a,400); Hf=(abs(H).*sigma).2; Htol(:,k)=10*log10(Hf); figure(1); plot(w,10*log10(Hf),'g');hold on; endsum=0;for j=1:20; sum=sum+Htol(:,j);endfigure(1);plot(w,s