function [C, iter] = chubM(z, q);
% [C, iter] = chubM(z, q);
% Computes the complex Huber's M-estimator of scatter
%         
% INPUT:
% z    : data matrix with n rows (observations) and k columns (variables).
% q    : user defined constant. A default choice is q=0.9.
% OUTPUT:
% C    : M(q)-estimator of scatter 
%         
% Esa Ollila, esollila@wooster.hut.fi 


if (nargin == 1)
   q=0.90;
end

[n k] = size(z);

csq = chi2inv(q,2*k)/2;  % cutoff-point c^2
b = chi2cdf(2*csq,2*(k+1))+(csq/k)*(1-q); % consistency factor

MAX_ITER = 250; % Max number of iteration
EPS = 1.0e-4;   % Iteration accuracy

invC0 = inv(z.'*conj(z)/n); % inverse of the initial estimate 
iter = 1;
while (iter<MAX_ITER)     
      s = real(sum(conj(z)*invC0.*z,2));
      w = (1/b).*(s<=csq) + ((csq/b)./s).*(s>csq);
      C = z.'*(conj(z).*repmat(w,1,k))/n;
      d = norm(eye(k)-invC0*C,1);
      if (d<=EPS) break; end
      invC0 = inv(C);
      iter = iter+1;     
end      

C(1:(k+1):(k^2)) = real(C(1:(k+1):(k^2)));

if (iter==MAX_ITER) 
     fprintf(1,'WARNING! Slow convergence: the error of the solution is %f\n',d);  
     fprintf(1,'after %d iterations\n',iter);
end      
           
   
