function [y, W,lambda] = LDA(X,classes)

    % gloabal mean of data
    globalMean = mean(X)';
    M = size(X,2);
    
    Sw = zeros(M,M);
    Sb = zeros(M,M);
    for i = [unique(classes)]'
        % find elements of each class
        index = find(classes == i);
        % compute the mean for class i 
        mui = mean(X(index,:))';
               
        Si= zeros(M,M);
        % go througth all elements of a class
        for j = 1:length(index)
            x = X(index(j),:)';
            Si = Si+ (x-mui)*(x-mui)';
        end
        Sw = Sw + Si;
        Sb = Sb + ((mui-globalMean)*(mui-globalMean)').*length(index); 
        
    end

    %%pseudo-inverse of Sw
    AAA = pinv(Sw)*Sb;
    
   % W: eigen-vectors	
   % lambda : eigen-values
   [  W , lambda ] = eig ( AAA);

   
   W = real(W);
   lambda = real(diag(lambda));
   
   %perform projection
   % sorts eigen vector w.r.t eigen values
    [ss,ii] = sort(lambda,1,'descend');
    W=  W(:,ii(1: min(length(unique(classes))-1,size(W,2))  ));
    y = X * W;
    
    
end