%%%雙十字陣2相干信源music和常規普
clc 
close all;
clear all;
c=1500;
f0=8000;
fs=65536;
snap=1024;
lamda=c/f0;
d=0.226;%陣元間距
M=8; %陣元數
p=2;
%構造信源
n=1:snap;
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]);