• 日期: 2025-07-23
  • 标签: 常规
clear all
close all
derad = pi/180;            %角度->弧度
radeg = 180/pi;            %弧度->角度
twpi=2*pi;
kelm = 16;                  %阵元数
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;                          %计算协方差矩阵
    doaes_tls(isnr,:)=sort(tls_esprit(dd,Rxx,iwave)); %调用TLS-ESPRIT子程序估计DOA
    doaes_capon(isnr,:)=sort(capon(Rxx,iwave));  %调用Capon子程序估计DOA
end
disp(doaes_tls);
disp(doaes_capon);

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

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

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;
end

function estimate = capon(Rxx,Le)
[N,~]=size(Rxx);
J=eye(N)-ones(N)/N;
R=J*Rxx*J/N;
[V,D]=eig(R);
EVA=real(diag(D)');
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(V(:,I));
EVA=EVA(1:Le);
EV=EV(:,1:Le);
w=pinv(EV)*ones(Le,1)/(ones(Le,1)'*pinv(EV)*ones(Le,1));
doa=asin(angle(w)/pi)*180/pi;
estimate=doa;
end

这段代码中使用了TLS-ESPRIT算法进行信号方向估计，并利用Capon算法进行信号方向估计，并绘制了两种算法的极坐标图，方便对比分析两种算法的性能。

  1. 在主程序中添加调用Capon算法的代码：
for isnr=1:5
    X0=A*S;                              %接收信号
    X=awgn(X0,snr(isnr),'measured') ;     %添加噪声
    Rxx=X*X'/n;                          %计算协方差矩阵
    doaes_tls(isnr,:)=sort(tls_esprit(dd,Rxx,iwave)); %调用TLS-ESPRIT子程序估计DOA
    doaes_capon(isnr,:)=sort(capon(Rxx,iwave));  %调用Capon子程序估计DOA
end
disp(doaes_tls);
disp(doaes_capon);
  1. 添加Capon算法的子程序：
function estimate = capon(Rxx,Le)
[N,~]=size(Rxx);
J=eye(N)-ones(N)/N;
R=J*Rxx*J/N;
[V,D]=eig(R);
EVA=real(diag(D)');
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(V(:,I));
EVA=EVA(1:Le);
EV=EV(:,1:Le);
w=pinv(EV)*ones(Le,1)/(ones(Le,1)'*pinv(EV)*ones(Le,1));
doa=asin(angle(w)/pi)*180/pi;
estimate=doa;
end
  1. 在主程序中添加绘制Capon算法的极坐标图的代码：
figure(2)
polarplot(doaes_capon(1,1)*pi/180,1,'*',doaes_capon(1,2)*pi/180,1,'square',doaes_capon(1,3)*pi/180,1,'d',doaes_capon(1,4)*pi/180,1,'o',doaes_capon(1,5)*pi/180,1,'+');
grid on;
title('Capon-SNR=30');
hold on;
drawnow;

完整代码如下：

clear all
close all
derad = pi/180;            %角度->弧度
radeg = 180/pi;            %弧度->角度
twpi=2*pi;
kelm = 16;                  %阵元数
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;                          %计算协方差矩阵
    doaes_tls(isnr,:)=sort(tls_esprit(dd,Rxx,iwave)); %调用TLS-ESPRIT子程序估计DOA
    doaes_capon(isnr,:)=sort(capon(Rxx,iwave));  %调用Capon子程序估计DOA
end
disp(doaes_tls);
disp(doaes_capon);

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

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

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;
end

function estimate = capon(Rxx,Le)
[N,~]=size(Rxx);
J=eye(N)-ones(N)/N;
R=J*Rxx*J/N;
[V,D]=eig(R);
EVA=real(diag(D)');
[EVA,I]=sort(EVA);
EVA=fliplr(EVA);
EV=fliplr(V(:,I));
EVA=EVA(1:Le);
EV=EV(:,1:Le);
w=pinv(EV)*ones(Le,1)/(ones(Le,1)'*pinv(EV)*ones(Le,1));
doa=asin(angle(w)/pi)*180/pi;
estimate=doa;
end
TLS-ESPRIT 和 Capon 算法方向估计对比

原文地址: https://www.cveoy.top/t/topic/nH6G 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录