esprit4.m

此算法是基于空间旋转不变技术（ESPRIT）编写的matlab程序

t=0:0.01:100;
n=length(t);
ys=100+50*cos(2*pi*2.5*t).*exp(-0.2*t)+20*cos(2*pi*t).*exp(0.1*t)+2*randn(1,n);
N=200;
M=200;
dt=0.1;
x=zeros(1,N);
Rxx=zeros(M,M);
for i=1:N
    for m=1:M
        x(m)=ys(10*(i+m-1));
    end
    z=x'*x;
    Rxx=Rxx+(1/N)*z;
end
[V1,D1]=eig(Rxx);
w1=diag(D1);
w2=zeros(M,1);
for i=1:M
    w2(i)=w1(M-i+1);
end
t0=norm(w2,2);
for k=1:M
    fk=zeros(1,k);
    for c=1:k
        fk(c)=w2(c);
    end
       tk=norm(fk,2);
       vk=tk/t0;
    if vk>0.999999999
       pb=k;
       break;
    end
end
S=zeros(M,pb);
S1=zeros(M-1,pb);
S2=zeros(M-1,pb);
for i=1:M
    for j=1:pb
    S(i,j)=V1(i,M-j+1);
    end
end
for i=1:(M-1)
    for j=1:pb
        S1(i,j)=S(i,j);
    end
end
for i=1:(M-1)
    for j=1:pb
        S2(i,j)=S(i+1,j);
    end
end
tls=zeros(pb,pb);   
tls=inv(S1'*S1)*(S1'*S2);
[V2,D2]=eig(tls);
eigv=diag(D2);    %求取特征值
f=zeros(pb,1);
T=zeros(pb,1);
for i=1:pb
    f(i)=atan(imag(eigv(i))/real(eigv(i)))/(2*pi*dt);
    T(i)=-log(abs(eigv(i)))/dt;
end
b=eigv';
c=zeros(N,pb);
for j=1:pb
    c(1,j)=1;
end
for i=2:N
    for j=1:pb
        c(i,j)=b(j).^(i-1);
    end
end
A=zeros(pb,1);
A=inv(c'*c)*c'*x';
a=zeros(1,pb);
o=zeros(1,pb);
for i=1:pb
    a(i)=abs(A(i));                      %求取幅值和初相
    o(i)=atan(imag(A(i))/real(A(i)));
end