AR模型和ARMA模型谱估计仿真

1、AR 模型和 ARMA 模型谱估计仿真一、问题重述有两个 ARMA 过程，其中信号1 是宽带信号，信号 2 是窄带信号，分别用AR 谱估计算法、 ARMA 谱估计算法和周期图法估计其功率谱。产生信号 1 的系统函数为：1 + 0.3544?-1 + 0.3508?-2+ 0.1736?-3+ 0.2401?-4H( z) =1.3817?-1 + 1.5632?-2- 0.8843?-3+ 0.4906?-41 -激励白噪声的方差为1.产生信号 2 的系统函数为：H( z) =1 + 1.5857?-1+ 0.9604?-21.6408?-1 + 2.2044?-2- 1.4808?-3+ 0

2、.8145?-41 -激励白噪声的方差为1.每次实验使用的数据长度为256.二、模型分析很多随机过程可以由或近似由均值为零、方差为2?的白噪声序列 u( n)输函数 H(z) 的 ARMA 线性系统来得到。称该随机过程为ARMA 过程。B(z)+U(z)-1?(?)经过具有有理想传X(z)?-?(?)H( z) =?=0?-?=?(?)?=0 ?()2| 2? =? ? ?(?)由上式可知只要估计出模型的参数（?和 ?），即可求出功率谱。?1.AR 模型的建立：AR 模型是一种特殊的ARMA 模型，利用AR（ p）模型，即：?x( n) = - ?( ?- ?) + ?(?)?=1逼近采样样本

3、 ,此时功率谱表达式为：2?2?=?-?|1+? |?=1?需要求解得未知量为参数?已知时，利用 x(n) 的自相关函数与AR 模型参数，当阶数 p?的关系，可建立 Y-W 方程，解该方程，即可得到AR 参数。?()()? +?- ?= 0?=1?( ?) + ?2?( ?- ?) = ?=1?= 1,2?= 0?(0)?(-1)?(-?)12?(1)?(0)?(1-?)?=01? ?(?) ?(?-1)?(0)?0对于 (p+1)元线性方程，若采用matlab 中的函数，则使用的是高斯消去法运算量为3?数量级。采用下面的 Levinson-Durbin 算法。它可以将运算量减少到2?数量级。

4、递推公式的获取方法如下：令 p=1，得一阶 AR模型对应的尤勒沃克方程为()+ ?1()()= ?1? 01? 1()+ ?1()(0)= 0? 11?可解得?( 1)?(1) = -1?( 0)?()2 (1)0-?1 -2(1)?1 = ?( 0)= ?0?1?令 p=2,得二阶 AR 模型对应的尤勒沃克方程为?(0) + ?(1)?(1) + ?(2)?(2) = ?2?2?2()+ ?2() ()+ ?2(2)(1)= 0 ? 11 ?0? ( 2) + ?( 1) ?( 1) + ?( 2)?(0) = 0? 2?2?解此方程组，得?2(2) =-? ?( 2) + ?1 (1)?(

5、1)/? 1?(1) =?( 1) +?( 2) ?( 1)21221?21 -(2)= ?1?2令 p=3,4, ,由递推规律可得到下式：?-1? = -?()(?-?+ ?)/?-1?-1(?)?=1( )( )()(?-)1,2, ,?- 1? ?= ?-1?+ ?-1?,?=?= ?21 - ?-1?式中， ?1,2, , ? ?=0 ?(0) .?= ?( ?) ,? =按照此递推式即可方便解出Ar 模型参数 a。2.ARMA 模型ARMA(p,q) 的差分方程如下：?x( n ) = - ?(?-?)+ ?(?- ?)?=1?=0由此可得出自相关的关系：?2? ( ?- ?) =

