function [l p] = glrtcirc(X)
% [l p] = glrtcirc(X,alpha)
% Computes the Generalized Likelihood Ratio Test (GLRT) of circularity 
%
% The progrman computes the test statistic and the p-value of the GLRT of 
% circularity (Ollila and Koivunen, 2004, 2009). Computation of the
% p-value requires STATISTICS toolbox (MathWorks).   
%
% INPUT:
% X                = data matrix with observations as rows
%
% OUTPUT 
% l                = the GLRT statistic  
% p                = p-value 
%
% REFERENCES.  See followin papers and the references therein: 
% Ollila and Koivunen (2004), "Generalized complex elliptical distributions", 
% in Proc. 3rd IEEE Sensor Array & Multichannel Signal Proc. Workshop(SAM'04), 
% pp. 460--464, Sitges, Spain, July 18-21, 2004.
%
% Ollila and Koivunen (2009), "Adjusting the generalized likelihood ratio 
% test of circularity robust to non-normality", in Proc. 10th IEEE Int'l
% Workshop on Signal Proc.Advances in Wireless Communications (SPAWC'09), 
% Perugia, Italy, June 21-24, 2009. 
%
% Version 1.0, June 2009, (c) Esa Ollila (esollila@wooster.hut.fi), 
% Signal Processing Laboratory, Helsinki University of Technology
% The program was developed under Matlab 7.6.0.32
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n k] = size(X);

if n <= k  
   error('glrtcirc: input data matrix has more variables (columns) as observations (rows)');
end

% Center the data
X = X - repmat(mean(X),n,1); 
C = cov0(X);   % Covariance matrix 
P = pcov0(X);  % Pseudo-covariance matrix 

if k > 1, 
   iC = inv(C); 
   % Gets rid of infinitesimal imaginary parts on the diagonal 
   iC(1:(k+1):k^2)=real(iC(1:(k+1):k^2)); 
   hiC = iC';
   % We force the matrix to be Hermitian (due to finite sample precision, this may
   % not always be the case) 
   iC(repmat(1:k,k,1) > repmat(1:k,k,1)')= hiC(repmat(1:k,k,1) > repmat(1:k,k,1)'); 
else 
    iC = 1/C;
end
R = iC*P;
glr = real(det(eye(k)-R*conj(R)));  % eq. (7) in Ollila and Koivunen (2009)
l = -(n-k)*log(glr);                % GLRT test statistic
if license('test','statistics_toolbox')
   p = 1- chi2cdf(l,k*(k+1));          % p-value   
else 
   p = [ ] ;
end



