[nb_source]=source_number_detection_radii(received_signal,criterion,display)

function [nb_source]=source_number_detection_radii(received_signal,criterion,display)

%HELP source_number_detection_radii
%
%Detection of the number of sources with gergoriin radii disks. The method
%assumes that the noise is spatially white.
%
%Input:     -received signal
%           -criterion ('AIC' or 'MDL')
%           -(optional) Display threshold if non-empty
%
%Ouput:     -nb_sources
%
%Example:   X=randn(3,500); %3 sources
%           H=randn(10,3)+i*randn(10,3); %10 receivers
%           B=2*(randn(10,500)+i*randn(10,500));
%           Y=H*X+B; %MIMO Mixing Model
%           source_number_detection_radii(Y,'MDL',1)
%
%Reference: [WYC95] H.T. Wu, J.F. Yang, and F.K Chen. Source number estimators using
%           transformed gerschgorin radii. IEEE Transactions on Signal Processing,
%           43(6):1325-1333, 1995.

[Nb_receivers,Nb_samples]=size(received_signal);

%compute the covariance matrix
covariance=(received_signal*received_signal')/Nb_samples;
%perform eigenvalue decomposition
[U,V]=eig(covariance);
eigenvalues=sort(diag(V),'descend');

%Modified covariance
cll=covariance(end,end);
C1=covariance(1:Nb_receivers-1,1:Nb_receivers-1);
[U1,D1]=eig(C1);
[value_p_order,index_order]=sort(diag(D1),'descend');
U1=U1(:,index_order);
U=[U1 zeros(Nb_receivers-1,1);zeros(1,Nb_receivers-1) 1];
modified_covariance=U'*covariance*U;

modified_center=abs(diag(modified_covariance));
modfied_covariance_without_diag=modified_covariance-diag(diag(modified_covariance));
modified_radii=sum(abs(modfied_covariance_without_diag),2);

%calcul de la fonction de vraisemblance
for k=0:Nb_receivers-1
    coef_alpha=(prod(modified_center(k+1:Nb_receivers-1)))^(1/(Nb_receivers-k-1));
    coef_beta=(1/(Nb_receivers-k-1))*sum(modified_center(k+1:Nb_receivers-1));
    coef_radii=sum((modified_radii(1:k).^2)./(modified_center(1:k)));
    
    likelihood(k+1)=Nb_samples*(Nb_receivers-1-k)*log(coef_alpha/coef_beta)...
                   -Nb_samples*log(cll-coef_radii);
         
    %Penaly factor          
    %Remark: There is a problem in the original publication on equation
    %38. GLE(k)=-phi(k)+Penalit !!!!! (and not + phi(k)+penalite)
    
    switch criterion
        case 'AIC'
            likelihood(k+1)=-likelihood(k+1)+(k^2+k);
        case 'MDL'
            likelihood(k+1)=-likelihood(k+1)+0.5*(k^2+k)*log(Nb_samples);
    end
end

[maximum_lik,nb_source]=min(likelihood);
nb_source=nb_source-1;    

if nargin==3    %affichage des figures   
    figure
    hold on;
    for k=1:Nb_receivers-1
        circle=rsmak('circle',modified_radii(k),[modified_center(k) 0]);
        fnplt(circle);
        plot(modified_center(k),0,'ro');
        plot(value_p_order,zeros(1,length(value_p_order)),'k'),
    end    
    grid
    hold off;
    figure
    hold on;
    plot(0:Nb_receivers-1,likelihood);
    plot(nb_source,likelihood(nb_source+1),'o');
    grid;
    hold off;
    
end