esprit_tls.m

这是现代信号处理领域最常用的信号分析方法！

%--------------TLS-ESPRIT---------------
clear;

%程序初始化
N=128;
M=30;
number=20;
p=0;
delta=0;
rank=zeros(1,number);
I=eye(M,M);
Z=zeros(M,M);
UU=zeros(M,M,number);
SS=zeros(M,M,number);
VV=zeros(M,M,number);
CCxy=zeros(M,M,number);

%产生矩阵Z
for j = 1:M-1
    Z(j+1,j)=1;
end

%开始第一次循环
for butter=1:number

%产生观测信号数据
w=randn(1,1000) ;   
n=1:128;         
x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n);

%计算自相关矩阵Rxx和互相关矩阵Rxy                           
Rxx=zeros(M,M);                       
Rxy=zeros(M,M);
for n=1:N-M+1
    X=zeros(M,1);
    X=x(n:n+M-1);
    Rxx=Rxx+X'*X;               
end
for n=1:N-M
    X=zeros(M,1);
    Y=zeros(M,1);
    X=x(n:n+M-1);
    Y=x(n+1:n+M);
    Rxy=Rxy+X'*Y;                
end
Rxx=(1/(N-M+1))*Rxx;                   %自相关矩阵
Rxy=(1/(N-M))*Rxy;                     %互相关矩阵


%计算 Cxx 和 Cxy
Cxx=zeros(M,M);
Cxy=zeros(M,M);
[Ux,Dx,Vx]=svd(Rxx);
delta=min(diag(Dx));
Cxx=Rxx-delta*I;
Cxy=Rxy-delta*Z;
CCxy(:,:,butter)=Cxy;
[Uxx,Sxx,Vxx]=svd(Cxx);
UU(:,:,butter)=Uxx;
SS(:,:,butter)=Sxx;
VV(:,:,butter)=Vxx;
sum=0;
for i=1:M
    sum=sum+Sxx(i,i)^2;
end
overhead=0;
for i=1:M
overhead=overhead+Sxx(i,i)^2;
if (overhead/sum>0.9992)
    rank(1,butter)=i;
    break;
  end
end

  end                             %第一次循环结束
 p=round(mean(rank));
 

%估计频率
Lambda=zeros(1,p);
freq_num=0;
%第二次循环开始
for butter=1:number
    if rank(1,butter)==p
        freq_num=freq_num+1;
U=zeros(M,p);
S=zeros(p,p);
V=zeros(p,M);
C=zeros(p,p);
Cxy=zeros(M,M);
Cxy=CCxy(:,:,butter);
U=UU(:,1:p,butter);
S=SS(1:p,1:p,butter);
V=VV(:,1:p,butter);
C=U'*Cxy*V;
Lambda=eig(S,C);                         %广义特征值
w_num=0;
for i= 1:p
    w=angle(Lambda(i));
    if w>0
       w_num=w_num+1;
        F(freq_num,w_num)=w/(2*pi);
    end
end
F(freq_num,:)=sort(F(freq_num,:));
end
end                                  %第二次循环结束

%统计结果
Freq=F
Mean_F=mean(F)
Var_F=std(F)