esprit_tls.m

张贤达老师 通信信号处理4－9章的程序 这本书是张老师的经典之作。

%%%%%%%%%%%%%%%%%%%%%%%  simulatation of the Pisarenko method          %%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%  using ESPRIT algorithms  P156:3.18       %%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%          initiatation        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
n=1:1280;                        %the length of data,it is larger than N,because of Rx
number=1000;                     %the iteration number
for pp=1:number
v1=sqrt(20)*sin(2*pi*0.2*n);   
v2=sqrt(2)*sin(2*pi*0.213*n);
w=randn(1,10000);
x=v1+v2+w(1:length(n));                %generate obervated data, which is in added gaussian white noise

m=30;                       %the number of the receivers
p=4;
rx=xcorr(x,m,'unbiased');
for i=1:m
    Rxx(i,:)=rx(m+2-i:2*m+1-i);
end
for i=1:m
    Rxy(i,:)=rx(m+3-i:2*(m+1)-i);
end

%%%%%%%%%%%%%%%%%%% the basic method %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for i=1:m
%     for j=1:m
%         if (i-j==1)
%             Z(i,j)=1;
%         else
%             Z(i,j)=0;
%         end
%     end
% end


eigenvalue=eig(Rxx);
%Cxx=Rxx-mean(eigenvalue(p+1:m));
%Cxy=Rxy-mean(eigenvalue(p+1:m))*Z;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%% the improved method %%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma=mean(eigenvalue(p+1:m)); 
[U,S]=eig(Rxx);
 for i=1:m                     % 经过特征值分解后，令后m-p个对角元素的值的平均为噪声方差，构造新的对角矩阵A
     if (i<=p)                 % 其前p个对角元素为S中的前p个对角元素减去噪声方差，后m-p个元素直接置零
         A(i,i)=S(i,i)-sigma;
     else
         A(i,i)=0;
     end
 end
 Cxx=U*A*U';                         % not use the mean of inferior eigenvalue, but use the left eig matrix
%%% 这样就得到了新的关于Cxx的计算公式

% 使用类似的方法构造出Cxy矩阵。由于sigma*Z表示的是E（w（n）w（n+1）’），所以可以首先构造对角阵，使其前
%% p个元素为平均值sigma，后m-p个元素为上面得到的S的后m-p个元素，然后分别用U左乘和U’右乘该对角阵得到E（w（n）w（n）’）
%% 的估计，再通过适当拼接（补零或者删去某一行）来得到E（w（n）w（n+1）’）的估计
 for i=1:m                    
     if (i<=p)
         k(i)=sigma;
     else
         k(i)=eigenvalue(i);
     end
 end
 B=U*diag(k)*U';
 B=[B(:,2:m) zeros(m,1)];   
 Cxy=Rxy-B;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U,S,V]=svd(Cxx);
B=U(:,1:p)'*Cxy*V(:,1:p);
A=S(1:p,1:p);
d=eig(A,B);
d_nonzero=0;                        % nonzero eigenvalue
for i=1:length(d)
    if (abs(d(i))>=0.9 & abs(d(i)<=1.1))
        d_nonzero=[d_nonzero d(i)];
    end
end
frequency=angle(d_nonzero(2:length(d_nonzero)))/2/pi;
frequency=sort(frequency);            % order the turn
for i=1:length(frequency)
    if (frequency(i)>0.19)
        k=i;
        break
    end
end
result(pp,:)=frequency(k:k+1);       % each of its row is the result of one run
end