%%%%%%信源數目估計----蓋而圓半徑法page46－48
clc;
clear ;
close all;
tic
M=16;%陣列的天線數
N=4;%信源數
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;
theta3=15;
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*u0*t.^2));
s2=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s3=exp(j*2*pi*f0*t);


%陣列流行矢量

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_theta3=exp(j*2*pi*d/lamda*a*sin(theta3/180*pi));

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

snr=10;%信噪比

%產生滿足一定信躁比的噪聲
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/20)*noise;

X=X0+noise;
R_x=X*X'/length(t);
R_pie=R_x(1:M-1,1:M-1);
[V,D]=eig(R_pie);
T=[V zeros(M-1,1);zeros(1,M-1) 1];%%%%%%%%構造矩陣T
Rt=T'*R_x*T;
Rou=abs(Rt(1:M-1,M));
temp0=sum(Rou);

for n=0:M-1    %%%%%%%%%M次迭代再求平均
    for k=1:M-1
        dl=0.5;
        GDE(n+1,k)=Rou(k)-dl/(M-1)*temp0;%%%%%%%%%式（2.5.25）
    end
    K0=find(GDE(n+1,:)>=0);%%%%%%%%%找出第一個大于零的數所對應的位置
    K(1,n+1)=M-K0(1,1);       %%%%%%%%%信源數目的估計
end
figure(1);plot(1:M-1,GDE(1:M,:));grid on;
figure(2);plot(0:M-1,K,'.');grid on;
disp('M次迭代后信源數目估計的平均值');
mean(K) %%%%%%%%%%%%平均的結果顯示出來