function B = bell(X)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Bell & Sejnowsky のアルゴリズムで混合データ x を分離する行列 %%%%%
%%%%%%                                         w を求めるプログラム %%%%%
%%%%%%                前処理の部分はCardosoのプログラムから一部拝借 %%%%%
%%%%%%                                                     樺島祥介 %%%%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% データ行列のサイズを調べる
[N,P]=size(X);

% 定行列の設定
Id=eye(N);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 前処理 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 各行の平均値を引く %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MX=mean(X'); 
X=X-mean(X')'*ones(1,P);                   

%%%%%%%%%%%%%%%%%% ホワイトニング %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 [U,D] 		= eig((X*X')/P)	; % 共分散行列の固有値と固有ベクトル
                                  % U: 固有ベクトル D: 固有値
 [puiss,k]	= sort(diag(D))	; % 固有値の大きさ順に並べ替え
 rangeW		= 1:N	; 
 scales		= sqrt(puiss(rangeW)); % 固有値の平方根
  W  		= diag(1./scales)  * U(1:N,k(rangeW))';
  % W=固有値の平方根で規格化＋固有値基底への変換を行うホワイトニング行列
  X		= W*X;            % ホワイトニング実行

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%     学習        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 分離行列の初期条件
B=eye(N);  %単位行列に選ぶ
L=0.1;    %学習率 

% Bell & Sejnowsky のアルゴリズム
% tanh(u) の部分は他の非線形関数でもよい
% for t=1:P
%  u=B*X(:,t); 
%  B=B+L*(Id -tanh(u)*u')*B;
% end;

% EPS で収束条件を設定 shiro
EPS =0.005;
DELTA=eye(N);
while norm(DELTA,'fro')>EPS,
  u   =B*X; 
  tanu=tanh(u); 
  DELTA=L*(-(tanu*u'-u*tanu')/P)*B;
  B=B+DELTA;
end;

% ホワイトニングを含む行列に変換
B	= B*W ;
return;