tls_esprit.m

旋转不变技术估计信号参数(ESPRIT)算法是子空间估计的两个经典方法,它们都能够有效地估计特征子空间。程序包里面实现了用LS估计和TLS估计的ESPRIT算法。

function [bias  var fnum] =TLS_ESPRIT(x, n, th)
%TLS_ESPRIT算法估计频率值

bias = zeros(1,4);
bias1 = zeros(100,4);
var= zeros(1,4);
var1= zeros(100,4);
ws=[-0.213 -0.2 0.2 0.213];
fnum =zeros(1,4);
th2=0.01;
th3=0.05;
p=0;
count = zeros(1,4);

for kj=1:100
    x1=x(kj,:);
    N=numel(x1);
    %计算自相关函数
    r = zeros(1,N);
    for k = 1:N
        for l = k : N
            r(k)=r(k) + x1(l)*x1(l-k+1);
        end
    end
    r = r / N;
    %构造自相关函数矩阵与互相关函数矩阵
    r1 = r(1,n+1:-1:2);
    r2 = [r1,r(1,1:n+1)];
    Rxx=zeros(n,n);
    Rxy = zeros(n,n);
    for i = 1:n
        Rxx(i,:) = r2(1,n-i+2:2*n-i+1);%确定自相关矩阵
        Rxy(i,:) = r2(1,n-i+3:2*n-i+2);%确定互相关矩阵
    end
    
    %对Rxx进行奇异值分解
    d=eig(Rxx);
    d=sort(d,'descend');
    q=0;
    for i=1:n
       if d(i)/d(1)<th2
            q=i-1;
            break;
       end
    end
    if q~=0
        noise=sum(d(q+1:n))/(n-q); %得到噪声功率
    else 
        noise=min(d(i));
    end
    
    %计算Cxx和Cxy
    Z = zeros(n,n);
    Z(2:n,1:n-1) = eye(n-1,n-1);
    Cxx = Rxx - noise * eye(n,n);
    Cxy = Rxy - noise * Z;
    
    %奇异值分解Cxx确定阶数
    [U,S,V] = svd(Cxx);
    for k = 2 : n
        v = S(k,k)/S(1,1);
        if v < th
            p = k - 1;
            break;
        end
    end
    if p==0
        p=n;
    end

    U1 = U(:,1:p);
    V1 = V(:,1:p);
    E1 = S(1:p,1:p);
    le = eig(E1,U1'*Cxy*V1);
    le=le';
    %计算频率
    num=1;
    for i=1:length(le)
        if abs( abs(le(i))-1 ) <th3
            if real(le(i))~=0
                temp=atan(imag(le(i))/real(le(i)))/(2*pi);
                f(kj,num)=temp;
            elseif imag(z(i))==0
                f(kj,num)=0;
            end
            num=num+1;
        end
    end
    %计算偏差
    ff=f(kj,:);
    k=find(ff==0);
    ff(k)=[];
    m0=numel(ff);
    
    %估计值分配给最佳逼近真实值
    if m0~=0 
       m1=m0;
       er=zeros(m0,4);
       for i=1:m0
           er(i,:)=abs(ff(i)-ws);
       end
   
       er0=reshape(er,1,4*m0);
       [y i]=min(er0);
       if floor(i/m0)<(i/m0)
           k=floor(i/m0)+1;
       else
            k=floor(i/m0);
       end
       j=i-(k-1)*m0;
       bias1(kj,k) =ff(j)-ws(k);
       var1(kj,k) =(ff(j)-ws(k))^2;
       count(k) = count(k) +1;  
       er(:,k)=10000;
       er(j,:)=10000;
       m1=m1-1;
  
       if m1~=0
          er0=reshape(er,1,4*m0);
          [y i]=min(er0);
          if floor(i/m0)<(i/m0)
              k=floor(i/m0)+1;
          else
              k=floor(i/m0);
          end
          j=i-(k-1)*m0;
          bias1(kj,k) =ff(j)-ws(k);
          var1(kj,k) =(ff(j)-ws(k))^2;
          count(k) = count(k) +1;   
          er(:,k)=10000;
          er(j,:)=10000;
          m1=m1-1;
      end
   
     if m1~=0
          er0=reshape(er,1,4*m0);
          [y i]=min(er0);
          if floor(i/m0)<(i/m0)
             k=floor(i/m0)+1;
          else
             k=floor(i/m0);
          end
          j=i-(k-1)*m0;
           bias1(kj,k) =ff(j)-ws(k);
          var1(kj,k) =(ff(j)-ws(k))^2;
          count(k) = count(k) +1;   
          er(:,k)=10000;
          er(j,:)=10000;
         m1=m1-1;
    end
   
    if m1~=0
         er0=reshape(er,1,4*m0);
         [y i]=min(er0);
         if floor(i/m0)<(i/m0)
              k=floor(i/m0)+1;
         else
              k=floor(i/m0);
         end
         j=i-(k-1)*m0;
         bias1(kj,k) =ff(j)-ws(k);
         var1(kj,k) =(ff(j)-ws(k))^2;
         count(k) = count(k) +1;  
         er(:,k)=10000;
         er(j,:)=10000;
         m1=m1-1;
    end
  end
end

%返回值
bias1=bias1';
var1=var1';
for k = 1 : 4
    if count(k)~=0
    bias(k) = sum(bias1(k,:))/ count(k);
    var(k)  = sum(var1(k,:)) / count(k);
    else
    bias(k) = 0;
    var(k)  = 0;
    end
end
fnum = 100 - count;

bias=bias;
var=var;