%Capon spectrum.
clear all;
%參數賦值
c=3e+8;                                               %光速
lmda=0.24;                                         %L波段
fc=c/lmda;                                          %載頻
N=16;                                                 %陣元數
DOA_sigA=10;                                 %信號源1來波方向與陣法線方向夾角 
DOA_sigB=30;                                 %信號源2來波方向與陣法線方向夾角 
d=0.12;                                               %陣元間距
fs=1e+7;                                     %采樣率
chirp_rate=2e+11;                                  %調頻率
%計算s(t)*exp(j*2*pi*fc*t)
t=[0:1/fs:0.0001]';                                     %采樣時間
s1=exp(j*pi*chirp_rate*t.^2);                               %信源1復包絡  
s2=exp(j*pi*chirp_rate*t.^2);                               %信源2復包絡  
% carw=(j*2*pi*fc*t);                                 %載波
% sc1=diag(s)*carw;  
% sc2=diag(s)*carw;  
%產生方向矢量1(steering vector)
sv_sigA=[];
for i=1:N
    sv_sigA=[sv_sigA;exp(j*2*pi*d*(i-1)*sin(DOA_sigA*pi/180)/lmda)];
end
%產生方向矢量2(steering vector)
sv_sigB=[];
for i=1:N
    sv_sigB=[sv_sigB;exp(j*2*pi*d*(i-1)*sin(DOA_sigB*pi/180)/lmda)];
end
%N元等距線陣接收信號矢量
x=sv_sigA*transpose(s1)+sv_sigB*transpose(s2)+randn(N,length(t));
%陣列輸出協方差矩陣
Rx=zeros(N,N);
for i=1:length(t)
    Rx=Rx+x(:,i)*x(:,i)';
end
Rx=Rx/length(t);
%產生capon譜圖
CaponSpec=[];
for sita=-pi/2:1/100:pi/2; 
    sv_sita=[];
    for i=1:N
        sv_sita=[sv_sita;exp(j*2*pi*d*(i-1)*sin(sita)/lmda)];
    end
    CaponSpec=[CaponSpec;1/(sv_sita'*inv(Rx)*sv_sita)];
end  
%繪圖
sita=-pi/2:1/100:pi/2;
sita_deg=sita*180/pi;
CaponSpec_log=10*log10(abs(CaponSpec/max(CaponSpec)));
plot(sita_deg,CaponSpec_log)
grid
axis([-90,90,-10,0])
title('波束形成測向 角譜圖')
xlabel('方位角（度）')
ylabel('輸出功率（dB）')
hold on
stem(DOA_sigA,-10,'*')
stem(DOA_sigB,-10,':')