%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reference
% Mike Novey and T. Adali, "On Extending the complex FastICA algorithm
% to noncircular sources" IEEE Trans. Signal Processing, vol. 56, no. 5, pp. 2148-2154, May 2008.
%
% Performs symmetric orthogonalization of the
% non-circular complex FastICA algorithm where
% xold is the mixtures and typeStr = 'log', 'kurt', or 'sqrt' is
% the nonlinearity.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [Ahat, shat] = nonCircComplexFastICAsym(xold,typeStr)
 type = 0;
 if strcmp(typeStr,'log') == 1
     type = 1;
 elseif strcmp(typeStr,'kurt') == 1
     type = 2;
 elseif strcmp(typeStr,'sqrt') == 1
     type = 3;
 end
 
 tol = 1e-5;
 a2=.05;
  defl = 1; % components are estimated one by one in a deflationary manner; set this to 0 if you want them all estimated simultaneously
 maxcounter = 50;
  [n,m] = size(xold);
% Whitening of s:
yyy = zeros(1,m);
[Ex, Dx] = eig(cov(xold'));
Q = mtimes(sqrt(inv(Dx)),Ex');
x = Q * xold;

%Pseudo-covariance 
pC = (x*transpose(x))/m;

% FIXED POINT ALGORITHM
 W = eye(n);
  Wold = zeros(n);
  k=0;
  
   while (norm(abs(Wold'*W)-eye(n),'fro')>(n*1e-5) && k < 15*n)
         k = k+1;
         Wold = W;
         
    for kk=1:n %Loop thru sources
        
        yy = W(:,kk)'*x;
       
        absy =abs(yy).^2; 
      %%Fixed point
      if type == 1 %%log
           g = 1./(a2 + absy);
           gp =  -1./(a2 + absy).^2;
      elseif type == 2  %Kurt
           g = absy;
           gp =  ones(size(absy));
      elseif type == 3  %sqrt
           g = 1./(2*sqrt(a2 + absy));
           gp =  -1./(4*(a2 + absy).^(3/2));
      end
       
       gRad =  mean(ones(n,1)*(g.*conj(yy)).*x,2);
       ggg = mean(gp.*absy + g);
       B = mean(gp .* conj(yy).^2)*pC;
       W(:,kk) =  Wold(:,kk)*ggg -(gRad) + (B*conj(Wold(:,kk)));
    
    end
    %Orthonormalization
    [E,D] = eig(W'*W);
    W = W * E * inv(sqrt(D)) * E';
   
  end; %Loop thru sources


shat = W'*x; %Estimated sources
Ahat = inv(Q)*W;