急求FastICA 的源程序 matlab，包括数据的预处理（中心化和白化），注释详细点，谢谢！

急求FastICA 的源程序 matlab，包括数据的预处理（中心化和白化），注释详细点，谢谢！

% function [Ahat2, shat, n_iteration Test] = nc_fastica_svd(xold,typeStr,N,A)
function [shat Ahat2] = nc_fastica_svd(xold,typeStr,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% non-circular complex FastICA算法，基于Newton迭代法，类似与fastICA
% ************************input***************************
% xold： 混合信号，m*n,m为阵元数，n为快拍数
% typeStr： 非线性函数，'log', 'kurt', or 'sqrt'
% **************************output**************************
% Ahat： 解混矩阵
% shat: 估计的源信号
% ********************************************************
% Reference
% Mike Novey and T. Adali, "On Extending the complex FastICA algorithm
% to noncircular sources" in
% (To appear 2007/2008) IEEE Journel on Signal Processing.,
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 = 0.1;
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] = svd(cov(xold'));
E = Ex(:,1:N);
D = Dx(1:N,1:N);
Q = mtimes(sqrt(inv(D)),E');
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-12) && 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
n_iteration = k;
shat = W'*x; %Estimated sources
% Ahat1 = inv(Q)*W;
Ahat2 = W'*Q;

这个是NC-fastica，可以用。稍微注释了些
原始程序，不知道是谁写的了