esprit-ls-2.m

现代信号处理中对基本的ESPRIT的一种改进算法的程序

%-----------------基本ESPRIT算法2(3.10.3)-----------%
clear
loops = 20;
N = 100;        % N是快拍次数
m = 100;         % m是阵元个数
for loop = 1:loops
    n = [1:200];
    w = randn(1,200);
    x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+w;    % x是[1,200]的行向量，生成观测数据
    for n = 1:N
        X(:,n) = x(n:(n+m-1));                        % 生成 X[X(1),X(2),...X(N)],X是(m*N)的矩阵
    end
    % 产生协方差矩阵
    Cov_Rxx = cov(X);
    % 计算有效秩，确定主特征值个数p
    [U,S,V] = svd(Cov_Rxx);
    Af = 0;
    for i = 1:m
        Af = Af + S(i,i)^2;
    end
    Akf = 0;
    for i = 1:m
        Akf = Akf + S(i,i)^2;
        v = sqrt(Akf/Af);
        if v>0.999
            p = i;
            break;
        end
    end
    Us = U(:,1:p);
    J1 = [eye(m-1),zeros((m-1),1)];              % 选择矩阵，选出特定的行
    J2 = [zeros((m-1),1),eye(m-1)];
    Us1 = J1*Us;
    Us2 = J2*Us;
    % 计算谐波频率
    delta = inv(Us1'*Us1)*Us1'*Us2;
    [A,B] = eig(delta);
    k = 0;
    for i = 1:p
        W = angle(B(i,i));                   % 因为特征根是以共轭对的形式存在，故只需选出正的频率即可
        if W > 0
           k = k+1;
           f(k) = angle(B(i,i))/(2*pi);
        end
    end
    f;
    F(:,loop) = f';
end
% 统计结果
p
F
F_mean = mean(F,2)
F_std =  std(F')'