function [nb_sources]=source_number_detection_IC(received_signal,criterion,display)
%HELP source_number_detection_IC
%
%Detection of the number of sources with information criteria. 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_IC(Y,'MDL',1)
%
%Reference: [WAX85] Wax, M. and Kailath, T., "Detection of signals by information
%theoretic criteria", IEEE Transactions on Acoustics, Speech and Signal
%Processing, 1985.
[Nb_receivers,Nb_samples]=size(received_signal);
%Computation of the covariance matrix
covariance=received_signal*received_signal'/Nb_samples;
%eigenvalue decomposition
[U,V]=eig(covariance);
%sort the eigenvalue in the decreasing order
eigenvalues=sort(diag(V),'descend');
%Compute criterion
for k=0:Nb_receivers-1
coef=1/(Nb_receivers-k);
a=coef*sum(eigenvalues(k+1:Nb_receivers));
g=prod(eigenvalues(k+1:Nb_receivers)).^(coef);
akaike_criterion(k+1)=-log(((g/a)^(Nb_samples*(Nb_receivers-k))))+k*(2*Nb_receivers-k);
MDL_criterion(k+1)=-log(((g/a)^(Nb_samples*(Nb_receivers-k))))+0.5*k*(2*Nb_receivers-k)*log(Nb_samples);
end
switch criterion
case 'AIC'
criterion_value=akaike_criterion;
case 'MDL'
criterion_value=MDL_criterion;
end
[criterion_value_min,nb_sources]=min(criterion_value);
nb_sources=nb_sources-1; %retrieve 1 because the value k=0 is stocked in the index 1.
if nargin > 2
plot(0:Nb_receivers-1,criterion_value);
grid;
hold on
plot(nb_sources,criterion_value_min,'ro');
hold off
ylabel('Criterion value')
xlabel('Number of sources')
end
%Created by: Vincent Choqueuse, PhD (vincent.choqueuse@gmail.com)