esprit.m

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

%-------ESPRIT算法-----------%
clear;

%程序初始化
N=128;
M=20;
number=20;
delta=0;                              %噪声方差变量
F=zeros(number,M/2);                          
I=eye(M,M);
Z=zeros(M,M);

%产生矩阵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;                 %互相关矩阵

%对Rxx特征值分解，并计算矩阵Cxx以及Cxy                                      
Cxx=zeros(M,M);                          
Cxy=zeros(M,M);
[Ux,Sx,Vx]=svd(Rxx);
[Uy,Sy,Vy]=svd(Rxy);
delta=min(diag(Sx));                    %最小特征值为噪声方差
Cxx=Rxx-delta*I;                        %计算Cxx
Cxy=Rxy-delta*Z;                        %计算Cxy

%计算谐波频率                                                                   
F_num=0;
[V,S]=eig(Cxx,Cxy);
for i=1:M
    w=angle(S(i,i)); 
    if w>0
        F_num=F_num+1;
        F(butter,F_num)=w;
        F(butter,F_num)=(F(butter,F_num))/(2*pi);
    end    
end
F(butter,:)=sort(F(butter,:));

end


%求谐波频率统计结果
Freq=F
mean_F=mean(F)
var_F=std(F)