function [X,whiteningMatrix, dewhiteningMatrix] = whitening(Y,rdim)

% Whitening data X. 

if(nargin<2)
    rdim = size(Y,1);
end

%Calculate the eigenvalues and eigenvectors of covariance matrix.
% Here, assuming that means of each rows in Y is are zero.
covMatrix = Y*Y'/size(Y,2);
[E, D]= eig(covMatrix);

%E: the matrix are composed by eigenvectors as the columns
%D: the diagonal components of this matrix are eigenvalue 

[~,order] = sort(diag(D),'descend');

E= E(:,order(1:rdim));

d=diag(D);
d=real(d.^(-0.5));
D=diag(d(order(1:rdim)));

X= D*E'*Y; 

whiteningMatrix = D*E';
dewhiteningMatrix = E*D^(-1);