信号処理の推定問題

講演内容

独立成分解析の基礎と応用

Ⅰ．ブラインド信号処理の枠組み■ 未知システムの同定 or 適応フィルタの設計

・ ブラインド信号処理取り扱う事例
・ モデルの設定及び定式化

xⅡ．独立成分解析
ソウ・ 独立成分解析とは何か
・ Kullback-Leibler情報量による独立成分解析
・ 独立成分解析の実装

ケンテイ

曹

建庭 （Jianting Cao）
cao@sit.ac.jp
埼玉工業大学 工学部

システム

x

d

+

e

-y

適応フィルタⅢ．独立成分解析による知的信号処理事例
・ 脳磁界(MEG)の信号処理
・ 脳波(EEG)によるBCI (Brain-Computer Interface)・ 脳波(EEG)によるBDD (Brain-Death Determination)

ブラインド信号処理の一推定問題s

未知システム

推定システムICAによるブラインド信号処理
y

信号源とシステムとも未知s

混合システムx

人工的な神経回路■ ブラインドシステム同定(blind system identification)の問題：● Blind Equalization● Blind Beamforming● Blind Array System● Blind Deconvolutiony■ Y. Sato (1975) [典型例]Two extensional applications of the zero-forcing equalization method,
IEEE Trans. on communication, Vol. com-23, pp. 679-682.
★ ブラインド同定とは原信号とシステムが未知と言う条件のもとで，システム或いは原信号を推定すること．推定法は様々なものがある．

ICAによるブラインド信号処理例■ 原信号の統計的な独立性を着目し，原信号の分離・復元，或いは混合システム推定を行う手法は独立成分解析(Independent
Component Analysis: ICA)と言う．
■ ICAの出典： P. Comon (1994)``Independent component analysis, a new concept?,'' Signal
Processing, vol. 36, no. 3, pp.287-314, 1994.
■ ICAでは，ブラインド信号分離(Blind Source Separation)の問題が始め，●Blind Equalization，●Blind Beamforming, ● Blind
Array System, ● Blind Deconvolutionにも多く適用されている．

ICAのデモ例s

未知システム

x

適応等化器

y

ICAに基づくBSS研究の流れ

信号源とシステムとも未知■ ブラインド(blind)とは信号源とシステム両方が未知を指している．

■ 逆システムの同定 or 逆フィルタ（等化器）の設計◎C. Jutten and J. Herault, (1991) ``Blind separation of sources, part I : An adaptivealgorithm based on neuromimetic architecture,'' Signal Processing, vol. 24, no. 1, pp.1-10, 1.
◎ P. Comon (1994) ``Independent component analysis, a new concept?,'' Signal Processing,vol. 36, no. 3, pp.287-314.
◎ J. F. Cardoso and A. Souloumiac (1996) : Jacobi angles for simultaneous diagonalization,
SIAM J. Mat. Anal. Appl., vol. 17, no. 1, pp.145-151.
◎ J. F.~Cardoso and B.~Laheld, (1996) ``Equivariant adaptive source separation,'' IEEE
Trans. on Signal Processing, vol. 44, no. 12, pp.3017-3030.
◎ A.T. Bell and T.J. Sejnowski, (1995) ``An information maximization approach to blindseparation and blind deconvolution,'‘ Neural Computation}, vol. 7, no. 6, pp.1004-1034.
◎ S. Amari, A. Cichocki and H.H. Yang, (1996) ``A new learning algorithm for blind signal
separation,'' Advances in Neural Information Processing System 8, MIT Press, pp.757-763.
◎ A. Hyvarinen and E. Oja, (1997) ``A fast fixed-point algorithm for independent component
analysis,'' Neural Computation, vol. 9, pp.1483-1492.

ICAによるノイズ画像の復元■ 聖徳太子型の音声の聞き分け
複数の会話が混じり合う騒音のなか，ある聴者が特定の話音に注意力を払えば，その音声だけを抜き出せることできる．このような聞き分け効果はカクテルパーティー効果と言う．本研究では，人間が持っている聞き分け機能が果たして機械で実現できるのかと言う問題を対象にする．

互いに独立

I
C
A

マイクで収録された混合音
声から，元個々の音声に復元すること．
数理的な手法だけとしての実現には：【独立成分解析 ICA】
情報理論，統計学
システム推定理論
ニューラルネットワーク
FIR or IIR フィルタ

互いに独立

Unknown

1

ICAによるブラインド等化システム■

動体通信における等化システ
ム
アンテナの信号から混信信号を分解し，元の個々の送信信号へ復元すること．

混合モデルの設定

ICAによる脳信号処理例■ EEG/MEG脳信号処理
電極から計測した脳信号から，脳がどこに、どのように活動しているかをとくていすること．

Noise sources :MEG signals

Power supply
interference

EF
BG

noise

Artifacts sources■ 混合モデルの設定については，取り扱っている実際の応用
問題に応じて行うべきである．次の3種類の方式で原信号が混合されることが多い．
▼ 空間混合モデル(spatial mixture)：● 線形行列● 脳活動から計測電極への転送路▼ 時間的混合モデル(temporal mixture)：● FIRフィルタ● 等化システムの転送路▼ 時空間的混合モデル(spatio-temporal mixture)：● FIRフィルタ行列● 会話音声からマイクへの音響転送路

時間的混合モデルの定式化■ 時間的混合モデルの表現

時空間的混合モデルの定式化■ 時空間的混合モデルの表現

K −1

x(t ) = ∑ ak s (t − k ), t = 1,2,!
k =0

離散時間： tインパルス応答： a0, a1,…, aK-1

フィルタ次数： K

観測信号時間列： x(t), x(t-1),…, x(t-K+1)■ 時間的混合モデルの表現は基本的にFIRフィルタを用いられている．
観測信号には原信号とインパルス応答の畳み込みで得られている．
このモデルはブランイド等化によく使われている．

観測信号の分解による解の表現■ 次の瞬時混合モデル

x(t) = As(t)に対して，各成分が独立である原信号と対応づける観測信号x(t)の分解

x(t)=a(1)(s (t))+!+a(n)(sn(t))
1を求める．

■ 多変量データ解析の手法と近い考え方．

x(t) = As(t)K −1x(t ) = ∑ Ak s(t − k ), t = 1,2 !
k =0

原信号時間列： s(t), s(t-1),…, s(t-K+1)

空間的混合モデルの定式化■ 空間的混合モデルの表現

原信号ベクトル： s(t-k)
行列インパルス応答： Ak
観測信号ベクトル： x(t)■ 時空間的混合モデルは複数個の信号源が現時点の値だけでなく，過去の値も加えて複数の行列により混合されている．多入力多出力(MIMO)コンボリューションモデルとも呼ばれている．このモデル
はカクテルパーティー効果を模倣した音声分離やMIMOブランイド
等化などによく使われている．

ニューラルネットワークによる解の表現■ 次の瞬時混合モデル

x(t) = As(t)に対して，ニューラルネットワークの重みWを用いて観測信号x(t)から復元信号y(t)をy(t) =Wx(t)によって構成したとき，復元信号の各成分が互いに独立になるように重みWを求める．
■ この手法を用いる場合には，復元信号の順序や大きさは必ずしも原信号と一致しない．

原信号： s = [ s1 ,!, sn ]T
混合行列： A = (aij ) i = 1,!, n, j = 1,!, mT

