
function [NewData,WhitenMatrix,DewhitenMatrix,Vsort, dsort] = GIG_Pca(data,numpc)
%%   
% compute covariance matrix by removing the mean and then calculate the
% eigen values and eigen vectors,then PCA.

% Input: 1. data - 2D array whose dimension is voxels by time points
% 2.numpc: retained number of components 
% Output:
% NewData:whitened data
% written by duyuhui, CAS,China,2011
    %% Calculate covariance matrix
    
    CovMat = GIG_Cov(data);
   
    %% Calculate the eigenvalues and eigenvectors of covariance matrix.
    disp('Calculating eigendecomposition');
    [V, Lambda] = eig(CovMat, 'nobalance'); %Lambda 

    %% Sort the eigenvalues - decending.
    disp('Sorting eigenvalues');
    [eigenvalues index] = sort(diag(Lambda)); %sort����Ľ��������
    
    cols=size(V,2);%  ?
    Vsort=zeros(size(V));
    dsort=zeros(size(eigenvalues));
    for i=1:cols
        Vsort(:,i) = V(:, index(cols-i+1) );
        dsort(i)   = eigenvalues(index(cols-i+1) ); 
    end   
    dsum = sum(dsort);
%     dsum_extract1 = 0;
%     p = 0;
%     while( dsum_extract1/dsum < 0.9)
%         p = p + 1;
%         dsum_extract1 = sum(dsort(1:p));
%     end
%      retained1=dsum_extract1/dsum;

    
    dsum_extract2 = sum(dsort(1:numpc));
    retained2=(dsum_extract2/dsum)*100;
%     fprintf('%g%% of (non-zero) eigenvalues retained.\n', retained2);
    Vpart=Vsort(:,1:numpc);
    dpart=dsort(1:numpc);
    Lambda_part=diag(dpart);
    WhitenMatrix = sqrtm(Lambda_part) \ Vpart'; % Use gaussian elimination approach to solve the equations
    DewhitenMatrix = Vpart * sqrtm(Lambda_part);

    fprintf('Whitening and PCA\n');

    %newVectors =  whiteningMatrix*data;
    NewData = data * WhitenMatrix';
  
    
    
