阵列信号处理——旋转不变技术（ESPRIT）

ESPRIT是借助旋转不变技术估计信号参数(estimating signal parameter via rotational invariance techniques)方法的英文缩写，最早由Roy等人[233]于1986年提出。ESPRIT方法现已成为现代信号处理中一种代表性方法，并得到了广泛的应用。
和MUSIC方法存在阵元空间和波束空间两类算法一样，ESPRIT方法也有阵元空间E-SPRIT算法和波束空间ESPRIT算法两种类型，并且还有专门针对复观测数据的酉ESPRIT方法。

算法一、基本ESPRIT算法步骤

步骤1、利用已知的观测数据x(1),…..,x(N)估计自协方差函数.

步骤2、由估计的自协方差函数构造MxM自协方差矩阵Rxx和MxM互协方差Rxy。

步骤3、求Rxx的特征值分解。对于M>p,最小特征值的平均值作为噪声方差的估计。

步骤4、利用计算Cxx=Rxx-I和Cxy=Rxy-Z。

步骤5、求矩阵束{Cxx，Cxy}的广义特征值分解，得到位于单位圆上的p个广义特征值，它们直接给出谐波频率wi的估计。然后，利用估计波
达方向θ1,…，θp。

Roy和Kailath指出[232],最小二乘算子会导致在求解广义特征值问题的某些潜在的数值困难。现在已广泛认识到，奇异值分解（SVD）和总体最小二乘（TLS）的应用可以将一个较大维数（MxM）病态广义特征问题转化为一个较小维数（p×p)的无病态广义特征问题。

算法二、TLS-ESPRIT算法步骤

TLS-ESPRIT算法matlab仿真实验

TLS-ESPRIT函数程序如下：

function estimate = tls_esprit(dd,cr,Le)
twpi =2.0*pi;
derad = pi / 180.0;
radeg = 180.0 / pi;
%对接收信号协方差矩阵进行特征值分解
[K,KK] = size(cr);
[V,D]=eig(cr);
EVA = real(diag(D)');
[EVA,I] = sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(V(:,I));
%构造E_{xy}和E_xys=E_{xy}HE_{xy}
Exy =[EV(1:K-1,1:Le) EV(2:K,1:Le)];
E_xys = Exy'*Exy;
%对E_xys进行特征值分解
[V,D]=eig(E_xys);
EVA_xys=real(diag(D)');
[EVA_xys,I] =sort(EVA_xys);
EVA_xys=fliplr(EVA_xys);
EV_xys=fliplr(V(:,I));
%将EV_xys分解
Gx = EV_xys(1:Le,Le+1:Le*2);
Gy=EV_xys(Le+1:Le*2,Le+1:Le*2);
%计算Psi=-Gx[Gy]{-1}
Psi = -Gx/Gy;
%对Psi进行特征值分解
[V,D]=eig(Psi);
EGS = diag(D).';
[EGS,I] = sort(EGS);
EGS=fliplr(EGS);
EVS=fliplr(V(:,I));
%估计DOA
ephi = atan2(imag(EGS), real(EGS));
ange = -asin( ephi / twpi / dd ) * radeg;
estimate(1,:)=ange;
%功率估计
T=inv(EVS);
powe = T*diag(EVA(1:Le)-EVA(K))*T';
powe = abs(diag(powe).')/K;
estimate(2,:)=powe;

主函数程序如下：

clear all
close all
derad = pi/180;            %角度->弧度
radeg = 180/pi;            %弧度->角度
twpi=2*pi;
kelm = 8;                  %阵元数
dd=0.5;                    %阵元间距
d=0:dd:(kelm-1)*dd;        
iwave = 3;                 %信源数
theta =[10 20 30];         %波达方向
snr = 10;                  %信噪比
n=500;                     %采样数（快拍）
A=exp(-1i*twpi*d.'*sin(theta*derad));  %方向矢量
S=randn(iwave,n );                     %信源信号
snr0=0:3:100 ;                         %信噪比
for isnr=1:10
    X0=A*S;                            %接收信号
    X=awgn(X0,snr0(isnr),'measured') ; %添加噪声
    Rxx=X*X'/n;                        %计算协方差矩阵
    [EV,D]=eig(Rxx);                   %特征值分解
    EVA=diag(D)';
    [EVA,I]=sort(EVA);                 %特征值从小到大排序
    EVA=fliplr(EVA);                   %左右翻转，从大到小排序
    EV=fliplr(EV(:,I));                %对应特征矢量排序
    estimates=(tls_esprit(dd,Rxx,iwave)); %调用子程序
    doaes(isnr,:)=sort(estimates(1,:));
end
disp(doaes);

运行结果如下：

参考文献：

阵列信号处理及MATLAB实现；张小飞，陈华伟，仇小锋（编著） 

现代信号处理（第三版）；张贤达（编著）