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

AR模型和ARMA模型谱估计仿真一、问题重述有两个ARMA过程，其中信号1是宽带信号，信号2是窄带信号，分别用AR谱估计算法、ARMA谱估计算法和周期图法估计其功率谱。产生信号1的系统函数为：H(z)=1+0.3544𝑧−1+0.3508𝑧−2+0.1736𝑧−3+0.2401𝑧−41−1.3817𝑧−1+1.5632𝑧−2−0.8843𝑧−3+0.4906𝑧−4激励白噪声的方差为1.产生信号2的系统函数为：H(z)=1+1.5857𝑧−1+0.9604𝑧−21−1.6408𝑧−1+2.2044𝑧−2−1.4808𝑧−3+0.8145𝑧−4激励白噪声的方差为1.每次实验使用的数据长度为256.二、模型分析很多随机过程可以由或近似由均值为零、方差为𝛿2的白噪声序列u(n)经过具有有理想传输函数H(z)的ARMA线性系统来得到。称该随机过程为ARMA过程。H(z)=∑𝑏𝑖𝑧−𝑖𝑞𝑖=0∑𝑎𝑖𝑧−𝑖𝑝𝑖=0=𝐵(𝑧)𝐴(𝑧)𝑃𝑥𝑥(𝑤)=𝛿2∙|𝐻(𝑤)|2由上式可知只要估计出模型的参数（𝑎𝑖和𝑏𝑖），即可求出功率谱。1.AR模型的建立：AR模型是一种特殊的ARMA模型，利用AR（p）模型，即：x(n)=−∑𝑎𝑖𝑥(𝑛−𝑖)+𝑢(𝑛)𝑝𝑖=1逼近采样样本,此时功率谱表达式为：𝑃̂𝑥=𝛿2|1+∑𝑎𝑖𝑒−𝑗𝑤𝑖𝑃𝑖=1|2需要求解得未知量为参数𝑎𝑖，当阶数p已知时，利用x(n)的自相关函数与AR模型参数的关系，可建立Y-W方程，解该方程，即可得到AR参数。{𝑅𝑥(𝑚)+∑𝑎𝑖𝑅𝑥(𝑚−𝑖)=0𝑚=1,2……𝑝𝑝𝑖=1𝑅𝑥(𝑚)+∑𝑎𝑖𝑅𝑥(𝑚−𝑖)=𝛿2𝑚=0𝑝𝑖=1𝐴(𝑧)−1B(z)+－X(z)U(z)[𝑅(0)𝑅(−1)𝑅(1)𝑅(0)…𝑅(−𝑝)𝑅(1−𝑝)⋮⋱⋮𝑅(𝑝)𝑅(𝑝−1)⋯𝑅(0)]∙[1𝑎1⋮𝑎𝑝]=[𝛿20⋮0]对于(p+1)元线性方程，若采用matlab中的函数，则使用的是高斯消去法运算量为𝑝3数量级。采用下面的Levinson-Durbin算法。它可以将运算量减少到𝑝2数量级。递推公式的获取方法如下：①令p=1，得一阶AR模型对应的尤勒—沃克方程为{𝑅𝑥(0)+𝑎1(1)𝑅𝑥(1)=𝜌1𝑅𝑥(1)+𝑎1(1)𝑅𝑥(0)=0可解得𝑎1(1)=−𝑅𝑥(1)𝑅𝑥(0)𝜌1=𝑅𝑥(0)−𝑅𝑥2(1)𝑅𝑥(0)=𝜌0[1−𝑎12(1)]②令p=2,得二阶AR模型对应的尤勒—沃克方程为{𝑅𝑥(0)+𝑎2(1)𝑅𝑥(1)+𝑎2(2)𝑅𝑥(2)=𝜌2𝑅𝑥(1)+𝑎2(1)𝑅𝑥(0)+𝑎2(2)𝑅𝑥(1)=0𝑅𝑥(2)+𝑎2(1)𝑅𝑥(1)+𝑎2(2)𝑅𝑥(0)=0解此方程组，得{𝑎2(2)=−[𝑅𝑥(2)+𝑎1(1)𝑅𝑥(1)]/𝜌1𝑎2(1)=𝑎1(1)+𝑎2(2)𝑎1(1)𝜌2=𝜌1[1−𝑎22(2)]③令p=3,4,…,由递推规律可得到下式：{𝑘𝑝=−[𝑅𝑥(𝑝)+∑𝑎𝑝−1(𝑖)𝑅𝑥(𝑝−𝑖)/𝜌𝑝−1𝑝−1𝑖=1]𝑎𝑝(𝑖)=𝑎𝑝−1(𝑖)+𝑎𝑝(𝑝)𝑎𝑝−1(𝑝−𝑖),𝑖=1,2,…,𝑝−1𝜌𝑝=𝜌𝑝−1[1−𝑘𝑝2]式中，𝑘𝑝=𝑎𝑝(𝑝),𝑚=1,2,…,𝑝𝜌0=𝑅𝑥(0).按照此递推式即可方便解出Ar模型参数a。2.ARMA模型ARMA(p,q)的差分方程如下：x(n)=−∑𝑎𝑖𝑥(𝑛−𝑖)+∑𝑏𝑖𝑞𝑖=0𝑢(𝑛−𝑖)𝑝𝑖=1由此可得出自相关的关系：∑𝑎𝑖𝑅𝑥(𝑚−𝑖)=𝜎2∑ℎ∗(𝑘)𝑏𝑘+𝑚∞𝑘=0𝑝𝑖=0𝑅(𝑚)={−∑𝑎𝑖𝑅(𝑚−𝑖)+𝜎2∙∑ℎ∗(𝑖)𝑏𝑘+𝑚𝑞−𝑚𝑖=0𝑚=0,1,2,…𝑞𝑝𝑖=1−∑𝑎𝑖𝑅(𝑚−𝑖)+𝜎2∙∑ℎ∗(𝑖)𝑏𝑖+𝑚𝑞−𝑚𝑖=−𝑚𝑝𝑖=1𝑚0−∑𝑎𝑖𝑅(𝑚−𝑖)𝑝𝑖=1𝑚𝑞当mq时的自相关函数关系可得,建立超定方程（Mp+q）：[R(q)𝑅(𝑞−1)𝑅(𝑞+1)𝑅(𝑞)…𝑅(𝑞−𝑝−1)𝑅(𝑞−𝑝+2)⋮⋱⋮𝑅(M−1)𝑅(M−p)⋯𝑅(𝑞)]∙[𝑎1⋮𝑎𝑝]=[−𝑅(𝑞+1)⋮−𝑅(𝑀)]R∙a=r利用最小二乘解得：𝑎̂=(𝑅𝐻𝑅)−1𝑅𝐻𝑟利用Kaveh谱估计方法，在计算出a后直接计算参数𝑐𝑘，则此时功率谱估计式变为：𝑃𝑥(𝑧)=∑𝑐𝑘𝑧−𝑘𝑞𝑘=−𝑞𝐴(𝑧)𝐴(𝑧−1)而参数𝑐𝑘可以由以下关系式求出：𝑐𝑘=∑∑𝑎𝑖𝑎𝑖∗𝑅𝑥(𝑘−𝑖+𝑗)𝑘=0,1,…𝑞𝑝𝑗=0𝑝𝑖=0三、实验内容信号一：按照题目要求要对信号一分别使用AR(4),AR(8),ARMA(4,4),ARMA(8,8)模型和周期图法进行估计，做20次独立实验结果绘制在一张图上，并绘制20次的平均谱和真实谱。信号一是均值为0，方差为1的白噪声通过如下系统函数产生的。1234123410.35440.35080.17360.2401()11.38171.56320.88430.4906zzzzHzzzzz得到如下实验结果：00.511.522.533.5-20-15-10-505101520离散系统功率谱AR(4)单次估计频谱平均频谱实际频谱00.511.522.533.5-20-15-10-50510152025离散系统功率谱AR(8)单次试验估计普平均估计谱实际功率谱00.511.522.533.5-120-100-80-60-40-20020离散系统功率谱ARMA(8,8)单次试验估计普平均估计谱实际功率谱00.511.522.533.5-20-15-10-505101520离散系统功率谱ARMA(4,4)单次试验估计普平均估计谱实际功率谱00.511.522.533.5-60-50-40-30-20-100102030离散系统功率谱谱估计法信号二：按照题目要求要对信号二分别使用AR(4),AR(8),AR(12),AR(16),ARMA(4,2),ARMA(8,4),ARMA(12,6)模型和周期图法进行估计，做20次独立实验结果绘制在一张图上，并绘制20次的平均谱和真实谱。信号一是均值为0，方差为1的白噪声通过如下系统函数产生的。12123411.58570.9604(z)11.64082.20441.48080.8145zzHzzzz实验结果如下图所示：00.511.522.533.5-50-40-30-20-10010203040离散系统功率谱AR(4)00.511.522.533.5-50-40-30-20-10010203040离散系统功率谱AR(8)00.511.522.533.5-50-40-30-20-10010203040离散系统功率谱AR(12)00.511.522.533.5-50-40-30-20-10010203040离散系统功率谱AR(16)00.511.522.533.5-60-50-40-30-20-10010203040离散系统功率谱ARMA(4,2)00.511.522.533.5-100-80-60-40-2002040离散系统功率谱ARMA(8,4)四、结果分析为比较算法的不同，现对两种不同的信道做以分析，信道零极点图如下：00.511.522.533.5-50-40-30-20-10010203040离散系统功率谱ARMA(12,6)00.511.522.533.5-60-40-200204060离散系统功率谱谱估计法由此可知信道二比信道一更不稳定。实验结果展示出采用不同方法有着不同的特点，周期图法做出的谱估计图最差，毛刺较多，对稳定性较差的系统估计不是很好。AR模型有较好的性能，估计谱比较平滑，但对稳定性较差的系统的估计也不是很好。ARMA模型估计谱也比较平滑，其性能的最主要优点体现在其对稳定性较差的系统的估计能力比较强，但是在每次估计时会出现某些极差的数据点，其算法稳定性较差。同时通过实验可以看出，AR和ARMA的性能并不是阶数越大越好，所以应该采用一定的方法确定其阶数，能得到较好效果。-1-0.500.51-2-1.5-1-0.500.511.52RealPartImaginaryPart信道一零极点图-1-0.500.51-2-1.5-1-0.500.511.522RealPartImaginaryPart信道二零极点图附件AR模型法clear;clc;echooff;closeall;h1={[10.35440.35080.17360.2401][1-1.38171.5632-0.88430.4906]};h2={[11.58570.9604][1-1.64082.2044-1.48080.8145]};sigma=1;h=h2;%%绘制信道功率谱[H,w]=freqz(cell2mat(h(1)),cell2mat(h(2)),400);Htol=zeros(400,20);%%产生500点数据fork=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');fori=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);%%奇异值分解法确定模型阶数%fork=1:20,temp=xcorr(x,'biased');cor_l=round(length(temp)/2);Cxx=zeros(cor_l,cor_l);fori=1:cor_lCxx(i,:)=temp(cor_l+1-i:length(temp)+1-i);end[X,B]=eig(Cxx);V=diag(B)./max(diag(B));%归一化特征值%figure;%plot(1:N,V);V(V0.5)=0;P=0;fori=1:N,ifV(i)0P=P+1;endend%%计算AR模型参数P=16;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);fori=2:P,%forj=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);end[Ha,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');holdon;endsum=0;forj=1:20;sum=sum+Htol(:,j);endfigure(1);plot(wa,sum/20,