function [ea esnr eainit] = getA(X,threshold,maxiter)

% --function ea = getA(X, threshold, maxiter)
%  Estimatie the mixing vector a from data matrix X.
%  
%  output -
%          ea:     estimated a (ea.^2 is estimated power)
%          esnr:   estimated SNR
%          eainit: estimated a before ALS
%  input -
%          X:         n-by-T data matrix, n:channel, T:sample size
%          threshold: Threshold of the alternating least square 
%          maxiter:   Maximum number of iteration in ALS.
%
% Author: Yoshikazu Washizawa (washizawa@bsp.brain.riken.jp)
%          Lab. for Advansed Brain Signal Procesing,
%          Brain Science Institute, RIKEN 
%
%

if nargin == 1
  threshold=10^(-6);
  maxiter=1000;
elseif nargin ~= 3
  error('Usage: ea = getA(X,threshold,maxiter)');
end

n=size(X,1);
T=size(X,2);

R=X*X'/(T-1);
R(1:n+1:n^2)=0;
ea=zeros(n,1);

%%% Get initial value
[tmp b]=sort(reshape(R,n^2,1));
sindex=sort([ceil(b/n) mod(b-1,n)+1]')';

t=sindex(end,1);
u=sindex(end,2);
indexset=setdiff(1:n,[t u]);
for s=indexset
  ea(s) = sqrt(R(s,t)* R(u,s)  / R(t,u) );
end

p=max(find(sindex(:,1)~=t & sindex(:,2)~=t));
s2=t;
t2=sindex(p,1);
u2=sindex(p,2);
ea(s2) = sqrt(R(s2,t2)*R(u2,s2)/ R(t2,u2));

s3=u;
p=max(find(sindex(:,1)~=u & sindex(:,2)~=u));
t3=sindex(p,1);
u3=sindex(p,2);
ea(s3) = sqrt(R(s3,t3)*R(u3,s3)  / R(t3,u3) );
ea=real(ea);

eainit=ea;
 
%%% Alternative Least Square 
% If you have complied mex file "ALS_mex", replace `1' to `0.'
if 1
  oldvec=zeros(n,1);
  iter=1;
  while (norm(oldvec - ea) > threshold) && (iter<maxiter)
    oldvec=ea;
    for ii=1:n
      index=setdiff(1:n,ii);
      ea(ii)=(ea(index)'*R(index,ii))/norm(ea(index))^2;
    end
    iter=iter+1;
  end
  
  if iter == maxiter
    warning('maxiter reached');
  end
else
  ea = ALS_mex(ea,R,threshold,maxiter);
end

if sum(isnan(ea)) > 0
  ea=eainit;
  warning('ea is inifnite');
end

esnr=abs( (ea.^2) ./ ( var(X')'-(ea.^2)));