6、? (?)? ?+?=0?=0?-?-()2?()? -?+ ? ?+?=1?=0?-?( )()2?( )? =-? -?+ ?+?=1?=-?= 0,1,2,?< 0- ?(?-?)>? ?=1当 m>q 时的自相关函数关系可得 ,建立超定方程（ M>p+q）：R(q)?(?- 1)?(?-?-1)?1-?(?+ 1) ?(?+1)?(?)?(?-?+2)?=?-?(?)?(M -1)?(M - p)?(?)R ?a = r利用最小二乘解得：?=?-1?(? ?)? ?利用 Kaveh 谱估计方法，在计算出a 后直接计算参数 ?，则此时功率谱估计式变为：?-?(?

7、?)?=-? ?=-1?(?)?(?而参数 ?可以由以下关系式求出：? ?()?= 0,1, ?= ?- ?+ ?=0?=0三、实验内容信号一：按照题目要求要对信号一分别使用AR(4),AR(8),ARMA(4,4),ARMA(8,8)模型和周期图法进行估计，做 20 次独立实验结果绘制在一张图上， 并绘制 20 次的平均谱和真实谱。 信号一是均值为 0，方差为 1 的白噪声通过如下系统函数产生的。H (z)得到如下实验结果：10.3544 z1 1.3817 z110.3508z1.5632 z220.1736 z0.8843z330.2401 z0.4906 z4420离散系统功率谱 AR

8、(4)15单次估计频谱10平均频谱实际频谱50-5-10-15-200.511.522.533.5025离散系统功率AR(8)谱20单次试验估计普15平均估计谱10实际功率谱50-5-10-15-200.511.522.533.50离散 系统功率谱 ARMA(8,8)200-20-40单次试验估计普平均估计谱-60实际功率谱-80-100-1200.511.522.533.50离散系统功率谱ARMA(4,4)2015105单次试验估计普平均估计谱0实际功率谱-5-10-15-2000.511.522.533.5离散系统功率谱谱估计法3020100-10-20-30-40-50-6000.511

9、.522.533.5信号二：按 照 题 目 要 求 要 对 信 号 二 分 别 使 用AR(4),AR(8),AR(12),AR(16),ARMA(4,2),ARMA(8,4),ARMA(12,6)模型和周期图法进行估计，做20 次独立实验结果绘制在一张图上，并绘制20次的平均谱和真实谱。信号一是均值为0，方差为 1 的白噪声通过如下系统函数产生的。H (z)11.5857 z 112.2044 z 21 1.6408 z实验结果如下图所示：0.9604z1.4808z230.8145 z 4离散系统功率谱AR(4)403020100-10-20-30-40-5000.511.522.533.

10、5离散系统功率谱AR(8)403020100-10-20-30-40-5000.511.522.533.5离散系统功率谱AR(12)403020100-10-20-30-40-5000.511.522.533.5离散系统功率谱 AR(16)403020100-10-20-30-40-500.511.522.533.50离散 系统 功率 谱ARMA(4,2)403020100-10-20-30-40-50-6000.511.522.533.5离 散 系 统 功 率 谱 ARMA(8,4)40200-20-40-60-80-10000.511.522.533.5离散系统功率谱ARMA(12,6)4

11、03020100-10-20-30-40-5000.511.522.533.5离散系统功率谱谱估计法6040200-20-40-6000.511.522.533.5四、结果分析为比较算法的不同，现对两种不同的信道做以分析，信道零极点图如下：信道一零极点图信道二零极点图traPryanigamI21.510.50-0.5-1-1.5-2traPryanigamI21.510.50-0.5-1-1.5-22-1-0.500.51Real Part-1-0.500.51Real Part由此可知信道二比信道一更不稳定。实验结果展示出采用不同方法有着不同的特点，周期图法做出的谱估计图最差，毛刺较多，对

12、稳定性较差的系统估计不是很好。 AR 模型有较好的性能，估计谱比较平滑，但对稳定性较差的系统的估计也不是很好。 ARMA 模型估计谱也比较平滑， 其性能的最主要优点体现在其对稳定性较差的系统的估计能力比较强， 但是在每次估计时会出现某些极差的数据点，其算法稳定性较差。同时通过实验可以看出， AR 和 ARMA 的性能并不是阶数越大越好，所以应该采用一定的方法确定其阶数，能得到较好效果。附件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=

13、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=zeros(400,20);% 产生 500 点数据for k=1:20,%做 20 次谱估计实验N=500;ha=cell2mat(h(1);hb=cell2mat(h(2);Na=length(ha);%噪声项系数长度Nb=length(hb);% 信号项系数长度x=zeros(1,N+Nb-1);e=wgn(1,N+Na-1,sigma,'line

14、ar','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);% 奇异值分解法确定模型阶数% 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=ei

15、g(Cxx);% figure;% plot(1:N,V);V(V<0.5)=0;P=0;for i=1:N,if V(i)>0P=P+1;endend% 计算 AR 模型参数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

16、)2);endHa,wa=freqz(1,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,

17、10*log10(Hf);Je=num2str(P);str1=strcat(' 离散系统功率谱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),cel

18、l2mat(h(2),400);Htol=zeros(400,20);% 产生 500 点数据for k=1:20,%做 20 次谱估计实验N=500;ha=cell2mat(h(1);hb=cell2mat(h(2);Na=length(ha);%噪声项系数长度Nb=length(hb);% 信号项系数长度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');h