独立成分分析 matlab 程序

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')

程序提示的错误是函数hua_fft()的输入参数个数不足。

你在定义函数的时候

function hua_fft(y,fs,style,varargin)

表明里面有3个固定参数，还有可变的参数，因此函数在调用的时候至少得3个输入参数。

但是，你在函数调用的时候，却是2个输入参数

subplot(212)hua_fft(y,fs)title('原始信号频谱')

并没有传递style参数的值，因此才会报错说找不到style。

MATALB属于半开源软件，其中很多函数可以通过“open/edit/type+filename”命令进行编辑和查看源代码。但是有一些函数仅仅可以找到它的帮助文档，却无法看到具体的源代码，比如min,fft,sum等函数，因为这些函数属于MATLAB的built-in function（内置函数），即MATLAB的built-in function的代码是不公开的，有人说这些函数的算法是最优化的，保证较低的时间复杂度提高效率。

要确定哪些函数可以看代码哪些函数不可以看代码，可以通过which命令：

上图中的函数fft, sum, min属于MATLAB的built-in function，无法看代码，只能通过“help/type+filename”来看帮助文档，通过这个路径一路查下去，M文件里面也只有帮助文档，没有源代码。最后一个函数axis.m就不是built-in函数，可以通过 “edit/open/type +filename”来查看源代码。

---