shuangshizizhen.m

双十字阵矢量水听器阵的指向性研究

%%%双十字阵2相干信源music和常规普
clc 
close all;
clear all;
c=1500;
f0=8000;
fs=100000;
snap=1024;
lamda=c/f0;
d=0.23;%阵元间距
M=8; %阵元数
p=2;
%构造信源
n=1:snap;
x1=[exp(j*2*pi*f0*n/fs);exp(j*2*pi*f0*n/fs)];%构造2个复信号
% x1=[exp(j*2*pi*f0*n/fs+randn(1,snap));exp(j*2*pi*f0*n/fs+randn(1,snap))];%构造2个复信号
%产生信号与噪声
theta=[90,130];%方位角
SN=[10,10];    %信噪比
rad=pi/180;%弧度
nr=randn(M,snap);
ni=randn(M,snap);
noise=nr+j*ni;%产生复噪声

nr1=randn(3*M,snap);
ni1=randn(3*M,snap);
noise1=nr1+j*ni1;%产生复噪声

Am=sqrt(2)*diag(10.^(SN/20));
S=Am*x1;
h=[[1,1];cos(theta*rad);sin(theta*rad)];
ca=[cos(theta*rad); cos(theta*rad+pi/2); cos(theta*rad+pi); cos(theta*rad+pi*3/2)];
cb=0.5.*[cos(theta*rad+pi/4); cos(theta*rad+pi*3/4); cos(theta*rad+pi*5/4); cos(theta*rad+pi*7/4)];
cc=[ca;cb];
A=exp(j*d*2*pi*cc/lamda);  %俯仰角相位差
A0=kron(A(:,1),h(:,1));%克罗内可积
A1=kron(A(:,2),h(:,2));%克罗内可积
Aa=[A0,A1];
X0=(Aa*S)+noise1;%第k个矢量水听器上的接收信号

Rx0=X0*X0'/snap;%接收信号的协方差
[E0,C0,V0]=svd(Rx0);%计算特征值
Vs0=E0(:,1:p);          
Vn0=E0(:,p+1:3*M);%小特征值对应的特征向量，即噪声子空间

X1=(A*S)+noise;%第k个声压水听器上的接收信号
Rx1=X1*X1'/snap;%接收信号的协方差
%Rx1=top_p(M,Rx1);
[E1,C1,V1]=svd(Rx1);%计算特征值
Vs1=E1(:,1:p);          
Vn1=E1(:,p+1:M);%小特征值对应的特征向量，即噪声子空间
for the=0:1:360
   th=fix(the*1+1);
        ca1=[cos(the*rad); cos(the*rad+pi/2); cos(the*rad+pi); cos(the*rad+pi*3/2)];
        cb1=0.5.*[cos(the*rad+pi/4); cos(the*rad+pi*3/4); cos(the*rad+pi*5/4); cos(the*rad+pi*7/4)];
        cc1=[ca1;cb1];
        ap=exp(j*d*2*pi*cc1/lamda);                     
        h1=[1;cos(the*rad);sin(the*rad)];
        aa=kron(ap,h1);       
        p0(th)=1/(aa'*Vn0*Vn0'*aa);
        p1(th)=1/(ap'*Vn1*Vn1'*ap);
        p2(th)=aa'*Rx0*aa;
        p3(th)=ap'*Rx1*ap;
        p4(th)=sum(aa)*sum(aa)';
 end
pp0=10*log10(abs(p0)/max(abs(p0)));
pp1=10*log10(abs(p1)/max(abs(p1)));
pp2=10*log10(abs(p2)/max(abs(p2)));
pp3=10*log10(abs(p3)/max(abs(p3)));
figure(1);
plot(0:360,pp0,0:360,pp1,0:360,pp2,0:360,pp3);
legend('矢量面阵MUSIC谱','声压面阵MUSIC谱','矢量面阵常规谱','声压面阵常规谱');
xlabel('方位角');
ylabel('谱/dB');
title('波束形成');
pp4=10*log10((p4)/max((p4)));
figure(2);
plot(0:360,pp4);
xlabel('方位角');
ylabel('谱/dB');
axis([0,360,-50,0]);