clc
clear;
f0=8000;   %信號頻率
fs=4*f0;   %采樣頻率 
Ns=64;   %采樣點數 
N1=8;      %陣元數
c=1500;    %聲速 
r1=0.226    %大子陣半徑
r2=0.113     %小子陣半徑
fai=0;     %信號入射方向 
w=2*pi*f0;
ph1=-180;                   
ph2=180;                   
phn=2;                  
ax=(ph2-ph1)/phn+1;

%************產生入射信號****************%
for n=1:N1;
       a=2;
    if mod(n,a)~=0
       for k=1:Ns;
     
        p(n,k)=sin(w*(k-1)/fs-w*r1*cos(fai*pi/180-(n-1)*2*pi/N1)/c);
       end
       
       p1(n,:)=hilbert(p(n,:));
    else
       for k=1:Ns;
        p(n,k)=sin(w*(k-1)/fs-w*r2*cos(fai*pi/180-(n-1)*2*pi/N1)/c);
       end
        p1(n,:)=hilbert(p(n,:));
    end
end
p0=1/8*(p1(1,:)+p1(2,:)+p1(3,:)+p1(4,:)+p1(5,:)+p1(6,:)+p1(7,:)+p1(8,:));
vx0=1/2*((p1(1,:)-p1(5,:))+(p1(2,:)-p1(6,:))*exp(-j*pi/4))*exp(-j*pi/2);
vy0=1/2*((p1(3,:)-p1(7,:))+(p1(4,:)-p1(8,:))*exp(-j*pi/4))*exp(-j*pi/2);
  for k=1:Ns;
 s(k)=sin(w*(k-1)/fs);
  end
pp=hilbert(s);
vxx=pp*cos(fai*pi/180);
vyy=pp*sin(fai*pi/180);
%**************對接收信號進行相移*******%
     Xp1=zeros(ax,Ns);
     Xvx1=zeros(ax,Ns);
     Xvy1=zeros(ax,Ns);
         Xvy1=zeros(ax,Ns);
         
for i=1:ax;     
    num=ph1+(i-1)*phn;
    fain=num*pi/180;  
    for n=1:N1;
        a=2;
      if mod(n,a)~=0 
        Xp01(n,:,i)=p1(n,:)*exp(j*w*r1*cos(fain-(n-1)*2*pi/N1)/c); 
        Xp1(i,:)=Xp1(i,:)+Xp01(n,:,i);
        Xvx01(n,:,i)=Xp01(n,:,i)*cos(fai*pi/180);
        Xvy01(n,:,i)=Xp01(n,:,i)*sin(fai*pi/180);
         Xvx1(i,:)=Xvx1(i,:)+Xvx01(n,:,i);
          Xvy1(i,:)=Xvy1(i,:)+Xvy01(n,:,i);
      else
        Xp01(n,:,i)=p1(n,:)*exp(j*w*r2*cos(fain-(n-1)*2*pi/N1)/c); 
        Xp1(i,:)=Xp1(i,:)+Xp01(n,:,i);
       Xvx01(n,:,i)=Xp01(n,:,i)*cos(fai*pi/180);
        Xvy01(n,:,i)=Xp01(n,:,i)*sin(fai*pi/180);
        Xvx1(i,:)=Xvx1(i,:)+Xvx01(n,:,i);
        Xvy1(i,:)=Xvy1(i,:)+Xvy01(n,:,i);
      end 
    end 
     vcc0(i,:)=(Xvx1(i,:)+vxx)*cos(fain)+(Xvy1(i,:)+vyy)*sin(fain);
   Xpp(i)=(sum(abs((Xp1(i,:)+pp)').^2));
    Xvcc0(i)=(sum(abs(vcc0(i,:)').^2));
    Rpp00(i)=(Xpp(i)+Xvcc0(i))*Xvcc0(i);
    vc0(i,:)=(Xvx1(i,:)+vx0)*cos(fain)+(Xvy1(i,:)+vy0)*sin(fain);
    Xvc1(i,:)=Xvx1(i,:)*cos(fain)+Xvy1(i,:)*sin(fain);
    Xp0(i)=(sum(abs((Xp1(i,:)+p0)').^2));
    Xvc0(i)=(sum(abs(vc0(i,:)').^2));
    Rp0(i)=(Xp0(i)+Xvc0(i))*Xvc0(i);
 Xp=(sum(abs(Xp1').^2));
Xvc=(sum(abs(Xvc1').^2));
Rp(i)=(Xp(i)+Xvc(i))*Xvc(i);
end 
i=1:ax;
  num=ph1+(i-1)*phn;
  fain=num*pi/180;
    maxRpp0=max(Rpp00);
Rppp0=10*log10(Rpp00(i)./maxRpp0);
  
  maxRp0=max(Rp0);
Rpp0=10*log10(Rp0(i)./maxRp0);
maxRp=max(Rp);
Rpp=10*log10(Rp(i)./maxRp);
figure
plot(ph1+(i-1)*phn,Rpp,ph1+(i-1)*phn,Rpp0,ph1+(i-1)*phn,Rppp0);
legend('交叉矢量陣(p+vc)*vc','jiaocha',4);
axis([-180 180 -100 0])
grid on
hold on
xlabel('角度')
ylabel('dB')