function B=BSS_LU_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 the LU 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,itr]:
%mu=step size(is the same for L and U part)
%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 loops in the L and U phases)
%itr= #of times the L-U phase is repeated
%===============================================================================
%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,itr=3 are typical values used)
%================================================================================
%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 LU gradient method a joint diagonalizer WW is found (note that by construction det(WW)=1)
%4.the overall unmixing matrix B=WW*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);
itr=pram(3);
%=========================================================================
%Whitening step(this is not the best implementaion possible though!)
%===================================================================================
[n,T]=size(X);
X=X'-ones(T,1)*mean(X');
X=X';
%Whitening
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
%the L-U phases are repeated 'itr' times
%============================================
 [n,N]=size(H0);
N=N/n;        
WW=eye(n,n);
j2=0;
for I=1:itr
Li=zeros(n,n);
Li1=eye(n,n);
%==============================================
%The L(=lower triangualr) phase 
%==============================================
while norm(Li-Li1,'fro')/norm(Li1,'fro') > mu*tol,
   			Li=Li1;
 				K=zeros(n,n);  
             for i=1:n:n*(N-1)+1
	                Hi=H0(:,i:i+n-1);
                 Gi=Li*H0(:,i:i+n-1)*Li';
                 K=K+(Gi-diag(diag(Gi)))*Gi;
              end
            grad=(tril(K-diag(diag(K))))*Li;
            Li1=Li-mu*grad;
            j2=j2+1;
         end
         %=================
         %The L part is computed and the cumulant matrices are updated
         %==================
  for i=1:n:n*(N-1)+1
   	H0(:,i:i+n-1)=Li1*H0(:,i:i+n-1)*Li1';
	end
   %================================
   %The U phase starts here
   %======================================
Ui=zeros(n,n);
Ui1=eye(n,n);
while norm(Ui-Ui1,'fro')/norm(Ui1,'fro') > mu*tol,
   			Ui=Ui1;
 				K=zeros(n,n);  
             for i=1:n:n*(N-1)+1
	                Hi=H0(:,i:i+n-1);
                 Gi=Ui*H0(:,i:i+n-1)*Ui';
                 K=K+(Gi-diag(diag(Gi)))*Gi;
              end

            grad=(triu(K-diag(diag(K))))*Ui;
            Ui1=Ui-mu*grad;
            j2=j2+1;
         end
         WW=Ui1*Li1*WW;
   				for i=1:n:n*(N-1)+1
   					H0(:,i:i+n-1)=Ui1*H0(:,i:i+n-1)*Ui1';
					end
			 
  end%end of loop over L-U phases
B=WW*B1;%the overall un-mixing matrix is computed