観測信号： x = [ x1 ,! , xm ]■ 空間的混合モデルは，時空間的混合モデルが時間遅れで生じた間接波を無視した場合のモデルと考えられる．このモデルは原信号がセンサーへ瞬間的に伝搬するため，瞬間混合モデルとも呼ばれている．通常脳波(EEG)や，脳磁気(MEG)などの生体信号が電極への伝搬速度が速く，瞬間混合モデルがよく適用されている．また，このモデルは非常に単純であるため，多くのアルゴリズムはそれに準じて開発されている．

Ⅱ．独立成分解析
１．独立成分解析とは何か
・ 独立性の定義
・ ブラインド信号処理への適用
・ 独立性の基準

２．Kullback-Leibler情報量による独立成分解析
・ 確率密度関数の近似法
・ 一般化された正規分布のモデル
・ t-分布モデル
・ 尖度及び単峰形分布との関係
・ スコア関数の決定
・ アルゴリズムの安定性・ 誤推定された場合の安定性

３．他の独立成分解析アルゴリズム
４．独立成分解析の実装
・ MATLABによる簡単な実装例
・ WEBに公開されているMATLABプログラム

2

原理：独立成分解析 (ICA)
n個互いに独立の信号源

s1 (t )s2 (t )

x(t) = As(t)x1 (t )x2 (t)

y(t) =Wx(t)y1 (t )y2 (t )

例題：原信号の分離・復元!
Asn (t )!

xm (t )!

yn (t )

続き：復元側の表現■ ニューラルネットワークの表現y (t ) = Wx(t )W

x(t ) = [ x1 (t ) x2 (t )]Ty(t ) = [ y1 (t ) y 2 (t )]T■ 復元信号と原信号の関係：■ 瞬時混合モデルの表現

W = A−1 ならば、 y=s

分離過程：
ｙ～ｓになるようにうまくWを見つける。

独立成分解析 (ICA)とは観測信号を利用するだけで、元の信号源にある独立性を再構成する方法

続き：独立性に基づく可能な解⎡w11 w12 ⎤W =⎢⎥⎣w21 w22 ⎦x1 (t ) = a11s1 (t ) + a12 s2 (t )y (t ) = WAs(t )t = 1, 2,!

x2 (t ) = a21s1 (t ) + a22 s2 (t )x(t ) = As(t )s(t ) = [ s1 (t ) s 2 (t )]Tx(t ) = [ x1 (t ) x2 (t )]T

⎡a11
A=⎢
⎣a 21

a12 ⎤
a 22 ⎥⎦

まとめ

重み行列は自由に調整・更新が可能

独立性の定義■ 復元信号は原信号と等しい：
W =A

−1

ならば y (t ) = s(t )■ 行列係数Aが既知の場合には解を求めることは簡単である．

原信号は観測信号から忠実に分離される■ 現実では、係数Aが未知であり、また重み行列Wを構成する係数の可能な組み合わせは無数にあるため，上述の理想的な解を求める
ことはほとんど不可能である．

■ 復元信号は原信号の大きさと順番が異なる場合 ：W = PDA −1
⎡1 0⎤
P=⎢
⎥ or
⎣0 1⎦

⎡0 1⎤
P=⎢
⎥
⎣1 0⎦

y1 (t ) = d1 s1 (t )y2 (t ) = d 2 s2 (t )ornp y (y ) = ∏ pi ( y i )i =1■ 自動的に重みWの調整・更新則（計算アルゴリズム）の設計が必要．

⎡d1 0 ⎤
⎡0 d1 ⎤
D=⎢⎥ or D = ⎢d 0⎥⎣0 d 2 ⎦⎣ 2⎦⎡ y1 (t ) ⎤ ⎡1 0⎤ ⎡d1 0 ⎤ ⎡ s1 (t ) ⎤
⎢
⎥=⎢
⎥⎢
⎥
⎥⎢
⎣ y 2 (t )⎦ ⎣0 1⎦ ⎣0 d 2 ⎦ ⎣ s 2 (t )⎦■ n個の信号 y = [y1, y2, …, yn]T（それぞれN個のサンプルが持つ）
を確率変数とし，同時(結合)分布の密度関数をpｙ(y)，周辺分布の密度関数をpi(yi)とする場合には，同時(結合)分布の密度関数が周辺分布の密度関数の積と等しくなると，独立と定義される．

" y1 (t) % "1 0 % "0 d1 % "s1 (t) %
$
'=$
'$
'$
'
# y 2 (t) & #0 1& #d2 0 & #s2 (t) &■ どのように計算アルゴリズムを設計するか？
どんな基準のもとで設計するか？
独立性って何か？

■ 原信号が互いに独立と仮定されている．原信号が混合行列によって混合され，独立性が崩れたまま伝達され，観測信号になる．そこで，観測信号を利用した復元信号の独立性を再構成すれば，原信号がもとの独立性に戻ると考えられる．
■ 独立性を再構成するために，ある目的関数での評価が必要．

€or

y1 (t) = d1s2 (t)€€y 2 (t) = d2 s1 (t)独立性の基準■ Kullback-Leibler情報量：分布の差を直接測る方法■ 最尤法：
重みの最尤推定量を求める方法■ 相互情報量の最小化：
2変数の場合には， Kullback-Leibler情報量と同じ．
■ エントロビー最小化：
正規分布から離れた度合いの評価■ 情報量の最大化：
伝送路に雑音がある場合できるだけ情報量を大きくする方法．
★ 上記の基準がそれぞれ異なるが，導いたアルゴリズムはほぼ同じ．

その他の独立性評価基準

Kullback-Leibler情報量■ 高次統計量による独立性の評価：
1．奇数次のクロスモーメントによる評価
2．クロスキュムラントによる評価■ Kullback-Leibler情報量とはある確率変数の同時(結合)分布の密度関数と周辺分布の密度関数の積との統計的な距離を測る
一尺度である．

■ 非線形相関による独立性の評価■ n個の復元信号 y = [y1, y2, …, yn]T（それぞれN個のサンプルが持つ）を確率変数とする場合のK-L情報量：

D(y W) = ∫ p y (y)log

p y (y )n∏ p (y )
i =1

i

dy

i

■ pｙ(y) ： 同時(結合)分布の密度関数
pi(yi) ： 周辺分布の密度関数

村田 昇： 独立成分解析 東京電機大学出版社，2004年

3独立条件とK-L情報量

エントロピーとK-L情報量

評価関数■ x=W-1yを利用したエントロピーH(y|W)■ 復元信号yのエントロピー■ 独立条件np y ( y ) = ∏ pi ( yi )D(y W) = ∫ p y (y )log

p y (y )n∏i=1 pi ( yi )■ Kullback-Leibler情報量の可計算な評価関数■ 各成分yiのエントロピー

dy = 0n■ エントロピーを用いてKullback-Leibler情報量の表現■ 偏微分により得られた逐次学習アルゴリズム：

ΔW(t ) = η[ W −T (t ) − ϕ (y (t ))xT (t )]■ 自然勾配法＊∂D(y | W ) T
W W =0∂W

D ( y W ) = − H ( y | W ) + ∑ H ( yi | W )スコア関数につて分散が一定である各成分yi (i=1,…,n)の分布のなかに，正規分布を持つ成分はそのエントロピーが最大であり，正規分布に離れている分布を持つ成分のほうがそのエントロピーが逓減される．また，合計エントロピーが最小になれば，各成分が互いに独立になると言える．

