esprit_hz.m

本代码的方法是借助旋转不变技术估计信号参数(Estimation of Signal Parameters via Rotational Invariance Techniques)

clear all
m=sqrt(-1);
delta=0.101043;
a1=-0.850848;
sample=32;      %采样点
p=10;           %number of sample spot in coef method;
%f1=0.05; f2=0.40; f3=0.42;
f1=0.05; f2=0.30; f3=0.42;
fstep=0.01;
fstart=-0.5;
fend=0.5;  
f=fstart:fstep:fend;
nfft=(fend-fstart)/fstep+1;

%构造噪声：高斯白噪声，正态分布，复数un=urn+juin
urn= normrnd(0,delta/2,1,sample);
uin= normrnd(0,delta/2,1,sample);
un=urn+m*uin;

%计算 zn
for n=1:sample-1
    zn(1)=un(1);
    zn(n+1)=-a1*zn(n)+un(n+1);
end
    
%求xn
    for n=1:sample
        xn(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

    pp=6;       %  假定值
x=xn';
M=length(x);
rxx=xcorr(x,'biased');    %求相关函数仅为向量
R=toeplitz(rxx(1:M+1));   %变相关距阵
Rxx=R(1:M,1:M);           %取一到M行,1到M列
Rxy=R(1:M,2:M+1);         %取一到M行,2到M+1列
[U,S,V]=svd(Rxx);         %求特征值
U1=U(:,1:pp);             %取所有行1到pp列
S1=S(1:pp,1:pp);          %取1到pp行1到pp列
V1=V';
V1=V1(1:pp,1:32);
A=U1'*Rxy*V1';
d=eig(S1,A);              %求特征向量
y=angle(d)/2/pi;
w=[1:6];
figure(3)
stem(w,y);
title('esprit算法f2＝0.30时');