【未完待续】
一、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与信源3 Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3 Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3 save 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);%无噪声信号功率--%矩阵维数=信源数*1 delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1 % delta1=ps./(2*10.^(sn/10)); %----单信号源 delta2=diag(delta1);%矩阵维数=1*1 delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1 Rev_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,:)); end Rev_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]‘;%矩阵维数=阵元数*1 tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1 a1=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); end sub_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与信源3 Rcos2=Rcos21;%构造相干信源--信源1、信源2与信源3 Rcos3=Rcos31;%构造相干信源--信源1、信源2与信源3 save 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);%无噪声信号功率--%矩阵维数=信源数*1 delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1 % delta1=ps./(2*10.^(sn/10)); %----单信号源 delta2=diag(delta1);%矩阵维数=1*1 delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1 Rev_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,:)); end Rev_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]‘;%矩阵维数=阵元数*1 tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1 a1=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与信源3 Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3 Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3 save 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);%无噪声信号功率--%矩阵维数=信源数*1 delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1 % delta1=ps./(2*10.^(sn/10)); %----单信号源 delta2=diag(delta1);%矩阵维数=1*1 delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1 Rev_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,:)); end Rev_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]‘;%矩阵维数=阵元数*1 tmp4=[0:L-1]‘;%子矩阵维数=子矩阵阵元数*1 a1=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; end end Rxx_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 end figure; % 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 end end figure; 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 end figure; 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); end end P_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); end end P_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); end for 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); end figure; 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); end for 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); end subplot(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]);
原文:http://www.cnblogs.com/xingshansi/p/7163605.html