esprit-ls-1.m

《现代信号处理》中关于利用基本的ESPRIT算法估计信号参数（频率和波达方向角）的仿真程序

%  -----------基本ESPRIT算法1（3.10.1）———————— %
clear all
loops = 20;     % loops 是独立实验次数
N = 100;        % N是快拍次数
m = 100;         % m是阵元个数
% for loop = 1:loops
    n = [1:200];
    w = randn(1,200);
    % w = randn(1,2000);
    % x = sqrt(20)*sin(2*pi*0.2*n) +
    % sqrt(2)*sin(2*pi*0.213*n)+w(1,n+50*loop);   % x是[1,200]的行向量，生成观测数据
    x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+w;
    for n = 1:N
        X(:,n) = x(n:(n+m-1));                        % 生成 X[X(1),X(2),...X(N)],X是(m*N)的矩阵,每一列为一次观察
    end
    for n = 1:N
        Y(:,n) = x((n+1):(n+m));                       % 生成y(n)=x(n+1)，Y为（m*N）的矩阵
    end
    %生成Rxx与Rxy
    Rxx = 0;
    for i = 1:N
        Rxx = Rxx+X(:,i)*X(:,i)';
    end 
    Rxx = Rxx/N;                                       % Rxx和Rxy都是（m*m）的矩阵
    Rxy = 0;
    for i = 1:N
        Rxy = Rxy+X(:,i)*Y(:,i)';                  
    end                   
    Rxy = Rxy/N;
    [A,B] = eig(Rxx);                          % 对Rxx进行特征值分解
    var = min(diag(B));                        % 取最小特征值作为噪声方差var估计
    % 计算Cxx与Cxy
    I = eye(m);
    for i = 1:(m-1)
        c(i)=1;
    end
    Z = diag(c,-1);                        
    Cxx = Rxx - I*var;
    Cxy = Rxy - Z*var;                       % Cxx与Cxy均为(m*m)的矩阵
    % 矩阵束(Cxx,Cxy)的广义特征分解
    [V,S] = eig(Cxx,Cxy);
    k = 0;
    for i = 1:m
        W = angle(S(i,i));                   % 因为特征根是以共轭对的形式存在，故只需选出正的频率即可
        if W > 0
           k = k+1;
           S_abs(k) = abs(S(i,i));
           f(k) = angle(S(i,i))/(2*pi);      % 相位差W为什么就是角频率w? sin(w(k+1))=sin(wk+w);
        end
    end
    % 求谐波频率
    S_abs_num = numel(S_abs);
    for i = 1:S_abs_num
        delta_S(i) = abs(1-S_abs(i));
    end
    delta_S = sort(delta_S);
    for i = 1:S_abs_num
        if abs(1-S_abs(i)) == delta_S(1);
            f_value1 = f(i);
        end
    end
    for i = 1:S_abs_num
        if abs(1-S_abs(i)) == delta_S(2);
            f_value2 = f(i);
        end
    end
    f_value1
    f_value2
    % F(:,loop) = [f_value1,f_value2]';                        % 统计的频率矩阵—f 
%end
% 计算统计结果
% F
% f_mean = mean(F,2)
% f_var = std(F')'