補足：高次モーメントについて■ スコア関数の各成分 ϕi (⋅) は信号の確率密度関数 pi (⋅) に依存している．
通常 pi (⋅) が未知であるので，何らかの方法で関数の近似が必要である．
■ 直感的なスコア関数の選択方法
１．裾が重い分布である音声信号に対して次の関数を適用する．

ΔW(t ) = η[I − ϕ (y(t ))yT (t )]W(t )ϕ i ( yi ) = −■ 右端第一項のみは行列Wに制約されているため， K-L情報量の最小化は，右端第一項の最小化と一致する．

■ 分布の形状を定める一つの規準としての尖度(kurtosis)は2次と4次のにより導かれたアルゴリズム：■ 最適なスコア関数 ：i =1★ エントロピー関数の性質により，エントロピーが最大となるのは不確定度が最大である．

n

i =1

学習アルゴリズムD(y W) = ∑ H( yi | W) −H(x) − log | W |H( yi W) = −∫ pi ( yi )logpi ( yi )dyi

■ Kullback-Leibler情報量は重み行列Wの関数とする場合，変数yの成分間の独立性の距離を最小化にする一つの評価基準であることが考えられる．

∂D(y | W )=0∂W

H(y W) = H(x) + log | W |H(y W) = −∫ p y (y)logp y (y)dyi =1が満足すれば，Kullback-Leibler情報量は

ϕ ( y) = tanh( y)２．裾が軽い分布である画像・通信信号に対して次の関数を適用する．

! (y )pdlog pi ( yi ) = − i idyipi ( yi )ϕ ( y) = y 3

m4 ( y )中心モーメントを用いて k ( y) = m ( y) − 3 で定義される．
2
2■ kが正ならば，信号の分布（裾が重い）はSuper-Gaussianと呼ばれる分布になる．音声
信号はSuper-Gaussianである
ことが多い．
■ kが負ならば，信号の分布（裾が軽い）はSub-Gaussianと呼ばれ，画像信号などはこの性質を持つことが多い．
■ 正規分布では k=0 となる．

＊

確率密度関数の近似法

一般化された正規分布のモデル■ モデルの定義：■ 直感的なスコア関数の選択方法(non- parametric法)に対して，
観測信号から計算された統計量(尖度：kurtosis)を用いて母数を決定するparametricな方法がある．

■ 一族分布の形を制御する係数を利用して，統計量の尖度(kurtosis)と結び付こうとしている．

pα ( y ) =

αλα
2Γ(

α

分散を正規化する係数
α ： 分布の形を制御する係数α)exp(− λα y )∞Γ( x) = ∫ t x −1 exp(−t )dt

一般化された正規分布と尖度の関係■ 一般化された正規分布とその高次モーメント0pα ( y ) =λα ：■ 2次中心モーメント∞■ 係数α変化により，単峰形の一族の分布が得られる．
α = 2 は正規分布である．

■ 単峰の形にしている一族分布のモデルを対象する．
１．一般化された正規分布のモデル
２．t-分布をモデル

1

α < 2 は一族の裾が重い分布
α > 2 は一族の時裾が軽い分布

m2 = ∫ y 2 pα ( y )dy =
−∞

αλα
2Γ(

α

)

exp(− λα y )■ 4次中心モーメント

3
Γ( )

α

1

α

∞

m42 = ∫ y 42 pα ( y )dy =1−∞

λα2 Γ( )
α

53
Γ( )

α

1

λα42 Γ( )α■ 正規化係数13 ⎤2⎡⎢ Γ( α ) ⎥1 ⎥⎢ m2 Γ( ) ⎥⎣⎦λα = ⎢

α■ 尖度と分布の形を制御する係数αの関係式
kα =

Γ(

5

α

)Γ (

Γ2 (

3

α

1

α

)
−3

)

★ 上記の結果により，信号の尖度kαを既知とした場合には係数αを求めることができる．

4スコア関数アルゴリズムの安定条件■ 最適なスコア関数
ϕ i ( yi ) = −■ Kullback-Leibler情報量に基づく学習アルゴリズムに対して，Hessianの正定性の評価による得られた安定性の必要十分条件[1]p! ( y )dlog pi ( yi ) = − i idyipi ( yi )安定性条件の導き■ 安定性条件の要素の計算

E[ϕ!i ( yi )] > 0■ 一般化された正規分布∞αλαpα ( y ) =2Γ(E[ yi2 ] = ∫ y 2 pα ( y )dy =

α

1

α

)

exp(− λα y )

−∞

2
i i

E[ y ϕ! ( yi )] + 1 > 0

αλα

3 ⎤
⎡
⎢ Γ( α ) ⎥1 ⎥⎢ m2 Γ( ) ⎥α ⎦⎣E[ϕ!i ( yi )] = ∫ ϕ!i ( yi ) pα ( yi )dyi =∞[ ]E yi2 E [ϕ!i ( yi )] =E[ϕ!i ( yi )] =λα2 α (α − 1)Γ(Γ(

E[ yi2 ]E[ϕ!i ( yi )] − 1 =

1

α

α −1
)
α
>0

)

3
α −1
α (α − 1)Γ( )Γ()αα −1 > 0
1
Γ2 ( )■Case 3 : α < 2super-Gaussianの信号に対して，条件①-③は常に満たされていることではない．例えば， α =0.2, 0.5,のとき，条件①-③にあるGamma関数が不連続である．

t-分布によるスコア関数■ 最適なスコア関数

2

)
(1 +

(λβ y ) 2 − 12 ( β +1))

ϕ i ( yi ) = −β2■ 2次中心モーメント∞m2 = ∫ y 2 pβ ( y )dy =
−∞

β Γ(

β −2
2
2

β

■ 4次中心モーメント)∞, β > 2,m4 = ∫ y 4 pβ ( y )dy =2λβ Γ( )2

3β 2 Γ (

−∞

4

β −4
2

β

4λ β Γ ( )2■ 正規化係数1), β > 4,β − 2 ⎤2⎡⎢ β Γ( 2 ) ⎥β ⎥⎢ 2m2 Γ( ) ⎥2 ⎦⎣p! ( y )dlog pi ( yi ) = − i idyipi ( yi )■ t-分布pβ ( y ) =λβ Γ (

m4
−3 =
m22

3Γ(

β −4

β

)Γ ( )
2 −3
β −2
Γ2 (
)
2

2

α

0

3

α

Γ2 (

)Γ (

1

α

ϕi ( y ) =pβ ( y ) =λβ Γ (

β +1
2

)

β
πβ Γ( )
2

(1 +

α −1
)
α

)

(λβ y ) 2 − 12 ( β +1)
)

β

α +1
λα +1α 2 (α − 1) Γ( α )× α +1= α −1
1
λα α
Γ( )
α

∞

Γ( x) = ∫ t x −1 exp(−t )dt
0

λβ ：
β

分散を正規化する係数： 分布の形を制御する係数■ t-分布族(0 < beta < ∞)は裾が重い形をしており，betaが無限大につれ正規分布に近づく．
■ 特別の場合として，beta=1のときの分布はcauchy分布と呼ばれる．
■ t-分布族は時裾が軽い分布がな
いこと．

