clc;
clear ;
close all;
tic
M=16;%陣列的天線數
N=3;%信源數
snap=1000;%快拍數目
L=snap;
C=3e8;
lamda=0.2;
f0=C/lamda;
d=0.5*lamda;
% k=d/lamda;
theta0=5;
theta1=20;
theta2=40;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]';%陣列矢量
u0=5;
u1=10;
u2=20;

%相干信號源 
s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s1=exp(j*2*pi*(f0*t+0.5*u1*t.^2));
s2=exp(j*2*pi*(f0*t+0.5*u2*t.^2));

%陣列流行矢量

a_theta0=exp(j*2*pi*d/lamda*a*sin(theta0/180*pi));
a_theta1=exp(j*2*pi*d/lamda*a*sin(theta1/180*pi));
a_theta2=exp(j*2*pi*d/lamda*a*sin(theta2/180*pi));

A=[a_theta0 a_theta1 a_theta2];%子陣的導向矢量，每個子陣完全相同
S=[s0;s1;s2];
X0=A*S;
R_s=S*S';

snr=[-20 -10 0 10 20]; %信噪比
for jj=1:length(snr)
    %產生滿足一定信躁比的噪聲
    randn('state',0);
    real_noise=randn(size(X0));
    randn('state',3);
    imag_noise=randn(size(X0));
    noise=(real_noise+j*imag_noise)/2^0.5;
    noise=10^(-snr(jj)/20)*noise;

    %%%%%%%%%####傳統的MUSIC算法
    X=X0+noise;
    R_x=X*X'/length(t);
    [V,D]=eig(R_x);
    d=diag(D);

    for n=1:M-1
        temp0=sum(d(n+1:M))/(M-n);
        temp1=1;
        for ii=n+1:M
            temp1=temp1*d(ii);
        end
        temp1=temp1^(1/(M-n));
        AIta=temp0/temp1;
        %%%%%%%%%%%%%%%%%%
        AIC(jj,n)=2*L*(M-n)*log(AIta)+2*n*(2*M-n);
        MDL(jj,n)=L*(M-n)*log(AIta)+1/2*n*(2*M-n)*log(L);
        HQ(jj,n)=L*(M-n)*log(AIta)+1/2*n*(2*M-n)*log(log(L));
        EDC(jj,n)=L*(M-n)*log(AIta)+n*(2*M-n)*1/2*log(log(L));
    end
end    
figure(1);
plot(1:M-1,AIC(1:length(snr),:));grid on;title('AIC');
legend('-20dB','-10dB','0dB','10dB','20dB');
figure(2);
plot(1:M-1,MDL(1:length(snr),:));grid on;title('MDL');
legend('-20dB','-10dB','0dB','10dB','20dB');
figure(3);
plot(1:M-1,HQ(1:length(snr),:));grid on;title('HQ');
legend('-20dB','-10dB','0dB','10dB','20dB');
figure(4);
plot(1:M-1,EDC(1:length(snr),:));grid on;title('EDC');
legend('-20dB','-10dB','0dB','10dB','20dB');