ls_esprit.m

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

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

th1=0.01;
th2=0.05;
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);
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');
    for i=1:n
       if d(i)/d(1)<th1
            p=i-1;%确定噪声
            break;
       end
    end
    noise=sum(d(p+1:n))/(n-p); %得到噪声功率
    
    %计算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;
    
    %广义特征值分解
    le= eig(Cxx,Cxy);
    le=le';
    %计算频率
    num=1;
    for i=1:length(le)
        if abs( abs(le(i))-1 ) <th2
            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;