%%%%%%%%%%%%%%%%%%%%%%%%%  估計信號源數目    %%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%% 
clear;clc
l=1.8;%波長
d=0.5*l;%陣元間距
xx=pi/12;yy=pi/6;tt=pi/4;  %三個信號的入射角
% 8 陣元
A=[     1 exp(-j*2*pi*d*sin(xx)/l) exp(-j*2*2*pi*d*sin(xx)/l) exp(-j*2*3*pi*d*sin(xx)/l) exp(-j*2*4*pi*d*sin(xx)/l) exp(-j*2*5*pi*d*sin(xx)/l) exp(-j*2*6*pi*d*sin(xx)/l) exp(-j*2*7*pi*d*sin(xx)/l);...
        1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*2*pi*d*sin(yy)/l) exp(-j*2*3*pi*d*sin(yy)/l) exp(-j*2*4*pi*d*sin(yy)/l) exp(-j*2*5*pi*d*sin(yy)/l) exp(-j*2*6*pi*d*sin(yy)/l) exp(-j*2*7*pi*d*sin(yy)/l);...
        1 exp(-j*2*pi*d*sin(tt)/l) exp(-j*2*2*pi*d*sin(tt)/l) exp(-j*2*3*pi*d*sin(tt)/l) exp(-j*2*4*pi*d*sin(tt)/l) exp(-j*2*5*pi*d*sin(tt)/l) exp(-j*2*6*pi*d*sin(tt)/l) exp(-j*2*7*pi*d*sin(tt)/l)]'; %陣列流型

%%%%%%%%%%%% 
num=100;snr=5;N_gde=zeros(1,num);Pr=0;
for i=1:100
    %%%信號源
    N=40;Fd=1e4; Fc=1e5; Fs=1e6; sig=[]; 
    seq=randint(1,N,2);  sig(1,:)=dmod(seq,Fc,Fd,Fs,'fsk',2,3*Fd); 
    seq=randint(1,N,4);  sig(2,:)=dmod(seq,Fc,Fd,Fs,'ask',4); 
    t=[1:4000];
    seq=sin(2*pi*0.005*t);  sig(3,:)=amod(seq,Fc,Fs,'amdsb-tc',2); 
    
    sig=sig(:,1:1000);  % 取采樣數據長度
    x=awgn(A*sig,snr); 
    %%%%%%%%%%%%   GDE     %%%%%%%%%%%% 
    [nr,nc]=size(x);
    R=x*x'/nc;
    %%%%%%%%%%%%%%%%%%
    numR=size(R,1);
    R1=R(1:numR-1,1:numR-1);
    [V1 D1]=eig(R1);
    VV=[V1 zeros(numR-1,1); zeros(1,numR)];
    VV(numR,numR)=1;
    %%%%%%%%%%%%%%%%%%
    R2=VV'*R*VV;
    pp=R2(1:numR-1,numR);
    pp=abs(pp);
    pp=sort(pp);
    pp=pp(end:-1:1);
    %%%%%%%%%%%%%%%%%%
    for k=1:numR-1
        GDE(k)=pp(k)-0.2*mean(pp);  %%% D(N)   ????????
    end
    for k=2:numR-1
        if GDE(k)<=0 & GDE(k-1)>0
            N_gde(i)=k-1;
            if N_gde(i)==size(sig,1)
                Pr=Pr+1;
            end
        end
    end
end

GDE