hushifenbianjiaodu.m

矢量水听器波束形成程序 分辨角度

%%%%%%%%%%------------准备湖试---%%%%%%%%%%%%%%
%%---------x和y轴测质点阵速的传感器,稀疏（均匀）布放阵元在z轴上----------%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        1个声压传感器加上2个阵速传感器 2维方位角分辨                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;close all;clc;
M=3;
c=1500;
f=6.3*10^3;
T=1;
fs=25.6*10^3*2.56;
theta=[90*pi/180,90*pi/180];  %%%%%%%%%%3D越接近90度越好
fi=[60*pi/180,200*pi/180];    %%%%%%%%%%3D越接近180度越好
%sita=0:pi/180:2*pi;
L=T*fs+1;
lamd=c/f;%%%%%%%%%%%%%%%%%%%%%%   wavelength
%d=lamd/2;%%%%%%%%%%%%%%%%%%%%%%%%  sensor spacing
d=0.3;
n=1:M;
Ln=1:L;
snr=10;
r=d*(n-1);
%r=[-d,0,d];
rou=10^(snr/10);
ap1=exp(j*2*pi*f*r*cos(theta(1))/c);%%%%---inverse （0，pi）-------%%%%%
ap2=exp(j*2*pi*f*r*cos(theta(2))/c);
Ap1=ap1.';
Ap2=ap2.';
u1=[cos(fi(1))*sin(theta(1)) sin(fi(1))*sin(theta(1))].';
u2=[cos(fi(2))*sin(theta(2)) sin(fi(2))*sin(theta(2))].';
%u=[cos(fi)*sin(theta) sin(fi)*sin(theta) cos(theta)].';%%%%%-- 3 compnent of  velocity sensor inverse
h1=[1 u1.'].';
h2=[1 u2.'].';
%h=1;%%%%%%%%%%%%%%%%%scalar sensor arrays
Av1=kron(Ap1,h1);
Av2=kron(Ap2,h2);
Am=sqrt(2*10^(snr/10));
ss1=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
ss2=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
ss=[];
ss=[ss;ss1;ss2];
Av=[];
Av=[Av,Av1,Av2];
y=Av*ss+(randn(3*M,L)+j*randn(3*M,L));%--2 compnent of  velocity vector sensor plus 1 scalar sensor
R1=y*y'/L;
R=R1;
%%%%%%%%%%%%%%%%%%%%%%5-------------music---------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U,V]=eig(R);
[VV,index]=sort(diag(V));
U=fliplr(U(:,index));
Un=U(:,3:(3*M));
%%%%%%%%%%%%%%%%%%%------------------eig-----------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Pm=[];
for fai=0:.1:360
    fai1=fix(fai*10+1);
%  for thita=0:.1:180
%      thita1=fix(thita*10+1);
     ap1=exp(j*2*pi*f*r*cos(theta(1))/c);
    Ap1=ap1.';
      u1=[cos(fai*pi/180)*sin(theta(1)) sin(fai*pi/180)*sin(theta(1))].';%%%--inverse （0，pi）
      %u1=[cos(fi(1))*sin(thita*pi/180) sin(fi(2))*sin(thita*pi/180)].';
     h1=[1 u1.'].';
%     %h1=1;%%%%%%%%%%%scalar sensor arrays
    Av1=kron(Ap1,h1);
     %Pm(fai1)=1/(Av1'*Un*Un'*Av1);
     Pb(fai1)=Av1'*R*Av1;
     Pm(fai1)=1/(Av1'*Un*Un'*Av1);
     %Pm(fai1,thita1)=1/(Av1'*Un*Un'*Av1);
 end
%end
%Pvb=10*log10(abs(Pbv)/max(max(abs(Pbv))));
 Pvb=10*log10(abs(Pb)/max(max(abs(Pb))));
 Pvm=10*log10(abs(Pm)/max(max(abs(Pm))));
 boshuDB=-3*ones(1,3601);
 figure
plot(0:.1:360,Pvb);hold on;
% plot(0:.1:180,boshuDB);hold off;
 %figure
plot(0:.1:360,Pvm,'r');hold on;
plot(0:.1:360,boshuDB,'b-.');%hold off;
legend('cbf','music','-3dB');
%  title('z轴放阵，3元矢量阵MUSIC谱，方位角[100,116]度，俯仰角[90,90]度，snr=10dB')
% ylabel('Elevation(deg)');
% xlabel('Azimuth(deg)');
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %       1个声压传感器加上2个阵速传感器 2维俯仰角分辨 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % clear all;close all;clc;
% % M=3;
% % c=1500;
% % f=6.3*10^3;
% % T=1;
% % fs=25.6*10^3*2.56;
% % theta=[90*pi/180,87*pi/180];  %%%%%%%%%%3D越接近90度越好
% % fi=[60*pi/180,60*pi/180];    %%%%%%%%%%3D越接近180度越好
% % %sita=0:pi/180:2*pi;
% % L=T*fs+1;
% % lamd=c/f;%%%%%%%%%%%%%%%%%%%%%%   wavelength
% % %d=lamd/2;%%%%%%%%%%%%%%%%%%%%%%%%  sensor spacing
% % d=0.3;
% % n=1:M;
% % Ln=1:L;
% % snr=10;
% % r=d*(n-1);
% % %r=[-d,0,d];
% % rou=10^(snr/10);
% % ap1=exp(j*2*pi*f*r*cos(theta(1))/c);%%%%---inverse （0，pi）-------%%%%%
% % ap2=exp(j*2*pi*f*r*cos(theta(2))/c);
% % Ap1=ap1.';
% % Ap2=ap2.';
% % u1=[cos(fi(1))*sin(theta(1)) sin(fi(1))*sin(theta(1))].';
% % u2=[cos(fi(2))*sin(theta(2)) sin(fi(2))*sin(theta(2))].';
% % %u=[cos(fi)*sin(theta) sin(fi)*sin(theta) cos(theta)].';%%%%%-- 3 compnent of  velocity sensor inverse
% % h1=[1 u1.'].';
% % h2=[1 u2.'].';
% % %h=1;%%%%%%%%%%%%%%%%%scalar sensor arrays
% % Av1=kron(Ap1,h1);
% % Av2=kron(Ap2,h2);
% % Am=sqrt(2*10^(snr/10));
% % ss1=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
% % ss2=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
% % ss=[];
% % ss=[ss;ss1;ss2];
% % Av=[];
% % Av=[Av,Av1,Av2];
% % y=Av*ss+(randn(3*M,L)+j*randn(3*M,L));%--2 compnent of  velocity vector sensor plus 1 scalar sensor
% % R1=y*y'/L;
% % R=R1;
% % %%%%%%%%%%%%%%%%%%%%%%5-------------music---------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % [U,V]=eig(R);
% % [VV,index]=sort(diag(V));
% % U=fliplr(U(:,index));
% % Un=U(:,3:(3*M));
% % %%%%%%%%%%%%%%%%%%%------------------eig-----------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %Pm=[];
% % % for fai=0:.1:360
% % %     fai1=fix(fai*10+1);
% %  for thita=0:.1:180
% %      thita1=fix(thita*10+1);
% %      ap1=exp(j*2*pi*f*r*cos(thita*pi/180)/c);
% %     Ap1=ap1.';
% %       %u1=[cos(fai*pi/180)*sin(theta(1)) sin(fai*pi/180)*sin(theta(1))].';%%%--inverse （0，pi）
% %      u1=[cos(fi(1))*sin(thita*pi/180) sin(fi(2))*sin(thita*pi/180)].';
% %      h1=[1 u1.'].';
% % %     %h1=1;%%%%%%%%%%%scalar sensor arrays
% %     Av1=kron(Ap1,h1);
% %      %Pm(fai1)=1/(Av1'*Un*Un'*Av1);
% %      Pb(thita1)=Av1'*R*Av1;
% %      Pm(thita1)=1/(Av1'*Un*Un'*Av1);
% %      %Pm(fai1,thita1)=1/(Av1'*Un*Un'*Av1);
% %  end
% % %end
% % %Pvb=10*log10(abs(Pbv)/max(max(abs(Pbv))));
% %  Pvb=10*log10(abs(Pb)/max(max(abs(Pb))));
% %  Pvm=10*log10(abs(Pm)/max(max(abs(Pm))));
% %  boshuDB=-3*ones(1,3601);
% %  figure
% %  %plot(0:.1:360,Pvb);hold on;
% % % plot(0:.1:180,boshuDB);hold off;
% %  %figure
% % plot(0:.1:180,Pvm,'r');hold on;
% % plot(0:.1:180,boshuDB,'b-.');%hold off;
% % %legend('cbf','music','-3dB');
% % %  title('z轴放阵，3元矢量阵MUSIC谱，方位角[100,116]度，俯仰角[90,90]度，snr=10dB')
% % % ylabel('Elevation(deg)');
% % % xlabel('Azimuth(deg)');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %       1个声压传感器加上2个阵速传感器 2维俯仰角分辨 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear all;close all;clc;
% M=3;
% c=1500;
% f=6.3*10^3;
% T=1;
% fs=25.6*10^3*2.56;
% theta=[80*pi/180,80*pi/180];  %%%%%%%%%%3D越接近90度越好
% fi=[60*pi/180,73*pi/180];    %%%%%%%%%%3D越接近180度越好
% %sita=0:pi/180:2*pi;
% L=T*fs+1;
% lamd=c/f;%%%%%%%%%%%%%%%%%%%%%%   wavelength
% %d=lamd/2;%%%%%%%%%%%%%%%%%%%%%%%%  sensor spacing
% d=0.3;
% n=1:M;
% Ln=1:L;
% snr=10;
% r=d*(n-1);
% %r=[-d,0,d];
% rou=10^(snr/10);
% ap1=exp(j*2*pi*f*r*cos(theta(1))/c);%%%%---inverse （0，pi）-------%%%%%
% ap2=exp(j*2*pi*f*r*cos(theta(2))/c);
% Ap1=ap1.';
% Ap2=ap2.';
% u1=[cos(fi(1))*sin(theta(1)) sin(fi(1))*sin(theta(1))].';
% u2=[cos(fi(2))*sin(theta(2)) sin(fi(2))*sin(theta(2))].';
% %u=[cos(fi)*sin(theta) sin(fi)*sin(theta) cos(theta)].';%%%%%-- 3 compnent of  velocity sensor inverse
% h1=[1 u1.'].';
% h2=[1 u2.'].';
% %h=1;%%%%%%%%%%%%%%%%%scalar sensor arrays
% Av1=kron(Ap1,h1);
% Av2=kron(Ap2,h2);
% Am=sqrt(2*10^(snr/10));
% ss1=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
% ss2=Am*exp(j*(2*pi*f/fs*Ln+2*pi*rand(1,L)));
% ss=[];
% ss=[ss;ss1;ss2];
% Av=[];
% Av=[Av,Av1,Av2];
% y=Av*ss+(randn(3*M,L)+j*randn(3*M,L));%--2 compnent of  velocity vector sensor plus 1 scalar sensor
% R1=y*y'/L;
% R=R1;
% %%%%%%%%%%%%%%%%%%%%%%5-------------music---------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [U,V]=eig(R);
% [VV,index]=sort(diag(V));
% U=fliplr(U(:,index));
% Un=U(:,3:(3*M));
% %%%%%%%%%%%%%%%%%%%------------------eig-----------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Pm=[];
% for fai=0:1:360
%     fai1=fix(fai*1+1);
%  for thita=0:1:180
%      thita1=fix(thita*1+1);
%      ap1=exp(j*2*pi*f*r*cos(thita*pi/180)/c);
%     Ap1=ap1.';
%       %u1=[cos(fai*pi/180)*sin(theta(1)) sin(fai*pi/180)*sin(theta(1))].';%%%--inverse （0，pi）
%      u1=[cos(fai*pi/180)*sin(thita*pi/180) sin(fai*pi/180)*sin(thita*pi/180)].';
%      h1=[1 u1.'].';
% %     %h1=1;%%%%%%%%%%%scalar sensor arrays
%     Av1=kron(Ap1,h1);
%      %Pm(fai1)=1/(Av1'*Un*Un'*Av1);
% %      Pb(thita1)=Av1'*R*Av1;
% %      Pm(thita1)=1/(Av1'*Un*Un'*Av1);
%      Pm(fai1,thita1)=1/(Av1'*Un*Un'*Av1);%music
%      %Pm(fai1,thita1)=(Av1'*R*Av1);%cbf
%  end
% end
% %Pvb=10*log10(abs(Pbv)/max(max(abs(Pbv))));
% % Pvb=10*log10(abs(Pb)/max(max(abs(Pb))));
%  Pvm=10*log10(abs(Pm)/max(max(abs(Pm))));
%  %boshuDB=-3*ones(1,3601);
%  figure
%  [x,y]=meshgrid(0:180,0:360);
% mesh(y,x,Pvm);