波束形成算法

%%%% MUSIC 蒙特卡罗实验%%%%%%%%%% clc;clear;close all; array_num=10; d=0.5;x=1; position=[0:(d/x):(array_num-1)*(d/x)].'; source_num=4;%%？？信源数 doa_ang=[0 -20 30 60];%%？？真实到达角,需为升序,单位：角度 A=DOA_mat(position,doa_ang);%%波达方向矩阵 f=[11 22 33 44];%%信号频率 fs=100; %%采样频率 % ai=1; % for sn=10:10:110 snap_num=256; %%%% 快拍数 s=exp(j*2*pi*f'/fs*[0:snap_num-1]); SNRexp=20:2:20; SNRI1=20; SNRI2=20; SNRI3=20; SNR_num=length(SNRexp); exp_num=1; error_squ_mean=zeros(SNR_num,source_num); for SNR_no=1:SNR_num; doa_est=zeros(exp_num,source_num); for k=1:exp_num %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Signal received generator %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%考虑通道随机误差%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Channel_error_considered=0; %%参数，可以修改它加上不同的通道误差 if Channel_error_considered amp_dis=2; %%参数，表示通道间随机幅度误差，不计固定误差,单位分贝,这里为正负0.5分贝 pha_dis=6; %%参数，表示通道间随机相位误差，不计固定误差，单位度 amp_error=1+(10^(amp_dis/20)-1)*(2*rand(1,array_num)-1); pha_error=exp(j*pha_dis*(2*rand(1,array_num)-1)*pi/180); error_matrix=diag(amp_error)*diag(pha_error); else error_matrix=eye(array_num); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%consider signal mode %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% signal_mode=0;%narrow_band signal or single_frequency signal envelope=ones(source_num,snap_num); if signal_mode%form narrow_band signal N1=[16 32 ]; %%？？信号包络宽度 for i=1:source_num envelope(i,:)=reshape(ones(N1(i),1)*(2*randint(1,snap_num/N1(i))-1),1,snap_num); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%control the power of the interferences%%%%%%%%%%%%%%%%%%%%%%%%%%% SNR_t=[SNRexp(SNR_no) SNRI1 SNRI2 SNRI3]; %%??信噪比 amp=10.^(SNR_t/10); noise_consider=1; color_noise_consider=0; if noise_consider Noise=sqrt(0.5)*(randn(array_num,snap_num)+j*randn(array_num,snap_num)); %%产生复高斯白噪声，噪声功率为1 elseif color_noise_consider Guass_noise=sqrt(0.5)*(randn(array_num,snap_num)+j*randn(array_num,snap_num)); Rn=diag(10*rand(1,array_num)); Rone=0.1*eye(array_num/2); Rone1=0.1*eye(array_num); for i1=1:array_num/2 for k1=1:array_num/2 Rn1(i1,k1)=(0.6^(abs(i1-k1)))*exp(j*0.77*pi*(i1-k1)); end end for i1=1:array_num for k1=1:array_num Rn2(i1,k1)=(0.9^(abs(i1-k1)))*exp(j*0.77*pi*(i1-k1)); end end Rn1=Rone+Rn1; Rn2=Rone1+Rn2; Rz=zeros(size(Rn1,1),size(Rn1,2)); Rcolor_noise1=sqrtm([Rn1 Rz;Rz Rn1])*Guass_noise; Noise=sqrtm(Rn2)*Guass_noise; else Noise=zeros(array_num,snap_num); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% S=envelope.*(diag(amp)*s);% %%%%%%%%%%%%%%%%%%%%%%%%%建立信号模型%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X=error_matrix*A*S+Noise; Rxx=X*X'/snap_num;%+Rcolor_noise; ang_range=[-90:0.1:90]; % Pmusic=my_music(Rxx,source_num,ang_range); [Wopt1,Gtheta1]=Diagloading_BF(Rxx,3,ang_range); % % % [Wopt2,Gtheta2]=Robust_Capon_LiuHongqing(Rxx,0,ang_range,source_num); % % [Wopt3,Gtheta3]=RMVB(Rxx,3,ang_range); [Wopt2,Gtheta2]=Eigenspace_Capon(Rxx,3,ang_range,source_num); [Wopt0,Gtheta0,DL(SNR_no),Spower1]=MY_RMVB(Rxx,3,ang_range,8); [Wopt3,Gtheta3,DL1(SNR_no),Spower2]=MY_RMVB_Double(Rxx,3,ang_range,3); % [Wopt0,Gtheta0,DL1,Spower]=MY_RMVB1(Rxx,3,ang_range,5); % % g=[0;1;0]; % % [Wopt,Gtheta]=LMVC_BF(Rxx,[-23 3 23],g,position,ang_range); % % [Wopt6,Gtheta6,Wopt_LD,Gtheta_LD,Wopt_GSC,Gtheta_GSC,doa_ang,itration,WW]=MY_LMVC(Rxx,[-17 3 24],g,position,ang_range); % % [Wopt7,Gtheta7]=RBF_MSWF(X,3,ang_range,20); % % [Wopt8,Gtheta8]=RBF_MSWF(X,3,ang_range,0); % [Wopt8,Gtheta8]=RBF_MSWF_PA(X,3,ang_range); % % [Wopt8,Gtheta8,DL1,Spower]=VDL_MSWF(X,3,ang_range,1000,8); % % [Wopt9,Gtheta9,DL2,Spower2]=VDL1_MSWF(X,3,ang_range,1000,7); % [Wopt11,Gtheta11,DL11(SNR_no),Spower6]=Worst_case_RMVB(Rxx,3,ang_range,2); % [Wopt0,Gtheta0,DL0(SNR_no),Spower0]=Worst_case_RMVB_Double(Rxx,1,ang_range,1/20); [Wopt4,Gtheta4,DL10(SNR_no),Spower1]=Worst_case_RMVB_Double(Rxx,3,ang_range,1/10); % [Wopt2,Gtheta11,DL11(SNR_no),Spower2]=Worst_case_RMVB_Double(Rxx,1,ang_range,1/2); % [Wopt3,Gtheta12,DL12(SNR_no),Spower3]=Worst_case_RMVB_Double(Rxx,1,ang_range,3/4); % [Wopt4,Gtheta13,DL13(SNR_no),Spower4]=Worst_case_RMVB_Double(Rxx,1,ang_range,1); % [Wopt11,Gtheta11,DL11,Spower6,ee1]=Worst_case_RMVB_ee(Rxx,3,ang_range,1); % [Wopt10,Gtheta10,DL3,Spower3]=VDL1_bi_lanczos(X,2,ang_range,1000,4); % [Wopt11,Gtheta11,DL4,Spower4]=VDL1_Tri_lanczos(X,2,ang_range,1000,4); % [Wopt8,Wopt9,Gtheta8,Gtheta9]=RBFNN_beamformer(Rxx,3,ang_range); % [Wopt8,Gtheta8]=LS_RMVB(Rxx,5,ang_range); figure(2);semilogy(ang_range,abs(Gtheta1),'k-',ang_range,abs(Gtheta4),'r--',ang_range,abs(Gtheta0)/max(abs(Gtheta0)),'k-.',ang_range,abs(Gtheta2),'k:');xlabel('Angle(degree)');ylabel('Array Beampattern(dB)'); legend('Loading-based beamformer','Proposed beamformer','RCB','Eigenspace-based beamformer'); Rs=amp(1)^2; Rin=A(:,2:source_num)*diag(amp(2:source_num))^2*A(:,2:source_num)'+eye(array_num); SINRopt1(k)=10*log10(abs(Rs*A(:,1)'*inv(Rin)*A(:,1))); SINR01(k)=10*log10(abs(Rs*(Wopt0'*A(:,1))^2/(Wopt0'*Rin*Wopt0))); SINR11(k)=10*log10(abs(Rs*(Wopt1'*A(:,1))^2/(Wopt1'*Rin*Wopt1))); SINR21(k)=10*log10(abs(Rs*(Wopt2'*A(:,1))^2/(Wopt2'*Rin*Wopt2))); SINR31(k)=10*log10(abs(Rs*(Wopt3'*A(:,1))^2/(Wopt3'*Rin*Wopt3))); SINR41(k)=10*log10(abs(Rs*(Wopt4'*A(:,1))^2/(Wopt4'*Rin*Wopt4))); end % end SINRopt(SNR_no)=sum(SINRopt1)/exp_num; SINR0(SNR_no)=sum(SINR01)/exp_num; SINR1(SNR_no)=sum(SINR11)/exp_num; SINR2(SNR_no)=sum(SINR21)/exp_num; SINR3(SNR_no)=sum(SINR31)/exp_num; SINR4(SNR_no)=sum(SINR41)/exp_num; % SINRopt(ai)=sum(SINRopt1)/exp_num; % SINR1(ai)=sum(SINR11)/exp_num; % SINR2(ai)=sum(SINR21)/exp_num; % SINR3(ai)=sum(SINR31)/exp_num; % SINR4(ai)=sum(SINR41)/exp_num; % ai=ai+1; SNRexp(SNR_no) end SNR=SNRexp; figure(1); plot(SNR,SINR0,'k-o',SNR,SINR1,'k-.',SNR,SINR2,'k-*',SNR,SINR3,'k-',SNR,SINR4,'k--',SNR,SINRopt,'k:');xlabel('SNR(dB)');ylabel('Output SINR(dB)'); legend('RCB','Loading-based beamformer','Eigenspace-based beamformer','DCRCB','Proposed beamformer','opt'); % aii=1:1:ai-1; % plot(aii,SINR1,'k-.',aii,SINR2,'k-',aii,SINR3,'k-',aii,SINR4,'k--',aii,SINRopt,'k:');xlabel('snapshots');ylabel('Output SINR(dB)'); % legend('Scalar Loading-based GSC beamformer','Proposed beamformer','GSC beamformer','Eigenspace-based beamformer','opt');