近场源定位---论文复现（RD-MUSIC）

weixin_44404510

已于 2023-10-24 09:33:14 修改

阅读量510

点赞数 3

文章标签：算法 1024程序员节

于 2023-08-29 21:49:14 首次发布

本文链接：https://blog.csdn.net/weixin_44404510/article/details/132570622

版权

论文名称：Localization of Near-Field Sources: A Reduced-Dimension MUSIC Algorithm

函数主要由以下构成：

1、信号生成（子函数）

function [Y] = singal(M_half,M,K,N_source,doa_real,dis_real,bochang,d,SNR)
wk=-2*pi*d*sin(doa_real*pi/180)/bochang;
ok=(pi*d^2*cos(doa_real*pi/180).^2)./(bochang*dis_real);
for idd=1:N_source
    for id=-M_half:M_half
        A(id+M_half+1,idd)=exp(1i*(id*wk(idd)+id^2*ok(idd)));
    end
end
Vj=10.^(SNR/10);
S=sqrt(Vj/2)*(randn(N_source,K)+1i*randn(N_source,K));
noise=sqrt(1/2)*(randn(M,K)+1i*randn(M,K));
Y=A*S+noise;
end

2、算法RD-MUSIC（子函数）

function [doa_est,dis_est]=RD_music(M_half,M,Y,K,N_source,d,bochang)
N_MAX=50;
N_doa=N_MAX*2;
Y_YY=Y*Y'/K;
[EV,D]=eig(Y_YY);
[~,I]=sort(diag(D)');
EV=EV(:,I); % 前 MM-N_source【噪声子空间】
U_n=EV(:,1:M-N_source)*EV(:,1:M-N_source)';
e1=zeros(M_half+1,1);
e1(M_half+1,1)=1;

T1=[1:1:M_half+1];
T2=fliplr(T1(1:M_half));
T(1:M_half+1)=T1;
T(M_half+2:M)=T2;

Grid_doa=[-N_MAX:(2*N_MAX)/N_doa:N_MAX-(2*N_MAX)/N_doa];

%     Grid_doa(find(Grid_doa==0))=[];
F_r=zeros(M,M_half+1);
rk=-2*pi*d*sin(Grid_doa*pi/180)/bochang;
for it=1:N_doa
    for id=-M_half:M_half
        F_r(id+M_half+1,T(id+M_half+1))=exp(1i*id*rk(it));
    end
    Q_r=F_r'*U_n*F_r;
    Pm_doa(it)=abs(e1'/(Q_r)*e1);
end

doa_estimaton=findmax(Grid_doa,Pm_doa,N_source,100/N_doa);
rk_estimation=-2*pi*d*sin(doa_est*pi/180)/bochang;
for it=1:N_source
    for id=-M_half:M_half
        F_r(id+M_half+1,T(id+M_half+1))=exp(1i*id*rk_estimation(it));
    end
    Q_r=F_r'*U_n*F_r;
    Bk(:,it)= (Q_r\e1/(e1'/Q_r*e1));
end
q=[0:M_half].';
q=q.^2;
p=ones(M_half+1,1);
p=[p,q];

gk=flip(angle(Bk));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%复指数修正%%%%%%%%%%%%%%
if max(max(gk<-0.00001))
    for ip=1:N_source
        imm=1;
        while imm<=M_half
            if gk(imm+1,ip)>gk(imm,ip)&&(gk(imm+1,ip)-gk(imm,ip))>=gk(2,ip)*(imm-1)
                imm=imm+1;
            else
                gk(imm+1,ip)=gk(imm+1,ip)+2*pi;
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CK=real((p.'*p)\p.'*gk);
dis_est=(pi*d^2*cos(doa_est*pi/180).^2)./(bochang*CK(2,:));
end

3、findmax（子函数）

function theta=findmax(areas,Pm,N,resolve1)
zz=[];

Pm=Pm(:);
d_Pm=diff(Pm);
d_Pm=[-1;d_Pm;1];

dpos=find(d_Pm>=0);
jj=1;
for ii=1:length(dpos)-1
    if d_Pm(dpos(ii)+1)<=0
      zz(jj)= dpos(ii);
      jj=jj+1;
    end
end

theta=zeros(N,1);

if length(zz)<N
    theta=zeros(N,1);
elseif  length(zz)==N
    theta(1:length(zz))=areas(zz); 
else
    [~,I]=sort(Pm(zz), 'descend');    
    temp_alpha=areas(zz(I));

    theta=[];
    theta(1)=temp_alpha(1);
    ii=2;
    while length(theta)~=N
        if min(abs(theta- temp_alpha(ii)))> resolve1
            theta=[theta; temp_alpha(ii)];
        end
        ii=ii+1;
    end

end

[theta,~]=sort(theta);  
theta=theta';

4、MTKL随SNR（主函数）

clear all
clc;

SNR_t=[-5:5:20];
K_t=[20 100 200 300 400 500];
H=SNR_t;
N_source=2;
runtimes=200;
lambda=1;
d=1/4*lambda;
M=11;
M_half=floor(M/2);

%%%%%%%%%%%%%%%%%%%%%   随SNR的MTKL    %%%%%%%%%%%%%%%
for ip=1:length(H)
    SNR=SNR_t(ip);
    K=500; 
for it=1:runtimes  
    [ ip  it  ]
     doa_real=[ 30 -20] ;
    dis_real=[ 0.2  0.6]*lambda;

   
 %%%%%%%%%%%%%%%% 信号生成%%%%%%%%%%%%%%%
 dL=(-(M-1)/2:(M-1)/2)'*d;
 Y=singal(M_half,M,K,N_source,doa_real,dis_real,lambda,d,SNR);
 

%%%%%%%%%%%%%%% 算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [doa_est,dis_est]=RD_music(M_half,M,Y,K,N_source,d,lambda);
  
    
    RMSE_RD_music_doa_iter(it,:)=doa_real-doa_est;
    RMSE_RD_music_dis_iter(it,:)=dis_real-dis_est;

end

    RMSE_RD_music_doa(ip)=norm(RMSE_RD_music_doa_iter,'f')/sqrt(runtimes*N_source);
    RMSE_RD_music_dis(ip)=norm(RMSE_RD_music_dis_iter,'f')/sqrt(runtimes*N_source);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    figure(1);

     hold on;
   grid on; 

set(gca,'yscale','log');
plot(H,RMSE_RD_music_doa,'color',[138,185,226]/255,'LineWidth',1.1,'Marker','s');
l1=xlabel('SNR');
 l2=ylabel('RMSE-DOA');
    set(l1,'Fontname','Times New Roman') ；
set(l2,'Fontname','Times New Roman')  ；  



    figure(2);
    hold on;
    grid on;  
set(gca,'yscale','log');
plot(H,RMSE_RD_music_dis,'color',[138,185,226]/255,'LineWidth',1.1,'Marker','s');

l1=xlabel('SNR');
l2=ylabel('RMSE-Distance');
   set(l1,'Fontname','Times New Roman') % %
 set(l2,'Fontname','Times New Roman') % % % %    

  
  save('SNR_3source')

就到这，如果实现有问题，随时联系。（侵权联系就删，仅供学习参考适用）