通信与信息系统专业张成根201112171933第二次作业.doc

报告二 AR模型谱估计一、 原理分析1、平稳ARMA过程：如果离散随机过程服从线性差分方程 （1）其中是一离散白噪声，则称为ARMA过程，式（1）所示的差分方程称为ARMA模型。该模型描述的是一个时不变的线性系统。2、平稳ARMA过程的功率谱密度：如果（1）式中的，则平稳ARMA过程的功率谱密度为 (2)3、Wold分解定理：任何一个具有有限方差的ARMA或MA过程都可以表示成唯一的、阶数有可能无穷大的AR过程；同样，任何一个ARMA或AR过程可以表示成一个阶数可能无穷大的MA过程。4、ARMA谱估计：该谱估计是使用个已知的观测数据计算出ARMA过程的功率谱密度估计值。若使用式（2）进行谱估计，需要事先辨识出整个ARMA模型及激励噪声的方差，而ARMA模型的辨识涉及到AR阶数、MA阶数的确定，以及AR参数和MA参数的估计。并且MA参数的估计需要求解非线性方程，不易实现，为此，我们可由Wold分解定理，将其简化为只有线性运算、需要辨识较少参量的谱估计过程。5、基于LS法的AR模型谱估计（1）给定AR的真实阶数，计算自相关函数。（2）由来生成维的Hankel矩阵 （3）（3）求解修正的YuleWalke方程： （4）得到 （5）其中 （4）对AR参数进行估计即可。6、基于SVD-TLS法的AR模型谱估计（1）对增广矩阵进行SVD分解，即，并存储奇异值和矩阵。其中 （6）（2）确定增广矩阵的有效秩。（3）由和来计算矩阵。其中是矩阵的第行、第列的元素，（4）解方程 （7）得到未知参数，进而实现总体最小二乘估计。二、 过程分析简化ARMA过程为AR过程AR阶数p是否已知计算计算由生成矩阵，并对进行SVD分解确定的有效秩由生成矩阵和解方程解方程求得AR参数实现AR谱估计计算矩阵结束开始YN三、结果分析1、LS法的AR模型谱估计图1 LS法和arccest函数法（AR阶数=4）图2 LS法和arccest函数法（AR阶数=100）2、SVD-TLS法的AR模型谱估计图3 SVD-LS法3、结果分析（1）对LS法和arccest函数调用法来说，当AR模型的阶数较小时，并且和比较接近时，两处的频谱混合在了一起，无法分辨出来，谱估计效果差。（2）对LS法和arccest函数调用法来说，适当增加AR模型的阶数，能够较好的将和两处的频谱分辨出来，谱估计效果好。（3）对SVD-TLS法来说，通过编程选择恰当的AR模型的阶数，也能够较好的将和两处的频谱分辨出来，谱估计效果好。四、实验心得通过这次实验，我对ARMA模型、AR模型以及AR模型的谱估计有了更深入的理解，掌握了实现AR模型谱估计的两种方法：LS法和SVD-TLS法，对这两种方法实现AR模型谱估计的区别有了更加清楚的认识，同时更加熟悉了MATLAB中的一些基本编程操作。附录：matlab程序代码1、LS法的程序代码：clcclear allclfp=4;N=128;n=1:N;w(n)=randn(1,N); %产生零均值方差为1的高斯白噪声x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n);rx=xcorr(x); %求信号x(n)的自相关序列for j=1:p R(j,:)=rx(j-p+N:j+N-1); %生成p阶Hankel矩阵endd=rx(1,N:2*N-1); %从信号x(n)的自相关序列提取后N个值for i=1:p r(i)=d(i+1);endls_a=-(R*R)(-1)*R*r; %LS法计算AR模型参数ls_a1=1 fliplr(ls_a); %fliplr讲行向量左右翻转h1,w1=freqz(1,ls_a1,1024,1); %AR模型的参数的频率特性估计subplot(211)plot(w1,abs(h1)xlabel(频率(Hz);title(LS法估计);legend(AR阶数=4);grid on%函数调用法avec=arrcest(x,p); %使用函数实现AR模型的参数估计h2,w2=freqz(1,avec,1024,1);subplot(212)plot(w2,abs(h2)xlabel(频率(Hz);title(arrcest函数法估计);legend(AR阶数=4);grid on2、SVD-TLS法的程序代码：clcclear allclfN=128; n=1:N;M=100; pe=30; w(n)=randn(1,N); %产生零均值方差为1的高斯白噪声xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n);rx=xcorr(xn,unbiased); %求信号x(n)的自相关序列 unbiased表示求的是无偏估计的相关函数for l=N-pe+1:N+M-1 Rx(l)=rx(l);endfor i=1:M for j=1:pe R(i,j)=Rx(N+i-j); endendr=rx(1,N+1:N+M);Re=r,R; %构造修正方程的增广矩阵ReU,S,V=svd (Re); %对矩阵Re进行svd分解for k=1:M a=S(k,k)/S(1,1); if(a0.05) break endend p=k; %确定矩阵Re的有效秩p sp=0; for j=1:p for i=1:(pe+1-p) sp=(S(j,j)2)*V(i:i+p,j)*V(i:i+p,j)+sp; %生成矩阵sp endendinv_sp=inv