现代谱估计

现代谱估计实验报告1实验目的功率谱估计在实际工程中有重要应用价值。如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域发挥了重要作用。本次实验的目的主要是深入理解现代谱估计的基本理论，包括ARMA模型、ARMA谱估计。掌握现代谱估计的基本方法，包括SVD-TLS算法等。利用ARMA功率谱估计中Cadzow谱估计子和Kaveh谱估计子来进行谱估计。2实验原理2.1背景若离散随机过程{x(n)}服从线性差分方程)()()()(11jneneinxnxqjjpiiba（1）式中e（n）是一离散白噪声，则称{x(n)}为ARMA过程，而式（1）所示的差分方程称为ARMA模型。系数a1,a2……ap,和b1,b2……bq,分别称为自回归参数和滑动平均参数，而p和q分别叫做AR阶数和MA阶数。式（1）所示的ARMA过程，其功率谱密度为)()()()()(22eePjwjwzxBBezAzBwjw（2）ARMA谱估计的目的是使用N个已知的观测数据x(0)，x(1)…..x(N-1)计算出ARMA过程{x(n)}的功率谱密度估计。在实际中，可以运用cadzow谱估计子和kaveh谱估计子来估计，cadzow谱估计子秩序确定AR阶数p和估计AR参数，而kaveh谱估计子也只需要确定AR阶数p和估计AR参数以及MA阶数。2.2相关算法AR阶数p的确定用奇异值分解（SVD）,AR参数的估计用总体最小二乘法(TLS),即应用（SVD—TLS）算法来完成ARMA谱估计。SVD—TLS算法：步骤1计算增广矩阵B的SVD，并储存奇异值和矩阵V;步骤2确定增广矩阵B的有效秩p；步骤3计算矩阵S;步骤4求S的逆矩阵S--,并计算出未知参数的总体最小二乘估计。3实验内容仿真的观测数据由下式给出：xn=square(W*n)+0.2*randn(1,N)（3）其中，fs=20000，n=0:1/fs:0.1，N=length(n)，W=2000*pi。1、采样周期图法进行谱估计2、假设AR阶数未知，用SVD-TLS方法确定AR阶数和参数，然后使用Cadzow谱估计子进行谱估计。4Matlab仿真仿真的观测数据时域信号如图1所示。图1观测数据时域信号1、经典功率谱估计周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列。直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。仿真图如图2所示。图2周期图法功率谱2、现代功率谱估计现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。按照上面介绍的步骤，编写程序对观测信号x(n)进行仿真，可以设置不同的M，qe，pe的值，以便分析对比。图3是设置了M=2001，qe=100，pe=50，后得出的x(n)的功率谱图形。图3ARMA模型功率谱5实验总结本次实验分别用了周期图法和ARMA模型的参数估计方法对方波信号进行了功率谱估计，通过实验和得到的仿真图对比可以发现：通过周期图法得到的功率谱估计频谱分辨率较低，不能适应高分辨率功率谱估计的要求，参数化的谱估计可以获得高频率分辨率的功率谱。经典功率谱估计的分辨率反比于有效信号的长度，但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N点有限长序列x(n)，虽然其估计出的自相关函数也是有限长的，但是现代谱估计的一些隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度，不像经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高，而且现代谱线要平滑得多，从上图可以清楚看出。6附录Matlab程序如下：main.mclear;closeallfs=20000;n=0:1/fs:0.1;N=length(n);W=2000*pi;x1n=square(W*n);x2n=randn(1,N);xn=x1n+0.2*x2n;figure;plot(n,xn);title('时域信号');Nfft=100;[Pxx,f]=period(xn,fs,Nfft);figure;plot(f,Pxx);title('周期图法功率谱');%ARMA谱估计pe=50;qe=100;NARMA=length(xn);M=length(n);[a,Rx,p]=ARMA(xn,qe,pe,M);%Cadzow谱估计子[Pw]=Cadzow(a,Rx,p,NARMA);%功率谱×figure;plot((0:length(Pw)-1)*fs/length(Pw),Pw);title('ARMA模型');period.mfunction[Pxx,f]=period(xn,fs,Nfft)Pxx=abs(fft(xn,Nfft).^2)/Nfft;f=(0:length(Pxx)-1)*fs/length(Pxx);ARMA.mfunction[a,Rxx,p]=ARMA(xn,qe,pe,M)Rxx=xcorr(xn,'unbiased');for(i=1:M)for(j=1:pe+1)Re(i,j)=Rxx(pe+i+1-j);endend[U,S,V]=svd(Re);Ak=0;for(i=1:pe+1)Ak=Ak+S(i,i)^2;end;Akf=0;v=Akf/Ak;p=0;whilev0.99995p=p+1;Akf=Akf+S(p,p)^2;v=Akf/Ak;endSp=0;forj=1:pfori=1:pe+1-pSp=Sp+S(i,i)^2*V(i:i+p,j)*V(i:i+p,j)';endendinvSp=inv(Sp);for(i=1:p)a(i)=invSp(i+1,1)/invSp(1,1);endCadzow.mfunction[Pw]=Cadzow(a,Rx,p,N)Ro=Rx;Ro(N)=0.5*Rx(N);p=p-1;n=zeros(1,p+1);fork=1:p+1n(k)=a(1:p+1)*Ro(k+N:-1:k-p+N)';endwfreq=linspace(0,2*pi,N);z=exp(sqrt(-1)*wfreq);Pw=polyval(n,z)./polyval(a,z)+polyval(n,z.^-1)./polyval(a,z.^-1);Pw=abs(Pw);