[Nb_source]=source_number_detection_hyp(received_signal,pfa,correction,display)

function [Nb_source]=source_number_detection_hyp(received_signal,pfa,correction,display)

%HELP source_number_detection_hyp
%
%Detection of the number of sources with hypothesis tests. The method
%assumes that the noise is spatially white.
%
%Input:   -received signal
%           -pfa: probability of false alarm (between 0 and 1)
%           -correction factor (0 or 1): 0 no correction, 1 correction
%           -(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_hyp(Y,0.005,1,1)
%
%Reference: [Law56] D. N. Lawley. Tests of signicance for the latent roots of covariance and
%           correlation matrices. Biometrika, 43(1/2) :128136, 1956.
%           [Jam69] A.T James. Tests of equality of latent roots of the covariance matrix.
%           Multivariate Analysis, pages 205217, 1969.

[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');

%Main Algorithm (See [Law56])

hyp_test=0;
for k=0:Nb_receivers-1
    %------------------------------%
    % Compute Likelihood function  %
    %------------------------------%
    coef=1/(Nb_receivers-k);
    a=coef*sum(eigenvalues(k+1:Nb_receivers));
    g=prod(eigenvalues(k+1:Nb_receivers)).^(coef); 
    likelihood(k+1)=-2*log((g/a)^(Nb_samples*(Nb_receivers-k)));    
    %If correction==1, we apply a correction term to the likelihood
    %function (see [Jam69])
    if correction==1
       correction_factor=1-(k/Nb_samples)-((2*(Nb_receivers-k)+1)/(6*Nb_receivers*(Nb_receivers-k)));
       likelihood(k+1)=correction_factor*likelihood(k+1); 
    end   
    %------------------------------%
    %   Perform Hypothesis test    %
    %------------------------------%
    %compute threshold                     
    degree_of_freedom=abs((Nb_receivers-k)^2-1);
    threshold(k+1)=chi2inv(1-pfa,degree_of_freedom);
    if (hyp_test==0) && ((likelihood(k+1) < threshold(k+1)) || (k==Nb_receivers-1))
            hyp_test=1;
            Nb_source=k;
    end
end

if nargin==4
   figure
   plot(0:Nb_receivers-1,likelihood);
   hold on;
   plot(0:Nb_receivers-1,threshold,'r');
   grid;
   xlabel('k');
   ylabel('Likelihood ratio');
   legend('Likelihood','threshold');
   plot(Nb_source,likelihood(Nb_source+1),'o');
  
end 

%Created by: Vincent Choqueuse, Ph.D, Website: http://www.vincentchoqueuse.com