clear
clc
f0=200*10^3;
B=20*10^3;
T=0.01;
fl=190*10^3;
fh=210*10^3;
fs=4*fh;
c=1500;
% source=1;
N=24;      %陣源數(shù)
Ns=T*fs;   %快拍數(shù)
d=0.003826;% 陣元間距 
sita=10;   %入射角度
K=20;
fm=fl:B/K:fh;
w=2*pi*fm;
ph1=-90;                   
ph2=90;                   
phn=0.5;  %掃描步長                
ax=(ph2-ph1)/phn+1;
%*******************產(chǎn)生信號********************************
for n=1:N;
       for k=1:Ns;
       p(k,n,:)=exp(j*(2*pi*f0*(k-1)/fs+pi*(((k-1)/fs)^2)*B/T+w*(n-1)*d*sin(sita*pi/180)/c));
       end
end 
%********************* 相位補償*************************************
for i=1:ax;     
    num=ph1+(i-1)*phn;
    fain=num*pi/180;
  for l=1:K+1;
      w=2*pi*(fl+(l-1)*1000);
      P=p(1:4,:,l);
      R=P'*P/4;  
%       [em,zm]=eig(R);                            
%       [zm1,pos1]=max(zm);                         
%       for l=1:source;
%          [zm2,pos2]=max(zm1);                    
%           zm1(:,pos2)=[];                          
%          em(:,pos2)=[];                           
%       end
      for n=1:N;
          x(n,i)=exp(-j*w*d*(n-1)*sin(fain)/c); 
      end
      v=x(:,i);
    pmusic(i,l)=v'*R*v;
  end
  pmusic(i)=sum(pmusic(i,:))/(K+1);
end
i=1:ax; 
p1=max(pmusic);
D1=pmusic/p1;
zz1=10*log10(D1);

% for l=1:K+1;
% D1(:,l)=pmusic(:,l)./p1(l);
% zz1(:,l)=10*log10(D1(:,l));
% end
figure;
i=1:ax;
plot(ph1+(i-1)*phn,zz1);
grid on 
title('信噪比為10dB，MUSIC')
xlabel('角度')
ylabel('dB')