esprit.m

来自「波束形成DOA估计中的ESPRIT算法matlab程序」· M 代码 · 共 73 行

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

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

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

%循环
for loops=1:loop

%产生观测信号数据
w=randn(1,128) ;   
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(loops,F_num)=w;
        F(loops,F_num)=(F(loops,F_num))/(2*pi);
    end    
end
F(loops,:)=sort(F(loops,:));

end


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