esprit-tls.m

用于谐波恢复的ESPRIT方法的SVD-TLS的程序实现

%-------------TLS—ESPRIT(3.10.2)---------------%
clear all
loops = 20;     % loops 是独立实验次数
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
    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估计
    % [U,S,V] = svd(Rxx);
    % var = S(1,1);
    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的有效秩
    [U,S,V] = svd(Cxx);
    Af = 0;
    for i = 1:m
        Af = Af + S(i,i)^2;
    end
    Akf = 0;
    % p = 0;
    % v = 0;
    % while v<0.999
      % p = p+1;
      % Akf = Akf + S(p,p)^2;
      % v = sqrt(Akf/Af);
    % end
    for i = 1:m
        Akf = Akf + S(i,i)^2;
        v = sqrt(Akf/Af);
        if v>0.999
            p = i;
            break;
        end
    end
    % 存储U1,V1,S1
    U1 = U(:,[1:p]);
    S1 = S([1:p],[1:p]);
    V1 = V(:,[1:p]);
    R1 = S1;
    R2 = U1'*Cxy*V1;
    % 对矩阵束进行广义特征值分解
    [A,B] = eig(R1,R2);
    % w = eig(R1,R2);
     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 = sort(f);
    F(:,loop) = f';
end
% 统计结果
p
F
F_mean = mean(F,2)                      % 参数2表示，对矩阵F_mean的各行求均值
F_std =  std(F')'