现代谱估计实验报告

杭州电子科技大学现代谱估计实验报告现代谱估计实验报告吴迪松 20040089 信息与信息处理一、实验题目：用MATLAB 编写一下算法的程序：(1)修正协方差法(2)多重信号分类（MUSIC ）算法(3)ESPRIT 算法(4)皮萨论科（Pisarenko ）谐波分解法并对各算法进行分析。

二、实验原理(1) 修正协方差法1、得出输入信号x(n),此输入信号是与白噪声相加的。

xx R ˆ()()()xn x n n t =+ 2、求自相关估计器(),xx c j k()()11**01,()(2()N N P xx n p n c j k x n j x n k x n j x n k N p −−−==)()⎡⎤=−−+++⎢⎥−⎣⎦∑∑ 3、求滤波器的系数 ˆ()ak1ˆ()(,)(,1)p xx xx k ak c l k c l ==−∑ 4、求白噪声方差的估计值2ˆσ21ˆˆ(1,1)()(1,)p xx xx k c a k c σ==+∑k 5、求功率谱()p f 2221ˆ()ˆ1()p j fk k p f ak e πσ−==+∑(2)多重信号分类（MUSIC ）算法1、估计样本自相关矩阵 ˆxx r2、对作特征值分解 ˆxx r3、确定的最小特征值分解的数目,求出这个最 ˆxx rE n E n 小特征值1,,p M λλ+"，令()2121p p M En σλλλ++=+++" 并求出相对应的特征向量,构造噪声矩阵1,,p v v +"M21,,p M V v v +⎡⎤=⎣⎦"4、构造函数：()211ˆMU M H i i p P f e v =+=∑(3)ESPRIT 算法1、利用已知信息求 {}(0),(1),()xx xx xx r r r M "2、由{}(0),(1),()xx xx xx r r r M "构造自相关矩阵xx R 和互相关矩阵xy R3、计算xx R 的特征值分解。

对于M>P ，最小特征值为噪声方差2σ4、利用2σ计算2xx xx C R I σ=−和2xy xy C R I σ=−5、计算矩阵对{},xx xy C C 的广义特征值分解。

从而确定谐波频率。

(4)皮萨论科（Pisarenko ）谐波分解法1、求x(n)的自相关函数xx r 2、求的toeplitz 矩阵xx r xx R3、求出xx R 的特征值，由此得出最小特征值4、求出最小特征值对应的一列特征向量5、求这个特征矢量形成的多项式的根。

由这些根求出它的角度，并且最后求出频率三、各算法的MATLAB 程序%AR model parameter estimateclearm=sqrt(-1);delta=0.1;a1=0.5;sample=32; %number of sample spotp=10; %number of sample spot in coef method;f1=0.05; f2=0.20; f3=0.45; %三个频率分量fstep=0.01;fstart=-0.5;fend=0.5;f=fstart:fstep:fend; %频率的范围nfft=(fend-fstart)/fstep+1;%un=urn+juinurn= normrnd(0,delta/2,1,sample);uin= normrnd(0,delta/2,1,sample);un=urn+m*uin;%calculate zn 噪声for n=1:sample-1zn(1)=un(1);zn(n+1)=a1*zn(n)+un(n+1);end%calculate xn 出入信号for n=1:samplexn(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1)) +sqrt(2)*real(zn(n));end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%calculate cxx(j,k) 公式3.102for j=1:pfor k=1:ps=0;for n=p+1:samples=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n-k)+xn(n-p+j)*conj(xn(n-p+k)));endcxx(j,k)=s;endend%calculate cxx0(j,1) 公式 3.103for j=1:ps=0;for n=p+1:samples=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n)+xn(n-p+j)*conj(xn(n-p))); endcxx0(j,1)=s;end%calculate a 公式3.104a=-inv(cxx)*cxx0;%a(k)*exp(-j2pifk)累加for i=1:length(f)sum=0;for k=1:psum=sum+a(k)*exp(-m*2*pi*f(i)*k);end%a(k)cxx(1,k)累加sun=0;for k=1:psun=sun+a(k)*cxx(1,k);enddelta1=cxx(1,1)+sun;%AR谱计的形式XIE(i)=delta1/(abs(1+sum))^2;endfiguresemilogy(f,XIE);title('修正协方差算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%xn= xn(:); x=xn';%calculate the value of rxxfor k=0:1:sample-1s=0;for n=1:sample-ks=s+conj(x(1,n))*x(1,n+k);endrxx(1,k+1)=(1/sample)*s;end%calculate the value of RxRx=zeros(sample,sample);Rx=toeplitz(rxx(1,1:32));[U,S,V]=svd(Rx); %特征值分解Pmusicf=zeros(1,1/fstep+1);ei=zeros(1,sample);for i=1:length(f)for j=1:sampleei(1,j)=exp(-2*pi*(j-1)*f(i)*m); %calculate the value of ei end;sum=0;for k=p+1:samplesum=sum+abs(ei*V(:,k))^2; %累加e*vi的绝对值endPmusicf(1,i)=(1/sum);endfiguresemilogy(f,Pmusicf);title('music算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%pp=6; %6 ge fen liangx=xn';M=length(x);rxx=xcorr(x,'biased'); % x de zi xiang guanrxx=[rxx(M:end),0]; %qu hou 32 weiR=toeplitz(rxx,rxx); %qiu Rxx,Rxy (公式5.158 5.159)Rxx=R(1:M,1:M);Rxy=R(1:M,2:end);D=eig(Rxx); %Rxx de te zheng zhi[delta2 i]=min(D); %qiu zhi xiao te zheng zhifor i=pp+1:Mey=eye(7,7);Z=ey(1:pp,2:end); %qiu ZCxx=Rxx(1:pp,1:pp)-delta2*ey(1:pp,1:pp); %p182 (4)Cxy=Rxy(1:pp,1:pp)-delta2*Z;D1=eig(Cxx,Cxy); %p182 (5)y=angle(D1)/2/pi;endfigurestem(y,ones(1,length(y)));title('esprit算法'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xt=xn';sin_num=3; %frequency fen liangN=length(xt);rxx=xcorr(xt,'biased');%x de zi xiang guangrxx=rxx(N:(2*sin_num)+N);%rxx de qu zhi%Frequencies estimationRxx=toeplitz(rxx); %qiu Rxxev=eig(Rxx); %qiu Rxx de te zheng zhi[S i]=min(ev); %qiu zhui xiao te zheng zhi[V D]=eig(Rxx); %te zheng zhi fen jiea=V(:,i); %qiu zui xiao te zheng zhi dui ying de yi lie te zheng xiang liangrts=roots(a); %qiu genw_est=[];for i=1:2*sin_numw_est(i)= angle(rts(i)); %qiu jiao duendF=(w_est/(2*pi))'; %qiu pin lvfigurestem(F,ones(1,length(F)),'*');title('Pisarenko算法');四、各算法MATLAB仿真后的分析四种谱估计算法中修正协方差法是对协方差法的一种改进，它可以得到高分辨率的统计稳定的谱估计值，但不能保证建立一个稳定的全极点滤波器。

而皮萨论科谐波分解法、MUSIC法都是噪声子空间频率估计，都是根据托布列兹自相关矩阵的噪声子空间特征矢量与信号矢量正交这一性质来求解的。

ESPRIT法是借助旋转不变技术估计参数，是一种主要特征分解法。

从上面四个图中可以明显的看出，它们对频率的估计都较精确。