esprit.m

谱估计中的ESPRIT算法仿真

NFFT=1024;
N=32;
M=10;
p=6;
n=0:1:31;
f1=0.05;
f2=0.40;
f3=0.42;
u=0.101043;
U1=wgn(1,N,10*log10(u/2));
U2=wgn(1,N,10*log10(u/2));
U=U1+U2*j;
a=[1 -0.850848];
b=[1];
z=filter(b,a,U);
x=2*cos(2*pi*f1*n)+2*cos(2*pi*f2*n)+2*cos(2*pi*f3*n)+z;
[X,R] = corrmtx(x,M);        %估计样本的自相关矩阵 M+1 × M+1
Rxx = R(1:M,1:M);             %计算Rxx
Rxy = R(1:M,2:M+1);           %计算Rxy
[U,S,V]=svd(Rxx);             %对Rxx进行奇异值分解
Sigma1=S(1:p,1:p);            %取出p个主奇异值对应的Sigma1，U1和V1
U1=U(:,1:p);                  
V1=V(:,1:p);
temp=U1'*Rxy*V1;              %计算p×p维矩阵U1H*Rxy*V1
[AA,BB,Q,Z]=qz(Sigma1,temp);  %进行广义特征值求解
alpha = diag(AA);             %取出广义特征值
f=angle(alpha)*pi/2;          %估计出谐波频率分量
f=f/10;
stem(f,(ones(size(f))));        %画图