clear all;%清除屏幕
%k=0.5;% d/Y,Y設為波長
N=10;%考慮１０個陣元
BWnn=1.732/N;
%迭代次數--波數，橫軸坐標
%K=[10 100 200 300 400 500 600 700 800 900 1000];
kn=100;
m_carlo=100;%Monte-Carlo數
dir_signal=0;
dir_gan=zeros(17,1);%干擾信號的到達角度  *******本題中把干擾信號當成第二個信號*******
for m=1:1:length(dir_gan)
   dir_gan(m)=(m-1)*0.025;
   dir_BW(m)=dir_gan(m)/BWnn;
end
%SNR
snr_db=20;%SNR,噪聲功率為N=1，矩陣形式
uu_b=zeros(length(dir_gan),m_carlo);
uu_c=zeros(length(dir_gan),m_carlo);
uu_m=zeros(length(dir_gan),m_carlo);
uu_r=zeros(length(dir_gan),m_carlo);
v = zeros(N,1);
vz=zeros(N,1);
u = zeros(N,1);%干擾的陣簇矢量
uz=zeros(N,1);
n = zeros(N,1);%噪聲
w = zeros(N,1);%隨角度變化的陣簇矢量
wr= zeros(N,1);
dir=zeros(1001,1);
P_b=zeros(length(dir),1);
P_c=zeros(length(dir),1);
P_m=zeros(length(dir),1);
P_r=zeros(length(dir),1);
%信號方向，從-1到1，間隔為0.01
for a=1:1:1001%u從-0.05到0.05
   dir(a)=(a-1)*(0.1/1000)-0.05;
end
V=zeros(N,N);%特征向量矩陣
D=zeros(N,N);%特征值矩陣
Vr=zeros(N,N);%特征向量矩陣
Dr=zeros(N,N);

zv=exp(j*pi*dir_signal);%
si_amp = sqrt(10^(snr_db/10));
in_amp = sqrt(10^(snr_db/10));%干擾的幅值

for i_gan=1:1:length(dir_gan)
   zu=exp(j*pi*dir_gan(i_gan));
    for i = 1:N
    v(i) = exp(j*pi*(-(N+1)/2+i)*dir_signal);
    u(i) = exp(j*pi*(-(N+1)/2+i)*dir_gan(i_gan));
    vz(i) = zv^(i-1);
    uz(i) = zu^(i-1);
    end 
         %計算協方差矩陣
      for j_monte=1:1:m_carlo
         signal = si_amp*exp(j*2*pi*randn);
         interupt = in_amp*exp(j*2*pi*randn);%相位隨機產生一個干擾
         for i_temp = 1:1:kn 
         signal_amp = signal*((randn+j*randn)/(sqrt(2)));
         interupt_amp = interupt*((randn+j*randn)/(sqrt(2)));  %幅度隨機產生一個干擾
         noise=(randn(N,1)+j*randn(N,1))/(sqrt(2)); 
         x=signal_amp*v+interupt_amp*u+noise;%信號+干擾+噪聲
         xr=signal_amp*vz+interupt_amp*uz+noise;
         if i_temp==1
            Cx=x*(x');  
            Cxr=xr*xr';
         else
            Cx=((i_temp-1)*Cx+x*x')/i_temp;
            Cxr=((i_temp-1)*Cxr+xr*xr')/i_temp;
         end 
         end
   %對協方差矩陣進行特征值、特征向量分解
       [V,D]=eig(Cx);
       Us=[V(:,1),V(:,2)];%最大特征值對應的特征向量
       [Vr,Dr]=eig(Cxr);
       Usr=[Vr(:,1),Vr(:,2)];
       %求隨角度變化的陣簇矢量
       for j_dir = 1:1:length(dir)
          z=exp(j*pi*dir(j_dir));
          for i = 1:N
             w(i) = exp(j*pi*(-(N+1)/2+i)*dir(j_dir));
             wr(i)=z^(i-1);
          end
          P_b(j_dir)=abs(w'*Cx*w);          
          P_c(j_dir)=abs(1/(w'*inv(Cx)*w));        
          P_m(j_dir)=abs(1/(w'*(eye(10,10)-Us*Us')*w));
          P_r(j_dir)=abs(1/(wr'*(eye(10,10)-Usr*Usr')*wr));
       end
       for i_dir= 1:1:length(dir)
          if P_b(i_dir)==max(P_b);
             uu_b(i_gan,j_monte)=dir(i_dir);
          end

          if P_c(i_dir)==max(P_c);
             uu_c(i_gan,j_monte)=dir(i_dir);
          end

          if P_m(i_dir)==max(P_m);
             uu_m(i_gan,j_monte)=dir(i_dir);
          end
          
          if P_r(i_dir)==max(P_r);
             uu_r(i_gan,j_monte)=dir(i_dir);
          end

       end
 end
end
%計算方差
%uu_m
for i=1:1:length(dir_gan)
   sum_b=0; 
   ave_b=0;
   vartotle_b=0;

   sum_c=0; 
   ave_c=0;
   vartotle_c=0;

   sum_m=0; 
   ave_m=0;
   vartotle_m=0;
   
   sum_r=0; 
   ave_r=0;
   vartotle_r=0;

   for j=1:1:m_carlo
      sum_b=uu_b(i,j)+sum_b;
      sum_c=uu_c(i,j)+sum_c;
      sum_m=uu_m(i,j)+sum_m;
      sum_r=uu_r(i,j)+sum_r;
   end
   ave_b=sum_b/m_carlo;
   ave_c=sum_c/m_carlo;
   ave_m=sum_m/m_carlo;
   ave_r=sum_r/m_carlo;
   for t=1:1:m_carlo
      vartotle_b=(uu_b(i,t)-0)^2+vartotle_b;
      vartotle_c=(uu_c(i,t)-0)^2+vartotle_c;
      vartotle_m=(uu_m(i,t)-0)^2+vartotle_m;
      vartotle_r=(uu_r(i,t)-0)^2+vartotle_r;

   end
   var_b(i)=10*log10(vartotle_b/m_carlo);
   var_c(i)=10*log10(vartotle_c/m_carlo);
   var_m(i)=10*log10(vartotle_m/m_carlo);
   var_r(i)=10*log10(vartotle_r/m_carlo);
end
figure
subplot(2,2,1)
plot(dir_BW,var_b,'-*r',dir_BW,var_c,'-.g',dir_BW,var_m,'-b',dir_BW,var_r,':ok');
title('Performance of versus Δu/HPBW(SNR=20db,K=100)','FontWeight','bold');
xlabel('Δu/HPBW','FontWeight','bold');
ylabel('10lg(Var)(db)','FontWeight','bold');
legend('Bartlett','Capon','MUSIC','Root MUSIC',2); 

grid