fastica.m

来自「FastICA算法」· M 代码 · 共 95 行

clear all
clc

%观察信号矩阵 X 为 K*N 矩阵
%源信号矩阵 S 为 K*N 矩阵
%混合矩阵 A 为 K*K 矩阵

K=4;
N=157;

%生成 K 个源信号, 存储到 S
t=1:N;
figure(1);

S(1,:)=rand(1,N);
subplot(2,2,1);
plot(t,S(1,:));
title('随机信号');
hold on;

S(2,:)=square(2*pi*t/64);  
subplot(2,2,2);
plot(t,S(2,:));
title('方波信号');
hold on;

S(3,:)=sin(2*pi*t/32);  
subplot(2,2,3);
plot(t,S(3,:));
title('正弦信号');
hold on;

S(4,:)=cos(2*pi*t/16); 
subplot(2,2,4);
plot(t,S(4,:));
title('余弦信号');
hold on;

%生成 K 个混合信号,存储到 X, A 中存储混合矩阵
A=randn(K);
X=A*S;
figure(2);
plot(t,X);
title('观察信号');

% X 数据中心化
m=mean(X,2);
for i=1:N
    X(:,i)=X(:,i)-m
end

% X 数据白化
covMat=cov(X');
[E,D]=eig(covMat);
V=E*D^(-0.5)*E';
X=V*X;
figure(3);
plot(t,X);
title('白化后的信号');

%FastIca算法
%求完后 W 是 混合矩阵 A 的近似
W=rand(K);
for p=1:K
    W(:,p)=W(:,p)/norm(W(:,p));
    e=1
    while e<800
        pred=W(:,p);
        
        %定点法
        W(:,p)=1/N*X*((W(:,p)'*X).^3)'-3*W(:,p);

        %采用收缩策略，类似于正交化的过程
        sum=zeros(K,1);
        for i=1:p-1
            sum=sum+W(:,p)'*W(:,i)*W(:,i);
        end
        W(:,p)=W(:,p)-sum;
        W(:,p)=W(:,p)/norm(W(:,p));
        
        %收敛条件
        e=e+1
    end
end

%推算源信号
out=W'*X

figure(4);
for k=1:K  
    subplot(2,2,k);
    plot(t,out(k,:));
    hold on;
end;