[nb_source]=source_number_detection_PET(received_signal,display)

function [nb_source]=source_number_detection_PET(received_signal,display)

%HELP source_number_detection_PET
%
%Detection of the number of sources with predicted eigenvalue threshold. The method
%assumes that the noise is spatially white.
%
%Input:     received_signal
%           criterion ('Akaike' or 'MDL')
%           (optional) Display criterion if non empty
%Sorties:   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_PET(Y,1)
%
%Reference: [CHE91] Chen. W and Wong.K.M. and Reilly. J.P. "Detection of the 
%number of signals: a predicted eigen-thresholdapproach", IEEE Transactions 
% on Signal Processing, 1991.


T=1.5;

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

%main algorithm (see [CHE91])
test=0;
nb_source=1;
nb_noise_eigenvalues=1;
threshold_save=zeros(1,Nb_receivers-1);
while test==0 && (nb_noise_eigenvalues ~= Nb_receivers);
   mean_eigenvalue=mean(eigenvalues(end-nb_noise_eigenvalues+1:end));
   threshold= ((nb_noise_eigenvalues+1)*...
          ((1+T*(Nb_samples*(nb_noise_eigenvalues+1))^(-0.5))/...
          (1-T*((Nb_samples*nb_noise_eigenvalues)^(-0.5))))...
          -nb_noise_eigenvalues)*mean_eigenvalue;
      
   threshold_save(1:end-nb_noise_eigenvalues+1)=threshold;   
   if eigenvalues(end-nb_noise_eigenvalues) <= threshold
       nb_noise_eigenvalues=nb_noise_eigenvalues+1;
   else
       test=1;
   end 
end    

nb_source=Nb_receivers-nb_noise_eigenvalues;

if nargin==2;
    stem(eigenvalues);
    hold on;
    plot(threshold_save,'r');
    hold off;
    grid;
    axis([0.5 Nb_receivers+0.5 -0.5 max(eigenvalues)+0.5]);
    xlabel('k');
    ylabel('eigenvalue');
end    
    
%Created by: Vincent Choqueuse, PhD ([e-mail])