function [Z, W, IW] = sphere(X)
% SPHERE: transforms the original data such that the mean is 0
% and the covariance is identity (de-correlates the variables).
%
% [Z,W,IW] = sphere(X)
%
% X - [n x m] matrix of n observations, each with m variables
% Z - [n x m] the transformed matrix
% W - [m x m] transformation matrix: Z = X * W
% IW- [m x m] the inverse transformation: X = Z * W
%
% If the covariance matrix of X is rank deficient (e.g. less
% samples than variables), the transformed matrix Z may have 
% less variables.

my_eps = 1e-8;

[n, m] = size(X); 
% Instead of:
% Sg = cov(X);             
% compute directly the covariance (sample-based):
Xm = X - repmat(mean(X,1), [n,1]);
Sg = Xm'*Xm / (n-1);  % Sigma = E[X * X']

[V, L] = eig(Sg);        % eigenvectors by columns in V

d = diag(L);
d(d <=my_eps) = 0;

if sum(d > 0) < m
    warning('Precision lost in transformation!');
end

d = 1./sqrt(d);

W = V * diag(d); % by formula: (diag(d) * V')';
Z = X * W;

IW = sqrt(L) * V'; % = inv(W) = inv(V*diag(d))=...

return
