ica.m

用于脑电信号特征提取的InfoMax Algorithm Based on ICA；也可以稍作改动用于其他信息提取。

%InfoMax Algorithm Based on ICA
%(C) All Copyrights Reserved By Yuheng Wang
%2005.5.21 revised
%--------------------------------------变量含义-----------------------------------
%m-------采样通道数
%n-------采样点数
%A-------变换矩阵
%W-------解变换矩阵
%I-------单位矩阵
%U-------仿真结果数据
%step_max-------训练最大步数
%step-------训练步数
%X-------训练观测数据
%S-------已知理想源数据
%lrate-------学习速率
%u-------学习系数
%errormax-------最大误差
%annealrata------------退火系数
%K4-------区超高斯和亚高斯信号的方阵
%NLU-------U的非线性矩阵
%t-------输出的独立分量数
%**************************************可执行代码*********************************
%--------------------------------------初始化数据-------------------------------
m=4;
n=2000;
t=4;
W=rand(t,m);                                                           
I=eye(m);
step_max=500;
lrate=0.0005;
u=0.0002;
K4=eye(m);
errormax=0.00001;
degconst = 180./pi;
%---------------------------------------------------------------------------载入训练的观测数据   X1:C3--01;X2:
load X.mat;                                          
%--------------------------------------画观测波形---------------------------
figure                                        
for k=1:1:4
    subplot(4,1,k);
    j=1:1:n;
    plot(j,X(k,j));
end
%---------------------------------------------------------------------------流程控制
button = questdlg('观测波形绘制完毕！现在开始训练解变换矩阵吗?',...  
'提示信息:','Yes','No','Yes');
if strcmp(button,'No')
   break;
end
%--------------------------------------训练变换矩阵--------------------------
 for step=1:1:step_max
%----------------------------------------------------------------------------计算信号峰度以区别高斯和超高斯信号,curt是求解峰度的函数
     for i=1:1:4                                                                                                                                                                                        
     K4(i,i)=kurt(U(i,:));    
     end
%----------------------------------------------------------------------------计算跌代差量   
     NLU=tanh(U);
     dW=(I-K4*NLU*U.'-U*U.')*W ;
     dW=u*dW;
%----------------------------------------------------------------------------dW归一化     
     mo=0;
     mo=sqrt(sum(sum(dW.^2)));                                               
     dW=dW./mo;
%----------------------------------------------------------------------------判断W是否收敛
     s=0;
     for i=1:1:4
         for j=1:1:4
             if (abs(dW(i,j))<=errormax)   s=s+1;                            
             else break;
             end
         end
         if (j<4) break;
         end
     end
     
     if(s==16)   break;
     end
                                                                             
      %if degconst*angledelta > annealdeg,                                    %用退火法调整学习速率
        %lrate = lrate*annealstep;          % anneal learning rate
        %olddelta   = delta;                % accumulate angledelta until
        %oldchange  = change;               %  annealdeg is reached
      %elseif step == 1                     % on first step only
        %olddelta   = delta;                % initialize 
        %oldchange  = change;               
      %end
      
     dW=lrate.*dW;
     W=W+dW;
 end    
%-------------------------------------------------------------------------------流程控制                                                                           
 button = questdlg('训练完毕！现在要仿真变换矩阵吗?',...
'提示信息:','Yes','No','No');
if strcmp(button,'No')
   break;
end

%--------------------------------------提取特征节律------------------------------
U=W*X;
%--------------------------------------------------------------------------------画仿真图形
figure                                       
for k=1:1:4
    subplot(4,1,k);
    j=1:1:n;
    plot(j,U(k,j));
end
%--------------------------------------------------------------------------------流程控制
button = questdlg('仿真结果波形绘制完毕！现在要绘制理想信号源波形吗?',...
'提示信息:','Yes','No','Yes');
if strcmp(button,'No')
   break;
end
%--------------------------------------画信号源的理想输出波形---------------------
%---------------------------------------------------------------------------------以下波形分别是典型的αβθδ波：
load Y.mat;
figure                                                                  
for k=1:1:4
    subplot(4,1,k);
    j=1:1:n;
    plot(j,Y(k,j));
end
%--------------------------------------解变换矩阵------------------------------
A=inv(W)