Capon-ESPRIT 算法实现：多信源 DOA 估计及功率谱分析

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 = 5;                 %信源数
theta =[10 20 30 40 50];   %波达方向
snr = [30, 25, 20, 15, 5]; %信噪比
n=500;                     %采样数（快拍）
A=exp(-1i*twpi*d.'*sin(theta*derad));  %方向矢量
S=randn(iwave,n );                     %信源信号

for isnr=1:5
    X0=A*S;                              %接收信号
    X=awgn(X0,snr(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=(capon_esprit(dd,Rxx,iwave)); %调用Capon-ESPRIT算法
    doaes(isnr,:)=sort(estimates(1,:));
end

disp(doaes);

figure(1)
polarplot(doaes(1,1)*pi/180,1,'*',doaes(1,2)*pi/180,1,'square',doaes(1,3)*pi/180,1,'d',doaes(1,4)*pi/180,1,'o',doaes(1,5)*pi/180,1,'+');
grid on;
title('Capon-SNR=30');
hold on;
drawnow;

function estimate = capon_esprit(dd,cr,Le)
twpi =2.0*pi;
derad = pi / 180.0;
radeg = 180.0 / pi;
%计算协方差矩阵的逆矩阵
invRxx = inv(cr);
%计算权矩阵
W = invRxx / trace(invRxx);
%对权矩阵进行特征值分解
[V,D]=eig(W);
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:Le,1:Le) EV(2:Le+1,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))*T';
powe = abs(diag(powe).')/Le;
estimate(2,:)=powe;
end

该代码实现了 Capon-ESPRIT 算法，并通过改变信噪比对 DOA 估计结果进行分析。代码中包含以下步骤：

初始化参数：定义阵元数、阵元间距、信源数、波达方向、信噪比、采样数等参数。
生成方向矢量：根据阵元位置和波达方向计算方向矢量。
生成信源信号：生成随机信源信号。
模拟接收信号：根据方向矢量和信源信号，模拟接收信号，并添加噪声。
计算协方差矩阵：根据接收信号计算协方差矩阵。
特征值分解：对协方差矩阵进行特征值分解，并对特征值和特征矢量进行排序。
调用 Capon-ESPRIT 算法：使用 Capon-ESPRIT 算法估计 DOA 和功率谱。
可视化分析：将不同信噪比下的 DOA 估计结果绘制在极坐标图上。

在代码中，capon_esprit 函数实现了 Capon-ESPRIT 算法，主要步骤如下：

计算协方差矩阵的逆矩阵：计算接收信号协方差矩阵的逆矩阵。
计算权矩阵：使用协方差矩阵的逆矩阵计算权矩阵。
特征值分解：对权矩阵进行特征值分解，并对特征值和特征矢量进行排序。
构造矩阵：构造 Exy 和 E_xys 矩阵。
特征值分解：对 E_xys 矩阵进行特征值分解，并对特征值和特征矢量进行排序。
分解矩阵：将 EV_xys 矩阵分解成 Gx 和 Gy 矩阵。
计算 Psi：计算 Psi 矩阵。
特征值分解：对 Psi 矩阵进行特征值分解，并对特征值和特征矢量进行排序。
估计 DOA：根据特征值计算 DOA。
功率估计：计算每个信源的功率谱。

该代码可以帮助理解和应用 Capon-ESPRIT 算法，用于多信源 DOA 估计和功率谱分析。可以通过修改参数来探索不同场景下的算法性能。