急求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 mannerset 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，可以用。稍微注释了些

原始程序，不知道是谁写的了

FastICA算法的基本步骤：

1. 对观测数据进行中心化，使它的均值为0；

2. 对数据进行白化，。

3. 选择需要估计的分量的个数，设迭代次数

4. 选择一个初始权矢量（随机的）。

5. 令,非线性函数的选取见前文。

6. 。

7. 令。

8. 假如不收敛的话，返回第5步。

9．令，如果，返回第4步。

二．MATLAB源程序及说明：

%下程序为ICA的调用函数，输入为观察的信号，输出为解混后的信号

function Z=ICA(X)

%-----------去均值---------

[M,T] = size(X)%获取输入矩阵的行/列数，行数为观测数据的数目，列数为采样点数

average= mean(X')' %均值

for i=1:M

X(i,:)=X(i,:)-average(i)*ones(1,T)

end

%---------白化/球化------

Cx = cov(X',1) %计算协方差矩阵Cx

[eigvector,eigvalue] = eig(Cx)%计算Cx的特征值和特征向量

W=eigvalue^(-1/2)*eigvector' %白化矩阵

Z=W*X %正交矩阵

%----------迭代-------

Maxcount=10000 %最大迭代次数

Critical=0.00001 %判断是否收敛

m=M %需要估计的分量的个数

W=rand(m)

for n=1:m

WP=W(:,n) %初始权矢量（任意）

% Y=WP'*Z

% G=Y.^3%G为非线性函数，可取y^3等

% GG=3*Y.^2 %G的导数

count=0

LastWP=zeros(m,1)

W(:,n)=W(:,n)/norm(W(:,n))

while abs(WP-LastWP)&abs(WP+LastWP)>Critical

count=count+1 %迭代次数

LastWP=WP %上次迭代的值

% WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP

for i=1:m

WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i)

end

WPP=zeros(m,1)

for j=1:n-1

WPP=WPP+(WP'*W(:,j))*W(:,j)

end

WP=WP-WPP

WP=WP/(norm(WP))

if count==Maxcount

fprintf('未找到相应的信号)

return

end

end

W(:,n)=WP

end

Z=W'*Z

%以下为主程序，主要为原始信号的产生，观察信号和解混信号的作图

clear allclc

N=200n=1:N%N为采样点数

s1=2*sin(0.02*pi*n)%正弦信号

t=1:Ns2=2*square(100*t,50)%方波信号

a=linspace(1,-1,25)s3=2*[a,a,a,a,a,a,a,a]%锯齿信号

s4=rand(1,N)%随机噪声

S=[s1s2s3s4]%信号组成4*N

A=rand(4,4)

X=A*S%观察信号

%源信号波形图

figure(1)subplot(4,1,1)plot(s1)axis([0 N -5,5])title('源信号')

subplot(4,1,2)plot(s2)axis([0 N -5,5])

subplot(4,1,3)plot(s3)axis([0 N -5,5])

subplot(4,1,4)plot(s4)xlabel('Time/ms')

%观察信号（混合信号）波形图

figure(2)subplot(4,1,1)plot(X(1,:))title('观察信号（混合信号）')

subplot(4,1,2)plot(X(2,:))

subplot(4,1,3)plot(X(3,:))subplot(4,1,4)plot(X(4,:))

Z=ICA(X)

figure(3)subplot(4,1,1)plot(Z(1,:))title('解混后的信号')

subplot(4,1,2)plot(Z(2,:))

subplot(4,1,3)plot(Z(3,:))

subplot(4,1,4)plot(Z(4,:))xlabel('Time/ms')

一、定义：帝国竞争算法(imperialist competitive algorithm, ICA )是一种受帝国竞争行为启发的新的智能优化算法，它与粒子群优化(PSO)、蚁群(BCO)等算法一样，都属于基于群体的随机优化搜索算法。

二、详细介绍：

受帝国主义殖民竞争机制的启发，Atashpaz-Gargari和Lucas于2007年提出了一种新的智能优化算法—帝国竞争算法 (ICA)。与GA, PSO, ABC等受生物行为启发的群智能算法不同，ICA受社会行为启发，通过摸拟殖民地同化机制和帝国竞争机制而形成的一种优化方法。ICA也是一种基于群体的优化方法，其解空间由称为国家的个体组成。ICA将国家分为几个子群，称为帝国。在每个帝国内，ICA通过同化机制使非最优的国家(殖民地)向最优国家(帝国主义国家)靠近，该过程类似于PSO。帝国竞争机制是ICA的关键，ICA通过帝国竞争机制将最弱帝国中的一个或多个殖民地移动到其他帝国，使帝国之间可以进行信息交互。

目前，国外已有许多学者对ICA的性能改进以及实际应用进行了大量的研究，也取得了一定的进展。ICA已被广泛用于解决各种实际的优化问题，如调度问题、分类问题、机械设计等。然而，该算法仍然存在多样性下降较快、易早熟收敛等缺陷。另外，ICA提出的时间较短，尚有很大的研究空间。

三、优点：

帝国竞争算法与PSO, GA相比，收敛速度快、收敛精度高，具有较强的全局收敛性。算法利用殖民地向帝国主义国家移动进行局部搜索，即在较优区域内进行大力度开采，保证了算法的局部搜索能力。同时，帝国竞争 *** 作使帝国内的殖民地可以向其他帝国移动，突破了原来的搜索范围，增加了种群多样性，在一定程度上起到克服“早熟”现象的作用。此外，帝国合并 *** 作大大加快了算法的收敛速度，对于低维度优化问题，具有较明显的优势。

四、存在的问题：

群智能优化算法的“开采”和“勘探”能力是互相制约的，“开采”能力较强时，群体的多样性会受影响，而“勘探”能力较强则算法的全局收敛速度会变慢。原始的ICA算法还不能很好地平衡这两点，其局部搜索能力较强，收敛速度快，因此优化高维多模问题时，容易陷入局部最优。

帝国合并以及帝国覆灭使ICA的帝国个数不断减少，导致群体多样性降低，算法的全局“勘探”能力受影响，易出现“早熟”现象。

帝国竞争 *** 作体现了帝国之间的信息交互，然而，帝国竞争在每一次迭代中只是将最弱的殖民地归于最强的帝国，该过程对每个帝国的势力大小影响很小，需要多次迭代才能体现出来，帝国之间缺乏更有效的信息交互，即群体多样性的体现并不明显。

---