function B=BSS_NH_1(X,pram)
%This function does BSS based on Joint Diagonalization of the fourth order cumulant matrix slices
%of data. The JD scheme used is a Non-holonomic gradient based one introduced in:
%Bijan Afsari, P.S.Krishnaprasad "Some Gradient Based Joint Diagonalization methods for ICA", ICA2004.
%====================================================================================================================
%X=data of size n*T, n=#of sources=#sensors, T=sample size
%pram=[mu,tol]:
%mu=step size
%tol=error or stop measure which is bascially how close is gradient to zero (in this implementation tol is an absolute number
%in the sense that the real measure used will be tol*mu, see the while loop in the JD phase)
%==========================================================================================================================
%The values of above parameters shold be found based on trial an error and depending on the data and
%#of sources. (for small n: mu=.01,tol=.1 are typical values used)
%Note that we need small mu to keep the determinant of the updated matrix Bi1 close to 1 and we need large mu
%to speed up convergence
%============================================================================================================================
%B=overall un-mixing matrix matrix, which is equal to the product of the non-orthogonal 
%joint daigonalizer and the whitening matrix
%================================================================================================================================
%the steps in the algorithms are:
%1.first the data X is whitened thorugh a whitening matrix B1
%2.the 4th order cumulant matrix slices of the whitened data are found
%3.Using the non-holonomic gradient based method a joint diagonalizer Bi1 is found which has a determinant close to 1
%4.the overall unmixing matrix B=Bi1*B1 is formed
%=========================================================================================================================
%We warn that this implementation is not the fastest possible!!
%===========================================================================================================================
%In case of any comments and questions please contact: bijan@glue.umd.edu
%
%====================================================================================================================================
mu=pram(1);
tol=pram(2);
%=========================================================================
%Whitening step(this is not the best implementaion possible though!)
%===================================================================================
[n,T]=size(X);
X=X'-ones(T,1)*mean(X');
X=X';
R=X*X'/T;
[G,D]=eig(R);
B1=sqrt(inv(D))*inv(G);
X=B1*X;
%===================================================================================
%Computing the cumulant slices of the whitened data
%=====================================================================================
H0=[];
for i=1:n
   for j=1:n
   	   Cij=(repmat(X(i,:),n,1).*X)*(repmat(X(j,:),n,1).*X)'/T;
        Cij=Cij-(X*X'*(X(i,:)*X(j,:)'))/(T^2)-(X*X(i,:).'/T)*(X*X(j,:).'/T)'-(X*X(j,:)'/T)*(X(i,:)*X'/T);
      H0=[H0,Cij];
      
   end
end
%===========================================================================================
%H0 contains the cumulant matrix slices
%=======================================================================================
%=============================================================================================
%The non-orthogonal JD starts here, Bi1 will be result
%===========================================================================================

[n,N]=size(H0);
N=N/n;

Bi=eye(n,n);
K=zeros(n,n);
			for i=1:n:n*(N-1)+1
	                Hi=H0(:,i:i+n-1);
                 Gi=Bi*H0(:,i:i+n-1)*Bi';
                 K=K+(Gi-diag(diag(Gi)))*Gi;
       end
        grad=(K-diag(diag(K)))*Bi;
         Bi1=Bi-mu*grad;  


while norm(Bi1-Bi,'fro')/norm(Bi,'fro') > tol*mu,
   			Bi=Bi1;
            K=zeros(n,n);  
             for i=1:n:n*(N-1)+1
	                Hi=H0(:,i:i+n-1);
                 Gi=Bi*H0(:,i:i+n-1)*Bi';
                 K=K+(Gi-diag(diag(Gi)))*Gi;
              end
            grad=(K-diag(diag(K)))*Bi;%make the diagonal of K zero!
            Bi1=Bi-mu*grad;  
end%End of while
         
B=Bi1*B1;
