本文書の著作権は株式会社シグナルアナリシスに帰属します．

Second-order blind source identiﬁcation(SOBI) アルゴリズム
株式会社1

シグナルアナリシス

関原謙介

データモデル
観測データがある時間間隔で次々と観測され，時系列データとして得られる状況を考える．したがって，観測データ y は時間変動するので y(t)：y1 (t) y2 (t) ,y(t) = ..
.
yM (t)(1)と表す．y(t) は離散時間点 t = t1 , . . . , tK で観測されるとする．K は観測時間点の総数である．この観測データがy(t) = Au(t),(2)と表されるとする．式 (2) において，A は混合行列（mixing matrix）とよばれる．当然ながら，

u(t) = A−1 y(t) = W y(t)も成り立ち，W = A−1 を分解行列（demixing matrix）と呼ぶ．ここで，u(t) は M 個の「独立成分」（independent component）からなる列ベクトルで，
u1 (t)..
u(t) = .

,(3)uM (t)である．各成分 u1 (t), . . . , uM (t) は統計的に独立であるとする．すなわち SOBI アルゴリズムは２次の統計量を用いて，観測データのタイムコース y(t)（t = t1 , . . . , tK ）を独立な変動 u(t)（t = t1 , . . . , tK ）に分解するアルゴリズムである [1, 2]．つまり，観測データ y(t)（t = t1 , . . . , tK ）から，u(t) の各成分が統計的に独立であるとの仮定をおき，独立成分 u(t)（t = t1 , . . . , tK ）と混合行列 A を推定する．

1

2SOBI アルゴリズム
要素 ui (t) と uj (t) のタイムラグ τ での相関は

1 X
ui (tk )uj (tk + τ ) = hui (t)uj (t + τ )iKK

φi,j (τ ) =(4)k=1で定義される．ここで，h·i は時間平均を意味する．タイムコース ui (t1 ), . . . , ui (tK ) と uj (t1 ), . . . , uj (tK )が独立であるとはφi,j (τ ) = 0,for all τ(5)が成り立つことである．SOBI アルゴリズムは，この条件をできるだけ満たす uj (t) を与える分解行列

W （= A−1 ）を見出す．
まず，ベクトル u(t) のタイムラグ τ での相関行列を Ru (τ ) と定義する．すなわち，

φ1,1 (τ )­®  φ2,1 (τ )T
Ru (τ ) = u(t)u (t + τ ) = ..
.

φ1,2 (τ )...

φ2,2 (τ ) . . .
..
..
.
.
φM,1 (τ ) φM,2 (τ ) . . .

φ1,M (τ )
φ2,M (τ ) 
..

.

(6)

φM,M (τ )である．i 6= j の場合において，タイムコース ui (t1 ), . . . , ui (tK ) と uj (t1 ), . . . , uj (tK ) が独立であれば全ての τ において Ru (τ ) が対角行列になる．
データベクトル y(t) に対するタイムラグ τ での相関行列を Ry (τ )：­®
Ry (τ ) = y(t)y T (t + τ )(7)

Ru (τ ) = W Ry (τ )W T(8)と定義する．ここで，の関係がある．SOBI アルゴリズムは複数個の Ry (τ ) を対角化する行列を見つけることで，式 (5) を満たす要素 uj (t) を与える分解行列 W を見出すアルゴリズムである．
まず，τ = 0 である相関行列 Ry (0) の固有値展開を

Ry (0) = V ΛV T(9)と表す．Λ は対角成分が固有値である対角行列であり，V は固有ベクトルを列ベクトルとする直交行列である．もし，hui (t)uj (t)i = 0（i 6= j ）を満たす uj (t) でデータベクトル y(t) を分解する分解行列 Wを求めるなら，W = V T とすればよい．実際，u(t) = W y(t) であるので，­®
­®
Ru (0) = u(t)uT (t) = V T y(t)y T V = V T Ry (0)V = V T V ΛV T V = Λ(10)となる．すなわち，得られた uj (t) は，hui (t)uj (t)i = 0（i 6= j ）を満たす．この場合には ui (t) と uj (t)は無相関ではあっても独立ではない．
少しでも独立性の条件に近づけるため，条件を１つ増やし，ある τ の値 1 個を用いて

hui (t)uj (t)i = 0およびhui (t)uj (t + τ )i = 02for i 6= j(11)を満たす uj (t) を与える分解行列 W を求めよう．この場合には，まず，次の行列 P を

