首页 > 代码库 > DOA——MUSIC算法
DOA——MUSIC算法
【未完待续】
一、MUSIC
1、
clear all;%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%一个符号对应的比特数Fbit=M*Fsym;%比特速率Nsour=3;%信源数Angle=[5,15,35];%信源的来波方向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因子Del=5;%群延迟因子% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相干信源--信源1、信源2与信源3Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsamp=1024;%采样点数或者快拍数 i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[2,2,2];%三信源的信噪比% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%子阵个数L=Ntx-p+1;%子阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产生背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).‘];%接收信号的采样点数%----单信号源s1=[Rcos1(t).‘;Rcos2(t).‘;Rcos3(t).‘];%矩阵维数=信源数*抽样点数ps=diag((s1*s1‘)/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta‘)*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR比条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n‘*n);for h=1:Nsamp pn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]‘;%矩阵维数=阵元数*1tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X‘)/Nsamp;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sub_FRxx=zeros(L,L);sub_BRxx=zeros(L,L);for i=1:p sub_FR=zeros(L,Nsamp); sub_BR=zeros(L,Nsamp); sub_FR=X(i:1:i+L-1,:); last=Ntx+1-i; first=last-L+1; sub_BR=conj(X(last:-1:first,:)); sub_FRxx=sub_FRxx+((sub_FR*sub_FR‘)./Nsamp); sub_BRxx=sub_BRxx+((sub_BR*sub_BR‘)./Nsamp);endsub_FRxx=sub_FRxx./p;sub_BRxx=sub_BRxx./p;sub_Rxx=(sub_FRxx+sub_BRxx)./2;[VFB,HFB]=eig(sub_Rxx);[HFB,IFB]=sort(diag(HFB),1);VFB=VFB(:,IFB);VnFB=VFB(:,1:L-Nsour);VsFB=VFB(:,L-Nsour+1:L);%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ScanAng=[-90:1:90];for i=1:length(ScanAng) tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda; tmp3=tmp2*tmp1; tmp5=tmp2*tmp4; A_Sita=exp(tmp3); Sub_Sita=exp(tmp5); Sub_FBsita(i)=(Sub_Sita‘*Sub_Sita)/(Sub_Sita‘*VnFB*VnFB‘*Sub_Sita);end figure(1);semilogy(ScanAng,real(Sub_FBsita),‘bo-‘);axis([-60 60 0.1 1e7]);xlabel(‘M_Angle(deg)‘);ylabel(‘M_Spectrum‘);grid on
2、
clear all;%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%一个符号对应的比特数Fbit=M*Fsym;%比特速率Nsour=3;%信源数Angle=[10,40,80];%信源的来波方向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因子Del=5;%群延迟因子% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos1=Rcos11;%构造相干信源--信源1、信源2与信源3Rcos2=Rcos21;%构造相干信源--信源1、信源2与信源3Rcos3=Rcos31;%构造相干信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsamp=1024;%采样点数或者快拍数 i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[10,10,10];%三信源的信噪比% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%子阵个数L=Ntx-p+1;%子阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产生背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).‘];%接收信号的采样点数%----单信号源s1=[Rcos1(t).‘;Rcos2(t).‘;Rcos3(t).‘];%矩阵维数=信源数*抽样点数ps=diag((s1*s1‘)/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta‘)*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR比条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n‘*n);for h=1:Nsamp pn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]‘;%矩阵维数=阵元数*1tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X‘)/Nsamp;[V,H]=eig(Rxx);%MUSIC算法---MUltiSIgnal Classification[H,I]=sort(diag(H),1);%特征值按照升序排列V=V(:,I);%特征值对应的特征向量也按照相应特征值的升序排列Vn=V(:,1:Ntx-Nsour);%噪声子空间---协方差的特征向量--最小特征值对应的特征向量Vs=V(:,Ntx-Nsour+1:Ntx);%信号子空间---协方差的特征向量--最大特征值对应的特征向量ScanAng=[-90:1:90];for i=1:length(ScanAng) tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda; tmp3=tmp2*tmp1; tmp5=tmp2*tmp4; A_Sita=exp(tmp3);MUSIC_Sita(i)=(A_Sita‘*A_Sita)/(A_Sita‘*Vn*Vn‘*A_Sita);end figure(1);semilogy(ScanAng,real(MUSIC_Sita),‘g*-‘);axis([-90 90 0.1 1e7]);xlabel(‘M_Angle(deg)‘);ylabel(‘M_Spectrum‘);grid on
3、
clear all;%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsym=500;%符号个数Fsym=1;%符号速率M=3;%一个符号对应的比特数Fbit=M*Fsym;%比特速率Nsour=3;%信源数Angle=[5,8,35];%信源的来波方向Fc=10;%载波频率Fs=100;%抽样频率R=0.5;%滚降因子Del=5;%群延迟因子% Nsamp=50;%采样点数或者快拍数S1=randint(Nsym,1,2^M);S2=randint(Nsym,1,2^M);S3=randint(Nsym,1,2^M);PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);Rcos11=rcosflt(PM1,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos21=rcosflt(PM2,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos31=rcosflt(PM3,Fsym,Fs,‘fir/sqrt/Fs‘,R,Del);Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相干信源--信源1、信源2与信源3Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3save xyc3 Rcos1 Rcos2 Rcos3%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波Nsamp=512;%采样点数或者快拍数 i=sqrt(-1);j=i;Ntx=8;%阵列数SNR=[2,2,2];%三信源的信噪比% sn=10; %----单信号源Lamda=2;%信号波长D=Lamda/2;%线性阵列的距离p=3;%子阵个数L=Ntx-p+1;%子阵阵元数nr=randn(Ntx,Nsamp);ni=randn(Ntx,Nsamp);n=nr+j*ni;%产生背景噪声load xyc3;t=1:Nsamp;% s1=[Rcos1(t).‘];%接收信号的采样点数%----单信号源s1=[Rcos1(t).‘;Rcos2(t).‘;Rcos3(t).‘];%矩阵维数=信源数*抽样点数ps=diag((s1*s1‘)/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1% delta1=ps./(2*10.^(sn/10)); %----单信号源delta2=diag(delta1);%矩阵维数=1*1delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1Rev_s1=(1./delta‘)*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数%计算各信源SNR比条件下，阵列接收到的信号幅度%Pn=zeros(Nsamp,1);pn=zeros(Ntx,Nsamp);Pn=diag(n‘*n);for h=1:Nsamp pn(:,h)=n(:,h)./sqrt(Pn(h,:));endRev_n=pn;%计算各阵列接收到的背景噪声下的信号幅度%tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----单信号源tmp1=[0:Ntx-1]‘;%矩阵维数=阵元数*1tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1a1=tmp1*tmp;%矩阵维数=阵元数*信源数A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数Rxx=(X*X‘)/Nsamp;Rxx_fb=zeros(L,L);Rxx_f=zeros(L,L);Rxx_b=zeros(L,L);J=fliplr(eye(L));for m=1:p for k=1:p Rxx_f=Rxx_f+Rxx(m:1:m+L-1,k:1:k+L-1)*Rxx(k:1:k+L-1,m:1:m+L-1); Rxx_b=Rxx_b+J*conj(Rxx(m:1:m+L-1,k:1:k+L-1))*conj(Rxx(k:1:k+L-1,m:1:m+L-1))*J; endendRxx_f=Rxx_f./p;Rxx_b=Rxx_b./p;Rxx_fb=(Rxx_f+Rxx_b)./p;[V_fb,H_fb]=eig(Rxx_fb);%特征分解---MUltiSIgnal Classification[H_fb,I_fb]=sort(diag(H_fb),1);%特征值按照升序排列V_fb=V_fb(:,I_fb);%特征值对应的特征向量也按照相应特征值的升序排列Vn_fb=V_fb(:,1:L-Nsour);Vs_fb=V_fb(:,L-Nsour+1:L);ScanAng=[-90:1:90];for i=1:length(ScanAng) tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda; tmp3=tmp2*tmp1; tmp5=tmp2*tmp4; A_Sita=exp(tmp3); Sub_Sita=exp(tmp5); fb_sita(i)=(Sub_Sita‘*Sub_Sita)/(Sub_Sita‘*Vn_fb*Vn_fb‘*Sub_Sita);end figure(1);semilogy(ScanAng,real(fb_sita),‘r-‘);axis([-60 60 0.1 1e7]);xlabel(‘M_Angle(deg)‘);ylabel(‘M_Spectrum‘);grid on
4、
%=========================================================================% UCA_multi_in_2D% %=========================================================================clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; phi = 60;%% ----------------入射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中心频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;dir8 =(62.5:.25:72.5)*pi/180; dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;dir10=(49:.25:59)*pi/180; ang=(50:.25:70)*pi/180;e_dir7 = zeros(1,length(sr_array));e_dir8 = zeros(1,length(sr_array));e_dir9 = zeros(1,length(sr_array));e_dir10= zeros(1,length(sr_array));for ss = 1:length(sr_array) snr = sr_array(ss); %--------------------7阵元--------------------------------------- sensor_num = 7; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; R = namda/(1-cos(d_angle(2))); x = R*cos(d_angle);y = R*sin(d_angle); theta = d_angle(2)*180/pi; for it = 1:iteration n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir7),length(ang)); for i=1:length(dir7) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir7(a)*180/pi; e_dir7(ss) = e_dir7(ss)+(aa-theta)^2; end %--------------------8阵元------------------------------------- sensor_num = 8; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2))); x = R*cos(d_angle);y = R*sin(d_angle); theta = 1.5*d_angle(2)*180/pi; for it = 1:iteration tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir8),length(ang)); for i=1:length(dir8) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir8(a)*180/pi; e_dir8(ss) = e_dir8(ss)+(aa-theta)^2; end %---------------9阵元--------------------------------------------------- sensor_num = 9; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; R = namda/(1-cos(d_angle(2))); x = R*cos(d_angle);y = R*sin(d_angle); theta = d_angle(2)*180/pi; for it = 1:iteration n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir9),length(ang)); for i=1:length(dir9) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir9(a)*180/pi; e_dir9(ss) = e_dir9(ss)+(aa-theta)^2; end %---------------10阵元--------------------------------------------------- sensor_num = 10; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2))); x = R*cos(d_angle);y = R*sin(d_angle); theta = 1.5*d_angle(2)*180/pi; for it = 1:iteration tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir10),length(ang)); for i=1:length(dir10) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir10(a)*180/pi; e_dir10(ss) = e_dir10(ss)+(aa-theta)^2; end endfigure;% subplot(121);% plot(sr_array+45,e_dir10(1,:)/iteration,‘-^k‘,sr_array+45,e_dir9(1,:)/iteration,‘-*k‘);% legend(‘8元UCA‘,‘9元UCA‘);% grid on;%axis([-10,20,-.01,.45]);% colormap gray;% xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);% title(‘半径相同精度试验‘);% % % subplot(122);plot(sr_array+45,e_dir7/iteration,‘-^k‘,sr_array+45,e_dir8/iteration,‘-*k‘);hold on;plot(sr_array+45,e_dir9/iteration,‘-sk‘,sr_array+45,e_dir10/iteration,‘-dk‘);legend(‘7阵元‘,‘8阵元‘,‘9阵元‘,‘10阵元‘);grid on;axis([-5.5,18,-.02,.65]);colormap gray;xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);title(‘不同阵元最大半径测向试验‘);
5、
%=========================================================================% UCA_multi_in_2D% %=========================================================================clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; phi = 60;%% ----------------入射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中心频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));R = 0.01;dir9(1,:)=(-50:.25:-40)*pi/180;dir9(2,:)=(75:.25:85)*pi/180; dir10(1,:)=(-51:.25:-41)*pi/180;dir10(2,:)=(107:.25:117)*pi/180; ang=(50:.25:70)*pi/180;e_dir9= zeros(2,length(sr_array));e_ang9= zeros(2,length(sr_array));e_dir10= zeros(2,length(sr_array));e_ang10= zeros(2,length(sr_array));for ss = 1:length(sr_array) snr = sr_array(ss); %---------------7阵元--------------------------------------------------- sensor_num = 9; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta_array = [-9*d_angle(2)/8,d_angle(3)]*180/pi; for it = 1:iteration n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); for dd = 1:2 theta = theta_array(dd); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir9(dd,:)),length(ang)); for i=1:length(dir9(dd,:)) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir9(dd,i))+y*sin(dir9(dd,i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir9(dd,a)*180/pi; bb = ang(b)*180/pi; e_dir9(dd,ss) = e_dir9(dd,ss)+(aa-theta)^2; e_ang9(dd,ss) = e_ang9(dd,ss)+(bb-phi)^2; end end %---------------8阵元--------------------------------------------------- sensor_num = 8; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta_array = [-9*d_angle(2)/8,2.5*d_angle(2)]*180/pi; for it = 1:iteration for dd = 1:2 theta = theta_array(dd); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir10(dd,:)),length(ang)); for i=1:length(dir10(dd,:)) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir10(dd,i))+y*sin(dir10(dd,i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir10(dd,a)*180/pi; bb = ang(b)*180/pi; e_dir10(dd,ss) = e_dir10(dd,ss)+(aa-theta)^2; e_ang10(dd,ss) = e_ang10(dd,ss)+(bb-phi)^2; end endendfigure;subplot(121);plot(sr_array+45,e_dir9(1,:)/iteration,‘-^k‘,sr_array+45,e_dir9(2,:)/iteration,‘-*k‘);legend(‘-9/8 360/M‘,‘3 360/M‘);grid on;%axis([-10,20,-.01,.45]);colormap gray;xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);title(‘9元阵方向角误差‘);% % subplot(222);% plot(sr_array+45,e_ang9(1,:)/iteration,‘-^k‘,sr_array+45,e_ang9(2,:)/iteration,‘-*k‘);% legend(‘-9/8 360/M‘,‘3 360/M‘);% grid on;axis([-10,20,-.1,1.45]);% colormap gray;% xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);% title(‘7元阵俯仰角误差‘);subplot(122);plot(sr_array+45,e_dir10(1,:)/iteration,‘-^k‘,sr_array+45,e_dir10(2,:)/iteration,‘-*k‘);legend(‘-5/4 360/M‘,‘5/2 360/M‘);grid on;%axis([-10,20,-.01,.45]);colormap gray;xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);title(‘8元阵方向角误差‘);% % % subplot(224);% plot(sr_array+45,e_ang10(1,:)/iteration,‘-^k‘,sr_array+45,e_ang10(2,:)/iteration,‘-*k‘);% legend(‘-5/4 360/M‘,‘3 360/M‘);% grid on;axis([-10,20,-.1,1.45]);% colormap gray;% xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);% title(‘8元阵俯仰角误差‘);
6、
%=========================================================================% UCA_multi_in_2D% %=========================================================================clc;clear all;close all;%------------------------常数表-------------------------------c = 3e8;namda = c/18e9;est_num = 1;iteration = 100;sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; phi = 60;%% ----------------入射信号模型-------------------------------N_x = 2^5; %快拍点数F0 = 18e9; %中心频率B = 20e6; %带宽Fs = 2*B; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/18e9;st = exp(1j*2*pi*(F0*t+.5*u*t.^2));dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;dir8 =(62.5:.25:72.5)*pi/180; dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;dir10=(49:.25:59)*pi/180; ang=(50:.25:70)*pi/180;e_dir7 = zeros(1,length(sr_array));e_dir8 = zeros(1,length(sr_array));e_dir9 = zeros(1,length(sr_array));e_dir10= zeros(1,length(sr_array));R = 0.1;for ss = 1:length(sr_array) snr = sr_array(ss); %--------------------7阵元--------------------------------------- sensor_num = 7; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta = d_angle(2)*180/pi; for it = 1:iteration n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir7),length(ang)); for i=1:length(dir7) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir7(a)*180/pi; e_dir7(ss) = e_dir7(ss)+(aa-theta)^2; end %--------------------8阵元------------------------------------- sensor_num = 8; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta = 1.5*d_angle(2)*180/pi; for it = 1:iteration tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir8),length(ang)); for i=1:length(dir8) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir8(a)*180/pi; e_dir8(ss) = e_dir8(ss)+(aa-theta)^2; end %---------------9阵元--------------------------------------------------- sensor_num = 9; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta = d_angle(2)*180/pi; for it = 1:iteration n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir9),length(ang)); for i=1:length(dir9) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir9(a)*180/pi; e_dir9(ss) = e_dir9(ss)+(aa-theta)^2; end %---------------10阵元--------------------------------------------------- sensor_num = 10; d_angle = (0:sensor_num-1)‘*2*pi/sensor_num; x = R*cos(d_angle);y = R*sin(d_angle); theta = 1.5*d_angle(2)*180/pi; for it = 1:iteration tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180)); A = exp(-1j*2*pi*tao./l); n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10)); xt = A*(sqrt(10^(snr/10))*st)+n; % -------------------2D-MUSIC算法----------------------- Rxx = xt*xt‘/N_x; [U,S] = svd(Rxx); [~,index] = sort(diag(S)); Un = U(:,index(1:sensor_num-est_num)); Gn = Un*Un‘; Pmusic = zeros(length(dir10),length(ang)); for i=1:length(dir10) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic(i,k)=1./abs((a_theta)‘*Gn*a_theta); end end [a,b]=find(Pmusic==max(max(Pmusic))); aa = dir10(a)*180/pi; e_dir10(ss) = e_dir10(ss)+(aa-theta)^2; end endfigure;plot(sr_array+45,e_dir7/iteration,‘-^k‘,sr_array+45,e_dir8/iteration,‘-*k‘);hold on;plot(sr_array+45,e_dir9/iteration,‘-sk‘,sr_array+45,e_dir10/iteration,‘-dk‘);legend(‘7阵元‘,‘8阵元‘,‘9阵元‘,‘10阵元‘);grid on;axis([-5.5,18,-.02,.25]);colormap gray;xlabel(‘信噪比/dB‘);ylabel(‘均方误差/°‘);title(‘不同阵元相同半径测向试验‘);
7、
%=========================================================================% Circular Array Classical-Music% %=========================================================================clc;clear all;close all;c = 3e8;phi = 60;namda = c/18e9;R = 7.5/100;R = 10/100;snr = -35; %信噪比N_x = 2^5; %快拍点数F0 = 18e9; %中心频率B = 20e6; %带宽Fs = 40e6; %采样频率Ts = 1/Fs; %采样时间T = (N_x-1)*Ts; %快拍持续时间u = B/T; %频率变化率t = -T/2:Ts:T/2; %时间轴点l = c/F0;st = sqrt(10^(snr/10))*exp(1j*2*pi*(F0*t+.5*u*t.^2));%% -----------------9------------------------sensor_num = 9;R = 7.1239/100;d_angle = (0:sensor_num-1)‘*2*pi/sensor_num;theta =d_angle(2)*180/pi;x = R*cos(d_angle);y = R*sin(d_angle); tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));A = exp(-1j*2*pi*tao/l);n = randn(sensor_num,N_x)*sqrt(10^(-45/10));xt = A*st+n;% -------------------2D-MUSIC算法-----------------------Rx = xt*xt‘/N_x;[U,S] = eig(Rx);est_sour = 1;[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);Gn = Un*Un‘;dir=(-180:.25:179.8)*pi/180; ang=(20:.25:91)*pi/180;Pmusic9 = zeros(length(dir),length(ang));for i=1:length(dir) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic9(i,k)=1./abs((a_theta)‘*Gn*a_theta); endendP_music9 = 10*log10(Pmusic9/min(min(Pmusic9)));%% -----------------9------------------------sensor_num = 10;R = 4.5879/100;%R = 10/100;d_angle = (0:sensor_num-1)‘*2*pi/sensor_num;theta =1.6*d_angle(2)*180/pi;x = R*cos(d_angle);y = R*sin(d_angle); tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));A = exp(-1j*2*pi*tao/l);n = randn(sensor_num,N_x)*sqrt(10^(-45/10));xt = A*st+n;% -------------------2D-MUSIC算法-----------------------Rx = xt*xt‘/N_x;[U,S] = eig(Rx);disp(est_sour);[~,index] = sort(diag(S));Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);Gn = Un*Un‘;dir=(-180:.25:179.8)*pi/180; ang=(20:.25:91)*pi/180;Pmusic10 = zeros(length(dir),length(ang));for i=1:length(dir) for k=1:length(ang) a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i))); a_theta = exp(-1j*2*pi*a_tao/l); Pmusic10(i,k)=1./abs((a_theta)‘*Gn*a_theta); endendP_music10 = 10*log10(Pmusic10/min(min(Pmusic10)));figure;% subplot(221);% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);% mesh(xx,yy,P_music9);% title(‘9元阵二维空间谱‘);% xlabel(‘俯仰角/°‘);ylabel(‘方向角/°‘);zlabel(‘空间谱/dB‘);% axis([20,91,-180,180,0,24]);%colormap gray;subplot(121);[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);mesh(xx,yy,P_music9);title(‘9元阵方向角空间谱‘);xlabel(‘俯仰角/°‘);ylabel(‘方向角/°‘);zlabel(‘空间谱/dB‘);axis([20,91,-180,180,0,24]);%colormap gray;% subplot(222);% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);% mesh(xx,yy,P_music10);% title(‘10元阵二维空间谱‘);% xlabel(‘俯仰角/°‘);ylabel(‘方向角/°‘);zlabel(‘空间谱/dB‘);% axis([20,91,-180,180,0,24]);%colormap gray;subplot(122);[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);mesh(xx,yy,P_music10);title(‘10元阵方向角空间谱‘);xlabel(‘俯仰角/°‘);ylabel(‘方向角/°‘);zlabel(‘空间谱/dB‘);axis([20,91,-180,180,0,24]);%colormap gray;
8、
% 波程差图clear all;close all;clc;sensor_num = 9;c = 3e8;d_angle = (0:sensor_num-1)‘*2*pi/sensor_num;phi = 90;namda = c/18e9; if mod(sensor_num,2); R = namda/(1-cos(d_angle(2)));%相邻最小间距最大值 %R = namda/(cos(floor(sensor_num/4)*d_angle(2))+cos(d_angle(2)*floor((sensor_num+1)/4)-.5));%相邻最大间距最大值 %R = namda/(1+cos(d_angle(2)/2));%任意阵元间距最大 else R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));%相邻最小间距最大值 %R = namda/(2*sin(d_angle(2)/2));%相邻最大间距最大值 %R = namda/2;%任意阵元间距最大值 end R = 1;dir = (0:.4:d_angle(2)*180/pi)*pi/180;%idx=nchoosek(1:sensor_num,2); % 16取2的组合idx1 = [1,2;2,9;9,3;3,8;8,4;4,7;7,5;5,6];idx2 = [2,1;1,3;3,9;9,4;4,8;8,5;5,7;7,6];r = zeros(length(dir),size(idx1,1));x = R*cos(d_angle);y = R*sin(d_angle);for i = 1:50 tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180)); r(i,:) = -diff(tao(idx1),1,2);endfor i = 51:101 tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180)); r(i,:) = -diff(tao(idx2),1,2);endfigure;subplot(211);plot(dir(1:5:length(dir))*180/pi,0.234*ones(1,length(dir(1:5:length(dir)))),‘.k‘,dir*180/pi,r(:,end),‘-.k‘,dir*180/pi,r(:,1:end-1),‘k‘);grid on;legend(‘0.234‘,‘k=1‘,‘k=其他‘);title(‘(a)9元阵投影间隔‘)xlabel(‘入射方向\°‘);ylabel(‘投影间隔/R‘);axis([0,40,-0.01,0.7]);% ---------------------------------------------------------------------sensor_num = 10;d_angle = (0:sensor_num-1)‘*2*pi/sensor_num;dir = linspace(0,d_angle(2),101);%idx=nchoosek(1:sensor_num,2); % 16取2的组合idx1 = [1,2;2,10;10,3;3,9;9,4;4,8;8,5;5,7;7,6];idx2 = [2,1;1,3;3,10;10,4;4,9;9,5;5,8;8,6;6,7];r = zeros(length(dir),size(idx1,1));x = R*cos(d_angle);y = R*sin(d_angle);for i = 1:50 tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180)); r(i,:) = -diff(tao(idx1),1,2);endfor i = 51:101 tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180)); r(i,:) = -diff(tao(idx2),1,2);endsubplot(212);plot(dir(1:5:length(dir))*180/pi,0.3633*ones(1,length(dir(1:5:length(dir)))),‘.k‘,dir*180/pi,r(:,2),‘-.k‘,dir*180/pi,r(:,[1,3:6,7]),‘k‘);grid on;legend(‘0.3633‘,‘k=1‘,‘k=其他‘);title(‘(b)10元阵投影间隔‘);xlabel(‘入射方向/°‘);ylabel(‘投影间隔/R‘);axis([0,36,-0.01,0.7]);
DOA——MUSIC算法
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。