close all
clear all
clc
%模擬信號
N=14;     %陣元數
K=200;    %快拍數
M=3;      %信源數
theta=[10,30,80]; %第一個是期望信號，后面的都是干擾
Xt=zeros(N,K);
Snr=0;     %信噪比
Sir1=-10;  %信干比1
Sir2=-20;  %信干比2
Pd=1;      %給出信號功率
Pn=Pd/(10^(Snr/10));    %計算噪聲功率
Pi1=Pd/(10^(Sir1/10));  %計算干擾1功率
Pi2=Pd/(10^(Sir2/10));  %計算干擾2功率
theta=theta*pi/180;
ai1=exp(-j*[0:N-1]'*sin(theta(2))*pi);       %陣列導向矢量（干擾信號1）
ai2=exp(-j*[0:N-1]'*sin(theta(3))*pi);       %陣列導向矢量（干擾信號2）
ad=exp(-j*[0:N-1]'*sin(theta(1))*pi);        %陣列導向矢量（期望信號）
Sd=sqrt(Pd/2)*(randn(1,K)+j*randn(1,K));     %構造期望信號包絡
Si1=sqrt(Pi1/2)*(randn(1,K)+j*randn(1,K));   %構造干擾信號1包絡
Si2=sqrt(Pi2/2)*(randn(1,K)+j*randn(1,K));   %構造干擾信號2包絡
Xt=ad*Sd+ai1*Si1+ai2*Si2;                    %構造輸入信號（期望信號+干擾）
Xt=Xt+sqrt(Pn/2)*(randn(N,K)+j*randn(N,K));  %構造輸入信號（期望信號+干擾+噪聲）
%以上是模擬信號產生
%Music DOA
Rx=Xt*Xt'/K;
[a,b,c]=svd(Rx);
Un=a(:,M+1:N);
Pn=Un*Un';
thetas=-90:90;
for m=1:length(thetas)
    tm=thetas(m)*pi/180;
    am=exp(-j*pi*sin(tm)*[0:N-1]');
    S_Music(m)=1/norm(Pn*am);
end
%Capon Beamformming
Rx=Xt*Xt'/K;
Wopt=inv(Rx)*ad;
thetas=-90:90;
for m=1:length(thetas)
    tm=thetas(m)*pi/180;
    am=exp(-j*pi*[0:N-1]'*sin(tm));
    S_Capon(m)=10*log(abs(Wopt'*am));
end
plot(thetas,S_Capon,'r--','linewidth',2)
grid
hold on
plot(thetas,S_Music,'b-','linewidth',2)
legend('Capon','Music',2)