P = Ry (0)−1/2 = V Λ−1/2 V T(12)と定義し，これを用いてデータベクトル y(t) をz(t) = P y(t)(13)と変換する．この操作は白色化 (prewhitening) と呼ばれる．白色化されたデータ z(t) に関するラグ τ の相関行列を Rz (τ ) として，Rz (τ ) の固有値展開をRz (τ ) = P Ry (τ )P = U Λz U T(14)と表記すれば，Λz は固有値を対角成分とする対角行列であり，U は直交行列である．そしてこれらを用いて分解行列は W = U T P と求まる．実際，­®
­®
Ru (0) = u(t)uT (t) = U T P y(t)y T (t) P U = U T P Ry (0)P U = U T IU = I(15)であるので，i 6= j に対して，まず，hui (t)uj (t)i = 0 は成立する．次に，­®
­®
Ru (τ ) = u(t)uT (t + τ ) = U T P y(t)y T (t + τ ) P U= U T P Ry (τ )P U = U T Rz (τ )U = U T U Λz U T U = Λz(16)が成り立つため，i 6= j に対して，hui (t)uj (t + τ )i = 0 も成立する．すなわち，ここで求まった uj (t) は単に無相関（hui (t)uj (t)i = 0）であるのみでなく，ある１つの τ において，hui (t)uj (t + τ )i = 0 を満たし，独立な信号に近づいたものとなっている．
さらに独立性の条件に近づけるため，複数個のタイムラグを用いて，hui (t)uj (t)i = 0および

hui (t)uj (t + τ` )i = 0for i 6= j` = 1, . . . , L(17)を満たす uj (t) を与える分解行列 W を求める1 ．τ が 1 個の場合と同様に，まず，白色化を行う．行列P を式 (12) で定義し，これを用いてデータベクトル y(t) をz(t) = P y(t)(18)と変換する．z(t) に関するラグ τ1 から τL までの相関行列を Rz (τ` )（` = 1, . . . , L）として，Rz (τ` ) を対角化する共通の直交行列 U を求める．このような U は存在し，文献 [4] に述べられている方法 (jointdiagonalization と呼ばれる) により求めることができる．このような U を用いて，TRz (τ1 ) = P Ry (τ1 )P = U Λ(1)
z U..
.
T
Rz (τL )= P Ry (τL )P = U Λ(L)z U(19)(20)(21)

1 τ が 2 個以上の SOBI アルゴリズムは TDZEP(temporal decorrelation source separation) と呼ばれる場合がある [3]．

3

と τ1 から τL までのラグを持つ相関行列を対角化できる．ここで，Λ(`)
z は対角行列であり，U は直交行列である．したがって，分解行列はやはり W = U T P と求まる．実際，式 (15) と全く同様に i 6= j に対して，hui (t)uj (t)i = 0 は成立する．次に，­®
­®
Ru (τ` ) = u(t)uT (t + τ` ) = U T P y(t)y T (t + τ` ) P UT(`)= U T P Ry (τ` )P U = U T Rz (τ` )U = U T U Λ(`)
z U U = Λz(22)が成り立つため，i 6= j に対して，hui (t)uj (t + τ` )i = 0 が ` = 1, . . . , L に対して成立する．したがって，ここで求まった成分 uj (t) は独立性の性質を強く満たすものである．

3

まとめSOBI アルゴリズムを用いたアーチファクトの除去は次の手順で行う．
1. データ y(t) に対するラグゼロの相関行列 Ry (0) を求め，この行列の固有値展開から白色化を行う
行列 P を求める．

2. 複数個のラグにおいて相関行列 R(τ` )（` = 1, . . . , L）を求め，これを共通対角化 ( joint diaginalization) することにより対角化直交行列 U を求める．
3. 分解行列（demixing）行列 W を W = U T P から求め，さらに，独立成分 u(t) を u(t) = W y(t)から求める．
4. タイムコース uj (t1 ), . . . , uj (tK )（j = 1, . . . , M ）を調べ，アーチファクトと思われる成分はゼロと置き換える．
e (t) とすれば，アーチファクト除去後の5. アーチファクト成分をゼロと置き換えた成分ベクトルを u−1e (t) を ye (t) = P U ue (t) から再構成する．
データベクトル y

参考文献[1] L. Molgedey and H. G. Schuster, “Separation of a mixture of independent signals using timedelayed correlations,” Phys. Rev. Let., vol. 45, pp. 434–44, 1997.
[2] A. Belouchrani, K. A.-Meraim, J.-F. Cardoso, and E. Moulines, “A blind source separation technique using second-order statistics,” IEEE Trans. Signal Process., vol. 45, pp. 434–44, 1997.
[3] A. Ziehe and K.-R. Müller, “Tdsep: an eﬃcient algorithm for blind separation using time structure,” in ICANN 98, pp. 675–680, Springer, 1998.
[4] J.-F. Cardoso and A. Souloumiac, “Jacobi angles for simultaneous diagonalization,” SIAM J.
Matrix Anal. Appl., vol. 17, pp. 161–164, 1996.

4東京医科歯科大学におけるMCGデータから冷凍機の振動ノイズ除去実験

MCGデータSVDによりdimension reduction:
120センサーデータ
20タームコースデータ
ICA アルゴリズムにより得られた２０個のコンポーネント（青いウィンドウで囲った成分を振動のイズ成分と推定）

特異値スペクトル

実験結果
波形拡大表示

処理前
処理前

処理後

処理後