はじめに
前回の「グラフで理解する独立成分分析」の続きです。
独立成分分析を使って混ざった信号を分離します。
scikit-learnに全く同じことをやったチュートリアルがありますが、
ここではscikit-learnなしでnumpyとscipyだけでも実装します。
参考
scikit-learnを使ったチュートリアルです。
実際のデータに使うときはライブラリに頼ったほうが簡単そうです。
scikit-learn Blind source separation using FastICA
データ化学工学研究室 独立成分分析 (Independent Component Analysis, ICA) ～PCAの無相関より強力な ”独立” な成分を抽出～
ブラインド音源分離の流れ
今回は2つの信号を混ぜることを考えます。
2つの信号が非ガウス分布で独立であるなら、独立成分分析を行えば分離できます。
1. 信号の混ざり方を決める混合行列
混ぜる前の信号を$S$とします。
S = \left[
\begin{matrix} s_1 & s_2 \end{matrix}
\right]
2つの信号をどんな割合で混ぜるかを決める行列を混合行列と言い、ここでは$M$とします。
M = \left[
\begin{matrix}
m_{11} & m_{12} \\
m_{21} & m_{22}
\end{matrix}
\right]
混ざった後の信号$X$は
X = SM
一般的に$X$が測定データです。
$M$の逆行列が求まれば元の信号がわかることになります。
S = XM^{-1}
2. 独立成分分析
独立成分分析では、
- 無相関化
- 分散の正規化
- 非ガウス性が最大になるように直交変換
の順に行います。
イメージを掴むために、以下のデータに操作を加えて可視化します。
2-1. 無相関化
主成分分析を使って無相関化します。データの平均値は予め0にしてあるとします。
手順は、分散共分散行列を計算、その固有値分解を行い、データに固有ベクトルをかけます。
cov = {\rm COV} (X) \\
UWU^T = {\rm SVD} (cov) \\
X_{2-1} = XU
$X$の添字は節と同じ数字にしておきます。
2-2. 分散の正規化
分散を正規化します。
固有値$W$は対角行列なので、逆行列は対角要素の逆数が並ぶ対角行列です。
それから対角成分の平方根を取ってかけます。
X_{2-2} = X_{2-1} W^{-1/2}
2-3. 非ガウス性が最大になるように直交変換
まだ回転方向に自由度があります。
非ガウス性が最大になるような直交変換を行う回転行列$Z$を探します。
どれくらいガウス分布から離れたかの指標には尖度などが用いられます。
S = X_{2-2}Z
これで信号の分離ができました。
3. 混合行列を計算する
測定データの中に分離できた信号がどんな割合で含まれているのか知りたくなります。
それは混合行列がわかればいいことになります。
先程までの操作をまとめます。
S = XUW^{-1/2}Z
測定データは混合行列$M$によってこのように生まれました。
S = XM^{-1}
両者を比較して、
M^{-1} = UW^{-1/2}Z \\
\therefore M = \left( UW^{-1/2}Z \right) ^{-1}
これで混合行列が求まりました。
ブラインド音源分離の実装
NumpyとScipyを使った実装
先程の式に従って実装します。
0. 使用するモジュールと関数
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.decomposition import FastICA
円を描くときに使います。
def circle(n,r=1):
rad = np.linspace(0.0,2.0*np.pi,n)
return np.asarray([np.cos(rad),np.sin(rad)]).T
データを回転させるときに使います。
def rotation_matrix(rad):
return np.array([[np.cos(rad),-np.sin(rad)],[np.sin(rad),np.cos(rad)]])
1. 混ざった信号を作る
2つの信号を作ります。最終的にこれが得られれば良いです。
t = np.linspace(0.0,1.0,5000)
s1 = np.sin(2.0*np.pi*3.0*t)
s2 = np.sign(np.sin(2.0*np.pi*4.0*t))
S = np.array([s1,s2]).T
S += 0.05 * np.random.normal(size=S.shape)
散布図はこのようになります。x軸とy軸が独立であることがわかります。
上の信号の混ざり具合の異なる2つの信号を作ります。
Mix = np.array([[0.2,0.4],[0.8,0.6]]) #混合行列
X = S.dot(Mix)
散布図はこのようになります。相関があるように見えます。
2. 無相関化
主成分分析を使って無相関化します。
cov = np.cov(X, rowvar=False) #分散共分散行列
U,W,_ = np.linalg.svd(cov) #固有値分解 Uが固有ベクトル Wが固有値
pca_X = X.dot(U)
#確率長円もどきを描きます
c = circle(200,r=1.0)
c = c*np.sqrt(W)
3. 分散の正規化
分散共分散行列を単位行列にします。
norm_pca_X = pca_X/np.sqrt(W)
c /= np.sqrt(W)
4. 非ガウス性の最大化
ここではガウス分布にどの程度従うかの指標に尖度を使い、総当たりで回転行列を探します。
rad_array = np.linspace(0.0,np.pi/2.0,1000)
kurt_array = []
for rad in rad_array:
rot_data = norm_pca_X.dot(rotation_matrix(rad))
kurt = stats.kurtosis(rot_data)
kurt_array.append(kurt)
kurt_array = np.asarray(kurt_array)
fig,axes = plt.subplots()
axes.plot(rad_array,kurt_array.T[0])
axes.plot(rad_array,kurt_array.T[1])
axes.set_xlabel("Rotation angle rad")
axes.set_ylabel("kurtosis")
上の2つのプロットを足していいのか知りませんが、足してうまくいくので良しとします。
Z = rotation_matrix(rad_array[np.argmax(np.abs(np.sum(kurt_array,axis=1)))])
S = norm_pca_X.dot(Z)
分離された信号$S$が求まりました。
fig,axes = plt.subplots()
axes.plot(t,S.T[0])
axes.plot(t,S.T[1])
5. 混合行列の計算
M = np.linalg.inv((U/np.sqrt(W)).dot(Z))
分離できた信号とこの混合行列を使って、計測データが再構成できるか試します。
np.allclose(X, S.dot(M))
>>> True
無事に再構成できました。
scikit-learnを使った実装
公式チュートリアル通りです。
ica = FastICA(n_components=2)
S_ = ica.fit_transform(X)
M_ = ica.mixing_
たった3行で済みます。
fig,axes = plt.subplots(figsize=(6,6))
axes.scatter(S_.T[0],S_.T[1],alpha=0.1)
axes.set_aspect("equal")
np.allclose(X, np.dot(S_,M_.T))
>>> True
こちらも同様の結果になります。
まとめ
独立成分分析を使ったブラインド音源分離を順を追って実装しました。
混ざる前の信号が非ガウス分布で独立であるとき、混合後の信号から分離が可能であることが示せました。
実際のマイクを使って試してみたいと思います。