informax_fastica.m

用matlab实现扩展Informax算法和固定点算法Fast ICA

clear all
clc

t=0:1/100:3;
s=zeros(3,length(t));

s(1,:)=sin(40*t);
s(2,:)=0.3*sawtooth(10*t);
s(3,:)=0.5*square(20*t);

figure(1),
for i=1:3,
    subplot(3,1,i),plot(t,s(i,:));
    axis([t(1) t(length(t)) min(s(1,:)) max(s(1,:))]),title('原始的信号');
end;


A=[0.3062 0.4668 0.7241;
   0.1122 0.0147 0.2816;
   0.4430 0.6641 0.2618];

s=A*s;
s1=s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%去均值
s=s-mean(s')'*ones(1,length(t));
s2=s

%白化处理
CX=cov(s',1);
[E namda]=eig(CX);
ss=E*inv(sqrtm(namda))*E'*s;    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%扩展Informax算法
s_informax=ss;

w= randn(3);
step=0.1; 

u=w*s_informax;

delta=0;
delta_w=0;
temp1=0;
while((delta_w==0)|(max(max(abs(delta_w))<=10^(-6)))),     
    temp2=temp1;
    temp1=delta_w;    
    delta_w=step*(eye(3)+tanh(u)*u'-u*u')*w;
    if temp2<temp1&temp1<delta_w,
        disp(['算法不收敛']);
        break
    end;
    w=w+delta_w;
    u=w*s_informax;
end

y_informax=w*ss;
figure(2),
for i=1:3,
    subplot(3,1,i),plot(t,y_informax(i,:)),title('恢复的信号');
end;



%固定点算法Fast ICA
s_ICA=ss;

ipslong=10^(-5); 

w_ICA=zeros(3,1);
for i=1:3,
    w_ICA(:,i)=rand(3,1); 
    w_ICA(:,i)=w_ICA(:,i)/norm(w_ICA(:,i));
    
    stop=0;
    
    while (stop==0),
        temp=w_ICA(:,i);
        for j=1:3,
            w_ICA(j,i) = mean(s_ICA(j,:).*(tanh((temp)'*s_ICA)))-(mean(1-(tanh((temp))'*s_ICA).^2)).*temp(j,1);
        end;
        sum=zeros(3,1);
        for k=1:i-1,
            sum=sum+(w_ICA(:,i)'*w_ICA(:,k))*w_ICA(:,k);
        end;
        w_ICA(:,i)=w_ICA(:,i)-sum;
        w_ICA(:,i)=w_ICA(:,i)/norm(w_ICA(:,i));
        if(abs(dot(w_ICA(:,i),temp))<1+ipslong)&(abs((dot(w_ICA(:,i),temp)))>1-ipslong),
            stop=1;
        end;
    end;
    
end;
w_ICA=w_ICA';

y_ICA=w_ICA*s_ICA;
figure(3),
for i=1:3,
    subplot(3,1,i),plot(t,y_ICA(i,:)),title('恢复的信号');
end;