X=randn(100,3);
[mx, nx] = size(X);
M = sum(X, 1)/mx; % centroid (mean)
Xc = X - M(ones(mx,1),:); % subtract centroid of X
W = (Xc' * Xc)/(mx - 1); % variance matrix
% The call to REAL here is only to remove ``numerical noise''.
D = real(sum((Xc / W) .* conj(Xc), 2)); % Mahalanobis distances