アルゴリズムの安定条件■ Kullback-Leibler情報量に基づく学習アルゴリズムに対して，Hessianの正定性の評価による得られた安定性の必要十分条件[1]E[ϕ!i ( yi )] > 0)(1 +■ t-分布のモデルによるスコア関数2★ 上記の結果により，信号の尖度を既知とした場合には係数betaを求める
ことができる．

β +1βπβ Γ( )2λβ = ⎢■ 尖度と分布の形を制御する係数betaの関係式
kβ =α (α − 1)Γ(

α −1
α −1
2
λαα +1α 2 (α − 1) Γ( α ) λα α (α − 1)Γ( α )× α −1
=
1
1
λα α
Γ( )
Γ( )
α
α

∫ y exp(−(λα y) )dy = α1Γ( )t-分布モデル■ t-分布モデルβπβ Γ( )λ α (α − 1) ∞ α■ 自由度betaを持つt-分布の確率密度関数：t-分布と尖度の関係
λβ Γ (exp(−(λα y )α )dy =安定性の検証★ 安定性を満たすための他のモデルがあるか？

pβ ( y ) =

α −2

0

2

■Case 1 : α > 2
sub-Gaussianの信号に対して，条件①-③は常に満たされている．
■Case 2 : α ＝ 2Gaussian信号に対して，条件①-②が常に満たされている．条件③ においては，i番目の信号に対して③i=1であるので，仮にもう一つの信号（j番目）も正規分布であれば，③j=1となるので，結局条件③＝１となり，安定条件が満たされていない．このことは，2個以上のGaussian分布の信号を分離することが出来ないことを示している． しかし，もう一つの信号がsub-Gaussian，或いはsuper-Gaussianの信号であれば，条件③ （>1）が満たされている．

α

β +1

α +1
α

3
Γ( )

α =
α
1
αλα3
λα2 Γ( )ααS. Amari, T. Chen and A. Cichocki,``Stability analysis of adaptive blind source separation,'' Neural
Networks, vol. 10, no. 8, pp.1345-1351, 1997.

▼ 条件②∫ (λα y)

1
Γ( )

−∞

■ 一般化された正規分布のモデルを利用する時の安定性条件：

3
Γ( )
×

α

λα = ⎢

E[ yi2ϕ!i ( yi )] + 1 = α > 0

λα3 α 2 (α − 1) ∞E[ yi2ϕ!i ( yi )] = ∫ yi2ϕ!i ( yi ) pα ( yi )dyi =導いた安定性の条件

αλα

1
Γ( )

α

−∞

★ 上記の安定性条件を利用して，一般化された正規分布のモデルの安定性の検証ができる．

1
2

▼ 条件①

y 2 exp( −(λα y )α )dy =αE[ yi2 ]E[ y 2j ]E[ϕ!i ( yi )]E[ϕ! j ( y j )] − 1 > 0

ϕi ( y) = αλα sgn( yi ) | λα yi |α −1 ,

∞

1 ∫0
Γ( )

∞

■ 一般化された正規分布のモデルによるスコア関数▼ 条件③

ϕ!i ( yi ) = α (α − 1)λα2 | λα yi |α −2■ スコア関数の微分(1 + β ) yi
yi2 +

β
λ2β
1

β − 2 ⎤2⎡⎢ β Γ( 2 ) ⎥β ⎥⎢ 2m2 Γ( ) ⎥2 ⎦⎣(λβ y ) 2 − 12 ( β +1))

βE[ yi2ϕ!i ( yi )] + 1 > 0

E[ yi2 ]E[ y 2j ]E[ϕ!i ( yi )]E[ϕ! j ( y j )] − 1 > 0,★ 上記の安定性条件を利用して，一般化された正規分布のモデルの安定性の検証ができる．

λβ = ⎢

S. Amari, T. Chen and A. Cichocki,``Stability analysis of adaptive blind source separation,'' Neural
Networks, vol. 10, no. 8, pp.1345-1351, 1997.

5t-分布の安定性条件の導き1t-分布の安定性条件の導き２■ t-分布によるスコア関数の微分■ 安定性条件の要素の計算

ϕ!i ( yi ) =(1 + β )(c − yi2 )( yi2 + c) 2

■ t-分布のモデルを利用する時の安定性条件：
1∞

E[ϕ!i ( yi )] = ∫ ϕ!i ( yi ) pβ ( yi )dyi = 2cβ c 2( β +1)−∞1= 2cβ c 2

■ 安定性条件の要素の計算( β +1)

∞
⎡ ∞
⎤
yi2
c
( β + 1) ⎢ ∫dyi − ∫
dy i ⎥
2
0.5 ( β + 5 )0 ( y 2 + c ) 0.5 ( β + 5 )i⎣ 0 ( yi + c )
⎦

⎡ 1
β +4
3
β +2 ⎤2) Γ ( )Γ ()⎥⎢ Γ ( )Γ (2
2 ⎥ = λβ ( β + 1)
( β + 1) ⎢ β2
− β2
β +3
⎢ 2c 2 +1Γ( β + 5 ) 2c 2 +1Γ( β + 5 ) ⎥
2
2 ⎦⎥
⎣⎢

∞

∞

∞

−∞

−∞

E[ yi2 ] = ∫ y 2 pβ ( y )dy = ∫ y 2cβ (1 +

1
q −1
⎤
( β +1) ⎡ 1
y 2 − 12 ( β +1)(−1) r +1y)dy = 2cβ c 2)⎥⎢ q −1 ∑q − 2 Cr −1 (β⎢⎣ c r =1 2r + 1
y 2 + c ⎥⎦ 0

1∞

E[ y 2ϕ!i ( yi )] = ∫ yi2ϕ!i ( yi ) pβ ( yi )dyi = 2cβ c 2
−∞

n

⎛n⎞

(−1) rnr =0r

Γ(n + 1)Γ(α )[ ]E yi2 =

⎝ ⎠

βΓ(
2

β −2
2

( β +1)λ2β ( β + 1)
>0
β +3

▼ 条件①

E[ϕ!i ( yi )] =▼ 条件②

E[ yi2ϕ!i ( yi )] + 1 =▼ 条件③

E[ yi2 ]E[ϕ!i ( yi )] − 1 =

∞
⎡ ∞
⎤
cyi2
yi4
( β + 1) ⎢ ∫dyi − ∫
dy i ⎥
2
0.5 ( β + 5 )0 ( y 2 + c ) 0.5 ( β + 5 )i⎣ 0 ( yi + c )⎦⎡ 3 β +25β ⎤) Γ ( )Γ ( ) ⎥1⎢ Γ ( )Γ (( β +1)2 −2
2 ⎥ = β −3= 2cβ c 2
( β + 1) ⎢ 2ββ
⎢ 2c 2 Γ( β + 5 ) 2c 2 Γ( β + 5 ) ⎥ β + 3
⎢⎣
2
2 ⎥⎦

∑ r + α C ⎜⎜ r ⎟⎟ = Γ(n + α + 1)

t-分布モデルの安定条件

β ( β + 1)Γ(

β −2
2

β

)
>0

2( β + 3)Γ( )
2

)

β

2β
>0
β +3

β ( β + 1)Γ( β 2−2 )E[ yi2 ]E[ϕ!i ( yi )] =β2( β + 3)Γ( )2λ β Γ ( )22t-分布モデルの安定性検証

t-分布＋裾が軽い一般化された正規分布のモデル

尖度と合成 (単峰)モデルの関係■ beta>2のsuper-Gaussianの信号に対して，条件①-③は常に満たされている．
■ betaが無限大につれ， 条件③ が1に近づく．この場合には信号が正規分布となる．
★ 安定性が満たされているモデルの構成
１．一般化された正規分布のモデル（ sub-Gaussian ）２．t-分布をモデル（ super-Gaussian ）正規分布族のモデル
sub-Gaussian

t-分布をモデル
super-Gaussian

合成モデル

p(α , β , y )J. Cao, N. Murata, S. Amari, A. Cichocki and T. Takeda : ``A robust approach to independent component analysis
with high-level noise measurements,'‘ IEEE Trans. on Neural Networks, Vol. 14, No. 3, pp.631-645, May 2003.

スコア関数

スコア関数の比較

α, βκ (α , β )ϕ (α , β , y )計算の流れ
ニューラルネットワークの結合重み W ,２,４次のモーメント m
ˆ , mˆ の初期値を設定する■ sub-Gaussian とGaussian信号に適用する関数2

分散を正規化する係数を
1

1

3
⎡
⎤2
β − 2 ⎤2⎡⎢ Γ(α ) ⎥⎢ βΓ ( 2 ) ⎥⎥⎥ λα = ⎢⎢ m Γ( 1 ) ⎥⎢ 2m 2 Γ ( β ) ⎥⎢ 2 α ⎦
⎥
⎣
⎢
⎥
2 ⎦
⎣

4

λβ = ⎢ϕi ( y) = αλα sgn( yi ) | λα yi | , kα = kˆi ≤ 0
α −1

⌢
⌢
j
Kurtosisの値を
m j (k ) = [1 − η ]m j (k − 1) + ηyi (k )⌢j = 2,4により計算する⌢mki = ⌢ 42 − 3 出力信号 y を y = Wx により計算する
m2■ super-Gaussian 信号に適用する関数
非線形関数を

ϕi ( y ) =(1 + β ) yiβ
yi2 + 2λβ

Kurtosis- α ,β, k β = kˆi > 0kα

ϕ i ( yi ) =■ 関数適用ための2次及び4次モーメントの推定ˆ j (t ) = [1 −η ]m
ˆ j (t − 1) + ηyij (t ), j = 2,4myi → ∞(1 + β ) yi
yi2 + β / λ2βϕ ( yi ) → 0

m4
=
−3 =
2
m2

kβ =

m4
m2

2

−3=

5

α

)Γ (

Γ2 (

3Γ (

ϕ i ( yi ) =

の関係式

Γ(

β −4
2

Γ2 (

3

α

1

α

2

−3

ϕi ( yi ) = αλα sgn( yi ) λα yi(1 + β ) yi
2

yi +
α −1

β
λ2β

⌢
ki > 0⌢ki ≤ 0 により計算する

)

)Γ (

β −2

)

β
2

)
−3

)

により,αとBetaを求める．
No

[

T

]

W を更新値を W(t + 1) = W (t ) + η(t ) I − ϕ (y(t ))y (t ) W(t )より求める

t = t last
Yes
出力信号y
算する

の最終更新値をy = Wxにより計6誤推定された場合の安定性■ 問題提起：１．次のモーメントが正しく推定されていない場合には安定条件が満たされるか？
２．満たされたとしても制御係数がある許容範囲に存在するか？

ˆ j (t ) = [1 −η ]m
ˆ j (t − 1) + ηyij (t ), j = 2,4m■ パラメータalpha，betaの推定値は真値からepsilonずれたとする場合

alphaの誤推定場合の安定条件■ 誤推定場合のスコア関数の微分

ϕi,αˆ ( yi ) = α + ελα sgn( yi ) | λα yi |α +ε −1■ alpha誤推定される場合の安定条件
E[ yi2ϕ!i ,αˆ ( yi )] + 1 =αˆ = α + ε((α + ε ) 2 − α − ε )Γ[(α + ε + 1) α ]+1 > 0Γ(1 α )

E[ϕ!i ,αˆ ( yi )] =λα2 (α 2 + αε )Γ (2α + ε − 1 α )>0Γ(1 α )αˆ = α + ε
α = 4βˆ = β + ε(α 2 + αε )Γ(3 α )Γ (2α + ε − 1 α )E[ y ]E[ϕ!i ,αˆ ( yi )] − 1 =−1 > 0Γ 2 (1 α )ε > −0.92ibetaの誤推定場合の安定条件

beta誤推定の許容範囲■ 誤推定場合のスコア関数の微分
ϕ i , βˆ ( yi ) =alpha誤推定の許容範囲

ˆ > 3.1
α

同類のICAアルゴリズム■ M. Girolami and C. Fyfe (1996)(1 + β + ε ) yi
β +ε
yi2 +
2

ΔW(t ) = η[I − sign (k ) tanh(y(t ))yT (t ) − y(t )yT (t )]W(t )λβｋは各成分ごとの尖度(kurtosis)■ beta誤推定される場合の安定条件

E[ϕ!i , βˆ ( yi )] =λ2β ( β + ε + 1)>0
β +3E[ yi2ϕ!i ,βˆ ( yi )] + 1 =( β + ε + 1)(β − 3)
+1 > 0β 2 + 4β + 3

E[ yi2 ]E[ϕ!i ,βˆ ( yi )] − 1 =β ( β + ε + 1)Γ( β 2 − 1)
−1 > 02( β + 3)Γ( β 2)

■ J. F.Cardoso and B.Laheld (1996)βˆ = β + ε
β = 4ε> −1.5

ˆ > 2 .5βΔW(t ) = −η[ϕ (y(t ))yT (t ) − y(t )ϕ (y(t ))T ]W(t )ϕ (⋅) はスコア関数
M. Girolami and C. Fyfe ``Negentropy and kurtosis as projection pursuit indices provide generalized ICA algorithms,''
NIPS'96 Workshop: Blind Signal Processing}, Snowmaas, Colorado, Dec, 1996.
J. F.Cardoso and B.Laheld,``Equivariant adaptive source separation,'' IEEE Trans. on Signal Processing}, vol. 44, no. 12,
pp.3017-3030, 1996.

ICAによる音声分離の実装例% Workspace変数をすべて削除
clear all% 音声信号の読み込み、信号源の生成
load birdchirps1;
load mosimosi1
s=[birdchirps1;mosimosi1];% 音声信号を試聴する
clc;input('鳥の鳴き声 … ')sound(birdchirps1);
clc;input('もしもし … ')sound(mosimosi1);% randn()を用いて混合係数を生成
A=randn(2,2);% 混合信号の生成
x=A*s;% 混合信号を試聴する
clc;input('混合信号① … ')sound(x(1,:));
clc;input('混合信号② … ')sound(x(2,:));% 混合信号の保存
save x x A% 信号波形の表示
subplot(2,2,1);plot(s(1,:));
xlabel('time');ylabel('s_{1}');title('(a) Birdchirps')
subplot(2,2,2);plot(s(2,:));
xlabel('time');ylabel('s_{2}');title('(b) Moshi-Moshi')
subplot(2,2,3);plot(x(1,:));
xlabel('time');ylabel('x_{1}');title('(c) Mixed signal ①')
subplot(2,2,4);plot(x(2,:));
xlabel('time');ylabel('x_{2}');title('(d) Mixed signal ②')% 学習回数、学習率、重みの初期値の設定
T=10000;eta=1.8;W=eye(2);% 逐次学習for t=1:T;home;y(:,t)=W*x(:,t);
W=W+eta/t*[eye(2)-tanh(y(:,t))*y(:,t)']*W;
W*A
end% 計算結果の表示
subplot(2,2,1);plot(x(1,:));
xlabel('time');ylabel('x_{1}');
title('(a) Mixed signal ①')
subplot(2,2,2);plot(x(2,:));
xlabel('time');ylabel('x_{2}');
title('(b) Mixed signal ②')
subplot(2,2,3);plot(y(1,:));
xlabel('time');ylabel('y_{1}');
title('(c) Separated signal ①')
subplot(2,2,4);plot(y(2,:));
xlabel('time');ylabel('y_{2}');
title('(d) Separated signal ②')% 分離した信号を試聴する
clc;input('分離した信号① … ')sound(y(1,:)/max(y(1,:)));
clc;input('分離した信号② … ')sound(y(2,:)/max(y(2,:)));

音声分離の実行結果

画像復元の実装例% Workspace変数をすべて削除
clear all
% 元画像データの読み込み
S1=double(imread('testpat1.tif'));
S2=double(imread('rice.tif'));% 元画像の表示
subplot(3,2,1);image(S1/2);title('(a) Original image ①')
subplot(3,2,2);image(S2/2);title('(b) Original image ②')
hold on% ２次元原画像データを１次元データに落とす[I J]=size(S1);s=[reshape(S1,1,I*J);reshape(S2,1,I*J)];% 混合係数と混合信号の生成
A=[-0.4689 0.3821; -0.3659 1.0933];
x=A*s;% １次元のデータを２次元配列データに戻す
X1=reshape(x(1,:),I,J);
X2=reshape(x(2,:),I,J);% 混合画像の表示
subplot(3,2,3);image(X1/2);
title('(c) Mixed image ①')
subplot(3,2,4);image(X2/2);
title('(d) Mixed image ②')
hold on% 学習データの作成
trainx=x-mean(x')'*ones(1,length(x));
trainx=trainx/max(max(x'));% 学習回数、学習率、重みの初期値の設定
T=10000;
eta=100;
W=eye(2);% 逐次学習for t=1:T;if rem(t,1000)==0;tend;home;
trainy(:,t)=W*trainx(:,t);
W=W+eta/T*[eye(2)-trainy(:,t).^3*trainy(:,t)']*W;
G=W*A
End
% 学習した復元データの計算
y=W*x;% ２次元の復元画像データへ変換
Y1=reshape(y(1,:),I,J);
Y2=reshape(y(2,:),I,J);% 復元画像の表示
scale1=2*G(1,1);scale2=2*G(2,2);
subplot(3,2,5);image(Y1/scale1);
title('(e) Recovered image ①')
subplot(3,2,6);image(Y2/scale2);title('(f) Recovered image ②')
hold off

7

画像復元の実行結果公開されているICAプログラム■ 情報量最大化算法(Infomax)
http://www.cnl.salk.edu/~tony/ica.html■ 行列同時対角化算法(JADE)
http://sig.enst.fr/~cardoso/jointdiag.html■ 不動点法(FastICA)
http://www.cis.hut.fi/~aapo/■ 多種多様な計算アルゴリズム
http://www.bsp.brain.riken.jp/ICALABＩＣＡＬＡＢの使い方■ ＭＡＴＬＡＢ（バージョン5.3以上）のＣＯＭＭＡＮＤ ＷＩＮＳＯＷＳでＩＣＡＬＡＢを入力すれば，右側図な画面を呼び出せます．
■ 先ず，信号源データを入力して，周波数を選択する上に，
select channelsを利用して，処理したいchannelsを選択するのは可能です．
■ mixing matrixの部分を利用すれば，信号源データを混ぜることはできます，更に雑音を添加するも可能です．

結果の可視化と復元

アルゴリズムの選択■ ＩＣＡＬＡＢの入手には，理化学研究所脳研脳信号処理研究チームのＨＰ
http://www.bsp.brain.riken.jp/ICALAB/でダウンロードするのは可能です．
■ ＩＣＡＬＡＢのインストール方法についでは，ＭＡＴＬＡＢ主メニューの中，バスの設定のボタンを押して，解凍されたＩＣＡＬＡＢのフォールダーをバスに増加すれば，インストール完了です．

アルゴリズムの調節と前処理■ アルゴリズムの調節には，
adv.optionsボタンを押せば，選択メニューが出られます（次のページ）．自分のアルゴリズムを利用するも可能です．
■ アルゴリズムの調節

Ⅲ_1 ICAによるMEG信号処理
・ 脳計測について・ 信号源と雑音源

２．無相関化及び雑音低減方法■ 独立成分を復原する場合は，do deflationのボタン押せば，選択された独立成分を復原されます．図の右側のように表示してます．

■ ＩＣＡＬＡＢツールボックスは理化学研究所脳研脳信号処理研究チーム開発され，画像処理，信号処理，EEG/MEG 信号解析からなる3つの独立した
MATLAB のソフトウエアーパッケージである．ICALAB には多くの強力なアルゴリズムが実装されている．ICAの他に，線型予測性やスムースネスや希釈さになどの基準に基づいたブラインド信号源分離などが実装されている．
一番の強みは，様々な長所や短所をもつ強力なバッチアルゴリズムを一同に公開している点である．

■ 信号源の入力と混合を完了されたら，信号処理を行います．ＩＣＡＬＡＢの中でおよそ２０
種類のアルゴリズムがあります．AMUSE, SOBI(SOS
algorithms), SANG, JADE,SIMBEC (HOS algorithms)な
ど処理方法の詳しい説明には，ＨＰで書いております．

１．脳磁界(MEG)■ 計算完了して，結果の表示方法はいくつがありますのに，ここには，deflationのボタン利用する際，結果としで図の左側のような表示せれてます．

ＩＣＡＬＡＢツールボックス

・ 主成分解析法
・ 因子解析法
・ アルゴリズム

３．MEG信号処理の事例
・ Phantomデータ解析例
・ ヒトの聴覚刺激による反応の動態情報の可視化■ 前処理の選択（IIR/FIR filter）

脳計測技術
1．コンピュータ断層画像(CT: Computer Tomography)
・ X射線を与え，生体外部からの計測により，生体内部の構造を画像で再構成する計測法である．
・ 画像処理にはこれまでの信号処理手法が多く適用されている．

2．陽電子崩壊断層画像(PET: Positron Emission Tomography)
・ 体内に放射性物質を注入し，放射線の計測により，体内の断層像を再構成する計測法である．
・ 独立成分解析による断層画像処理の報告がある．

3．機能的核磁気共鳴画像(f-MRI: functional-Magnetic Resonance Imaging)・ 強い磁場を与えて，血液の中にある酸化ヘモグロビン及び還元ヘモグロ
ビンの量を計測する方法である．
・ 独立成分解析によるfMRI画像の分解・再構成が多数報告されている．

4．脳波(EEG: Electroencephalography)
・ 神経細胞の活動により流れている電流を，頭皮上に置く電極を通して記録する．
・ 独立成分解析によるデータ解析の典型例として盛んに研究されている．

5．脳磁界(MEG: Magnetoencephalography)
・ 神経細胞電流が引き起こした超微弱磁場を， SQUID (超伝導量子干渉素子)を通して記録する．
・ 独立成分解析によるデータ解析の典型例として盛んに研究されている．
武田常広： 「脳工学」 コロナ社 2003.4

8

脳磁気の計測原理

MEG計測システム

MEG計測装置■ 神経細胞が活動すると，電流が流れる．
■ 群内電流は単一電流源として観測される．
■ 電流の変化により磁場が生じる．
■ 右ねじ法則に従って磁場の検出が可能である．

■ 脳磁場が超微弱であるため，できるだけすべての磁場エネルギーを検出する．その場合検出コイルのマグネトメー
タ(magnetometer)を利用する．
■ 一方，脳磁場に比べて強い雑音的な磁場(例えば，地球磁場)の影響を低減するため，磁場の空間的な変化(傾き=勾配)，を検出コイルのグラディオメー
タ(gradiometer)を利用する．

・ 地球の磁場強度(磁束密度)： 10・ 脳の磁場強度(磁束密度)：

−4

～ 10

−5

10 −12 ～ 10 −13

Tesla
Tesla

・ SQUIDの感度 < 10 −15 Tesla

カナダＣＴＦ社製ＭＥＧ計測装置 http://www.ctf.com/home.html

MEG実測データ例

ある成分はx(t ) = As(t ) + e(t ), t = 0,1,!
Math. definition

Specificationx = [ x1 , ! , xm ]T

CTF 64-Channel SQUID 観測信号

近隣センサの信号の強度少なくても2個以上のセンサを寄与すれば、信号源といい．
多くても1個のセンサをだけ寄与すれば、雑音といい．

Noise sources :Power supply
interference

MEG signals

EF
BG

noise

Artifacts sources

Brain sources

・ チャンネルごとの計測信号

ICAによるMEG信号解析

A = (aij )信号源と雑音源

MEGの計測環境では雑音の対策が非常に重視される．例えば，非磁性材料で作られた磁気シールド室やMEG計測装置自体などを使用したり，また磁気検出コイルにも多次グラディオメータの使用などの色々な工夫がされている．ただそれだけ実施しても，原理上でも実測した信号から見ても限界があるので，最終的な手段としてデー
タ処理を工夫することがどうしても必要になる．

・ 13Ｈｚ正弦波が信号源として計測されている場合：

・ 64チャンネル全体の計測信号

脳信号処理技術が必要データから雑音の識別例
実測データ

センサの配置図
センサの配置図

Averagedi-番目の信号源からj-番目センサへの距離に関する係数行列s = [ s1 , ! , sn ]T

共通成分 (信号源：雑音源を含む)

xi = ai1s1 + ! + ainsn + ei

独自成分

e = [e1 ,!, em ]T

独自成分 (雑音)xj = aj 1s1 + ! + ajnsn + ej

独自成分

近隣センサ信号共通成分

9主成分解析による前処理

信号解析の方式

雑音がある場合の信号分離

前処理、ICA、後処理
Data model and its covariance matrix :X( m × N ) = A( m × n )S( n × N ) + Ξ( m × N )■ 雑音がある場合の復元信号と原信号の関係：

C = AAT + Ψy (t ) = WAs(t ) + Wξ (t )重み行列の推定が難しい．
Decomposing
super-& subGaussian sources.

noise power,
correlation,
dimensionality.

■ ICA処理の前に，予め雑音の取り除くあるいは雑音パワーを低減させる方法は有効（現実的）と考えられる．

L(A) = tr[AAT − C]2

Cost function :Reducing

Projection
Localization

Eigenvalue decomposition :Pre-processing

主因子解析による前処理

計算の流れ

Ψ >> 0Eigenvalue decomposition :Factor analysis:

Pseudo-inverse matrix :−1 −1
雑音の共分散行列 Ψの初期値をΨ(0) = diag (C ) により与える

C のn個の最大固有値 Λ n と固有ベクトル U n を求める

ˆA
ˆ T = UnΛnUTAnPrincipal-Factory
Analysisˆ = diag (C − AA T )
Ψ
ˆ Ψ
ˆ A
ˆ] A
ˆ Ψ
ˆ
Q = [A

Result II

T

−1

−1

T

Result I

交差確認法による次元の縮約

観測信号 X の共分散行列 C を計算

L(A, Ψ) = tr[AAT − (C − Ψ)]2

Cost function :Post-processing

ICAˆ = Un ΛnˆAˆ T = UnΛnUT
A
A
n
1
UTn y
Λn
E[zz T ] = Λnz = Q PCA y =Principal scores:

Case II : Pre-whitening with noise power reduction

Ψ→0Case I: Pre-whitening by standard PCA⌢⌢

The cross-validatory technique
1

2

3

4

5

1 --- 4

5

determine data characteristics4

verify1

A 伝送路行列の推定値 A を A = U n Λ n2 により推定する

⌢
⌢ ⌢T
⌢
Ψ の推定値 Ψ を Ψ = diag (C − AA ) により推定する

General Error =}Averaging

Data sets :1 ,…, i-1, i+1 ,…, K

1 K∑ tr[Ci − ( A j ATj + diag (Ψˆ j ))]2
K i =1i≠ j

A novel criterion for estimating the source number :

−1
収束する

E[zz T ] = Λn + CΨCTˆ) =Error ( n

1 K
ˆ ( nˆ ) ) − diag ( Ψˆ ( nˆ ) )]2∑ tr[diag (Ψij
K i =1ˆ TΨˆ −1 Aˆ −1 ]−1 Aˆ TΨ
ˆ −1より計算する
A の一般化逆行列 Q を Q = [AJianting Cao : ``A robust approach for ICA with high-power noisy measurements and unknown number of sources'', in Neural Information Processing, Edts:
Liming Zhang and Fanji Gu, ISBN 7-309-03012-5, Fudan University Press, vol. II, pp.1012-1017, Sept. 2001.

雑音分散の推定

次数縮約の最適値
1 K
ˆ ( nˆ ) ) − diag ( Ψˆ ( nˆ ) )]2ˆ) =Error ( n∑ tr[diag (Ψij
K i =1diag(Ψ)SNR = 10log

E[ si2 ](dB)
E[ξi2 ]

ˆ)diag(Ψ

1.2499-0.97

1.2594101.0902 0.2577

-20.1

5.9

101.0757 0.2877

0.00010.00011.00780.9823

40

40

-0.03

0.08

0.1125 0.0073

0.9931 1.0120

PCAとFAの比較
Standard
Pre-whitening

RobustPre-whitening

10

ICAアルゴリズム

Phantom 実験

後処理：推定信号の射影と位置推定

Typical ICA algorithms(1) KL information based ICA algorithms （Amari-Cichocki–Yang; Bell-Sejnowski）(2) Fixed-point based ICA algorithm (Hyvarinen-Oja )(3) JADE based ICA algorithm (Cardoso-Souloumiac)Applied ICA algorithm(1) Source projection:
Estimated sources (ICs):
Projecting individual IC
to sensor space:

Natural gradient algorithm (Amari-Cichocki-Yang; 1995)

ΔW = η[I − ϕ (y )y T ]W

A single current dipole was put in the phantom head where thephysiological saline solution was filled.

sˆ(t ) = y(t ) = WQx(t )Source : Sine wave
current; f = 13Hz, p = 10 µAxˆ (k ) = aˆ i sˆi (k )Sensors : 64

(2) Source localization:κˆi ≤ 0α −1ϕ i ( yi ) = αλα sgn( yi ) | λα yi |Duration: 4’16”
Trial number : 256lsκˆ i > 0

The equivalent current dipoles are found by minimizing
the error function between solution of the forward modeland measured signals.

Sampling rate : 250 Hz.

tria(1 + β ) yi
ϕ i ( yi ) = 2
yi + β / λ2β256t-distribution based function (Cao-Murata; 1999)

Samples : 250 per trial.

We use spatio-temporal fit method is developed by CTF system Inc..

1 sec.

同期平均加算

平均処理 vs. ICAによる単試行データ処
理

前処理＋ICA

Averaged

256

256

tria
ls

tria

ls

Averaged

Error=3.61%

ICA

Cumulative error 2.7 %聴覚誘発(AEF)実験

ICAによる単試行データの処理

聴覚誘発の平均化処理

単試行処理の方式：
330 fT
Averaged AEF

630

Sensors : 64
Duration: 6’19”
The number of trials : 630
Sampling rate : 312.5 Hz.
Samples : 188 per-trial

tria
ls

The data set of auditory evoked
fields (AEF) experiment wasrecorded by testing normal adult.
Ave.

N100 component

N100

Alpha-wave 11Hz

630

tria
ls

平均化処理についてP100 component
Environment noisewith a high frequency

11結果比較I：刺激による反応強度

結果比較II：位置推定z

Two dipoles for averaged data

First Stimulus

Decomposed

ls

63

0

tria

330 fT

184 fT

IC 2

IC 3

4th. trial

721 fTIC3 : left component

721 fT

Relationship :184 fT

4th. trial
516 fT

330 fT

IC1 : right component

1[ xˆ , yˆ , zˆ ]TR
IC 3 = [ −1.2,3.9,4.7][ xˆ , yˆ , zˆ ]LAve,R ≈1st.trial

201 fT[ xˆ, yˆ , zˆ]RAve = [1.7,−5.6,6]y
IC 1

1st.trial[ xˆ, yˆ , zˆ]LAve = [−0.6,4.9,6]

x

ICA

結果III：反応強度の変化と現象観察[ xˆ , yˆ , zˆ ]TrIC11 = [0.5,−1.8,7.3]

1 T∑ [xˆ, yˆ , zˆ](ICi )1, IC 3
T i =1(1) Relationship:[ Amp ]LAve
,R ≈

1 T
∑ [Amp ](ICi )1, IC 3
T i =1(2) Some phenomena `familiar', `novel' and `loss of attention' can be observed.

Ⅲ_2. EEGによるBCIBCIシステムとは

非侵襲的なBCI■ 脳の電気的な活動信号を利用して外部装置を制御するシステム。
■ もともと人間の手足で動かしたものを，「脳」の念じることで動かすことが可能であるか？
■ 意義：生体内に新たな「通信チャンネル」を創出することになる。

・ EEGに基づいたBCI生理実験
・ EEGデータ採集
・ EEGデータの前処理と雑音除去
・ 特徴抽出と分類
・ 装置制御インタフェース
・ 装置の制御■ BCIシステムの三大要素：●信号採集：(非)侵襲, EEGなど●信号処理：雑音除去，特徴抽出等●装置制御：車椅子，携帯，ロボット，
PC，スイチなど
104

BCI実験課題(P300)

105

P300 聴覚刺激例■ P300反応 （oddball課題）

BCI実験課題(SSVEP)■刺激：1KHz と2KHz純音刺激■刺激呈示頻度：2 vs.8
■試験者： 2KHz純音を数える■ P300とは2種類以上の感覚刺激による引き起こした陽性電位の反応■反応時間：約250～500msec■刺激：聴覚・視覚・体性感覚・臭覚・味覚などの感覚刺激■ 定常的視覚誘発電位 （Steady-State Visual Evoked Potential：SSVEP）■ SSVEP とはフリッカ刺激による視覚野へ引き起こした電位(Morgan et al., 1996)．
■刺激：フリッカ(flicker)：ディスプレイに文字などをちらつかせる。
■刺激周波数：10—40Hz■反応周波数：刺激周波数に対応する■ 目標に注意力を払い、誘発電位の振幅が増大■刺激呈示頻度：2 vs.8(低頻度の刺激を被験者に選択的に注意させる)■刺激間隔：約0.5秒■一試行間隔：1—2 or 3 秒

106■加算回数：20—30回■FIRフィルタ：1-40Hz■試験者： 個人差がある

107
http://www2.oninet.ne.jp/ts0905/human/p300.htm ビギナーのためのERP入門

108
Tsinghua University

12

BCI実験課題(MI)Ⅲ_3. EEGによるBDD
1. 脳死の定義
全脳機能が完全,不可逆で失う．
２．脳死判定の項目と流れ■ 運動想起 (Motor Imagery: MI課題）■ MIとは左右手足の運動想起による反応■反応場所：C3, C4

右図の❶ → ❺３．問題点：
1) 呼吸器外しのリスク
2) 平坦脳波の長時間確認
４．EEG予備システム導入の利点：
生存者の早期発見，誤診防止
５．EEG予備システムの構成：

Primary Motor Cortex

EEG recording
・ Recording place : Shanghai Huashan Hospital in affiliation with Fudan University
・ Recording environment : ICU, at the bedside of patient・ Recording duration : June 2004 to March 2006
(JSPS & NSFC joint research Project: Prof. Fanji GU, & Jianting Cao)・ Recorded number : 35 patients with 64 sessions (patients in a deep-coma state)
・ Equipment : a portable NEUROSCAN ESI system
・ Electrodes : 6 electrodes, GND, A1+A2 Ref.
・ Sampling Frequency : 1000 Hz
・ Filtering : none
・ Resistances : < 5k, or < 8k
F3

F4
F8

Primary Somato Sensory Cortex

Fp2

GND

Fp1

A2
109

110

Jiangting Cao and Zhe Chen: "Chapters 15: Advanced EEG Signal Processing in Brain Death Diagnosis", in the Book "Signal
Processing Techniques for Knowledge Extraction and Information Fusion", pp.275-297, Springer. March 2008

Patient (A)’s symptom and EEG recording
・ 18-year-old male patient (May 20, 2004)
・ with a primary cerebral disease・ Coma test : lost consciousness, in a deep-coma state・ Pupil test : dilated.
・ Brainstem test : no reflexes.
・ Respiratory machine was used, but Apnea test hasn’t done
* The symptom was very similar to a brain-death case.

Patient (A): (June 22, 2004)
Result I: used FA+ICA algorithms

112

A１
111

Patient (B): (Oct. 19, 2005)

The symptom was very similar to a brain-death case.
Result Ⅱ: used FA+ICA algorithms and others
41 sec.

Decomposed 3 components in the time and/or frequency domain.

Under the family’s agreement,
EEG pre-examination was madein ICU at the bedside of patient.
・ Three times EEG recording,each one approximately
in 5 min.
・ A time window of 10 sec.
(190s-200s) EEG data.

111

F7

42 sec.

Heart beating

A 12.5 Hz alpha-wave component has been extracted
112

Jianting
Cao "Analysis of the quasi-brain-death EEG data based on a robust ICA approach", Springer-Verlag, ISBN113
3-540-465421, Lecture Notes in Artificial Intelligence, Vol. 4253, No. 3, Springer, pp.1240-1247, Oct. 2006.

114

114

Summary of ICA result for 35 patients

1. ICA result is almost identical to the clinical diagnosis.
2. The patients are categorized into two groups.
3. Coma group: 19 patients, quasi-brain-death 16 patients (18cases)
115
115

13