Capon 算法信号方向估计 MATLAB 代码

clear all close all derad = pi/180; %角度->弧度 radeg = 180/pi; %弧度->角度 twpi=2pi; 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(-1itwpid.'sin(thetaderad)); %方向矢量 S=randn(iwave,n ); %信源信号

for isnr=1:5 X0=AS; %接收信号 X=awgn(X0,snr(isnr),'measured') ; %添加噪声 Rxx=XX'/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.0pi; 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:Le2); Gy=EV_xys(Le+1:Le2,Le+1:Le2); %计算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 = Tdiag(EVA(1:Le))*T'; powe = abs(diag(powe).')/K; estimate(2,:)=powe; end