第3章 音声・画像の信号分離
3.1 実験目的
電子情報処理分野において、音声、画像、通信データ等の信号処理が主な課題としてしばしば取り上げられる。本実験の目的は、近年発展の著しいブラインド信号処理(Blind SignalProcessing)分野の応用として信号分離(Source Separation)問題を取り上げ、MATLAB のプログラミングを通して実際に信号分離を行うことである。また、音声・画像データ等の取り扱い方を習得する。更に混合信号分離問題への数式化を行い、その数理モデルの解析および計算アルゴリズムの設計方法について学ぶ。

3.2 実験原理
3.2.1 混合信号の分離方法
1. 信号分離問題の背景

電子情報工学などの分野では、複数の信号成分が互いに混じり合っている混合信号から、目的とする信号成分のみを取り出すという問題がしばしば現れる。例えば、何人かが同時に発話した音声をマイク（センサー）によって収録した混合信号から、各話者に対応する音声信号に分離する問題がある。またノイズに埋もれた画像から元画像に復元する問題もある。
これらの問題は、近年、ブラインドセパレーション(Blind Source Separation)の枠組みでその解決が試みられている。また、90 年代の中期頃に独立成分解析(Independent ComponentAnalysis)が提唱され、混合信号分離問題の有力な一つの手法として注目されている。
2. 信号分離問題の定式化

信号分離の目標は、混合された信号を個々の原信号（または信号源）に分離することである。この信号分離の問題を効率よく検討するためには、まず数式で問題を記述すること（すなわち問題の定式化）が必要である。今、２人同時に発話された音声信号の分離問題について考えてみよう。ある場所で A さんと B さんが同時に喋っている。ここで、A さんと B さんの声を信号源としてそれぞれ s1 と s2 とする。その発話された音声信号を記録するために、二つのマイクロフォン１と２が設置されていると仮定し、各々のマイクロフォンによって記録された信号を x1 と x2 とする。また、信号源から各マイクロフォンまでの伝達特性は、係数 a11、a12、a21、a22 のみを用いて簡略化されるとすると、マイクロフォン 1 で記録された混合信号 x1 は、式(1)で表現することができる。

x1 (t ) = a11 s1 (t ) + a12 s 2 (t )t = 1, 2, !

(1)

上式において、係数 a11 と a12 は信号源 s1 と s2 からマイクロフォン 1 までの物理距離を示していると考えればよい。すなわち、信号源とマイクロフォンの距離が近ければ、係数の値が大きくなり（つまり、音声信号の振幅が増加する）、信号源がマイクロフォン信号 x1 を占有している度合いも大きくなる。ここで、t=1,2,… は離散時間とする。
同様にして、マイクロフォン 2 で記録された混合信号 x2 は式(2)で表現できる。

3-1

x 2 (t ) = a 21 s1 (t ) + a 22 s 2 (t )(2)

以上、音声分離の問題を簡単な数式(1)と(2)で記述することができた。つまり、今我々が解こうとしている信号源の分離問題は、上述した問題の定式化により、式(1)と(2)の s1 とs2 の解を求める、という数式上の問題に帰着することができる。しかし、実際には係数 a11、a12、a21、a22 は未知であり、式(1)と式(2)のみを利用して信号源 s1 と s2 を求めることは不可能である。そのため、何らか方法を用いて s1 と s2 の（少なくとも）近似解を求める必要がある。これについては次節に詳しく説明する。
簡潔のため、式(1)と式(2)をまとめて次式のように行列形式で書き直す。

x(t ) = As(t )(3)

ここで、マイクロフォン信号は縦ベクトル x(t ) = [ x1 (t ) x 2 (t )]T（T は転置の記号である。）、

a12 ùéaú となる。
ëa 21 a 22 û

信号源は s(t ) = [ s1 (t ) s 2 (t )]T 、また、伝達特性を表す係数行列は A = ê 11

例題 1：ホームページ(http://www.sit.ac.jp/user/cao/Teaching/Teaching.htm 第６章
音声データ)から「moshimoshi」と「birdchirps」という音声信号（各々信号源信号として用いる）をダウンロードし、ハードディスクD:にあるフォルダ MYWORKに保存する。次のプログラムを実行し、実際に元の音声と混合された音声信号を視聴して確かめよう。また、以下のプログラムでは、二つのマイクロフォン信号にあたる混合信号が保存される。

プログラム samp1.m% Original speech signals and mixed signals% filename: samp5_1.m% Jianting Cao, Aug. 10, 2003
% Workspace 変数をすべて削除
clear all% 音声信号の読み込み、信号源の生成
load birdchirps1;
load mosimosi1
s=[birdchirps1;mosimosi1];% 音声信号を試聴する
clc;input('鳥の鳴き声 … ')sound(birdchirps1);
clc;input('もしもし … ')sound(mosimosi1);% randn()を用いて混合係数を生成
A=randn(2,2);3-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 ②')
上記のプログラムを実行すると、元の音声信号と混合された信号の波形が次のように表示される。

3. モデルの解析

3-3

前節では「マイクロフォン信号 x(t)は既知であるが伝達特性係数 A が未知である、という条件下で、源信号 s(t)の近似解を求めるには何らかの手法が必要である」と述べた。
本節では、その近似解を求めるために、下図に示すようなニューラルネットワークを利用することを考えてみる。
今 ニ ュ ー ラ ル ネ ッ ト ワ ー ク の 入 力 を x(t ) = [ x1 (t ) x 2 (t )]T と し 、 出 力 をw12 ùéwy (t ) = [ y1 (t ) y 2 (t )]T とする。また、ニューラルネットワークの重み行列を W = ê 11úë w21 w22 ûとする。

上図から、ニューラルネットワークの入出力の関係は次式で表現される。

y (t ) = Wx(t )(4)

式(3)を式(4)に代入すると、式(5)が得られる。

y (t ) = WAs(t )(5)

ここで、式(5-5)において、係数行列 W は自由に調整できることを仮定している。
もし、偶然 W の値が W = A -1と調整された場合、式(5-5)から y (t ) = s(t ) となり、分離信号 y (t ) は信号源 s(t ) と等しくなる。つまり，源信号はマイクロフォン信号から忠実に分離されたことになるが、実際、常にこのような場合になるとは限らない。
今、もう少し現実的に、仮に W の値が W = PDA -1となることを考えよう。ここで、 P は順番置換行列(permutation matrix)と呼ばれ、各行と各列に 1 に等しい成分が 1 個ずつあり、

é1 0ù
é0 1ù
かP=êúú
ë0 1û
ë1 0û

その他の成分が 0 の行列である。例えば、2 次元の順番置換行列は P = êのどちらかである。

3-4

éd 0 ùú または
ë0 d 2 û

また、 D は対角行列(diagonal matrix)で、2 次元の場合には D = ê 1

é0 d 1 ùD=êú と表現される。
d0ë 2û
上記仮定を考慮して、 W = PDA -1を式(5)に代入すると、分離された信号 y (t ) と信号源s(t ) との関係は

é y1 (t ) ù é1 0ù éd1 0 ù é s1 (t ) ùê y (t )ú = ê0 1ú ê0 d ú ê s (t )ú
ûë
2 ûë 2ë 2 û ëû(6)

é y1 (t ) ù é0 1ù é0 d1 ù é s1 (t ) ùê y (t )ú = ê1 0ú êd 0ú ê s (t )ú
ûë 2
ë 2 û ëûë 2 û

(7)

またはの二通りであることが分かる。
ここで、式(6)に着目すると、 y1 (t ) = d1 s1 (t ) また、 y 2 (t ) = d 2 s 2 (t ) であり、分離された信号 y (t ) の大きさは対応する対角行列 D の成分によってスケーリングされている可能性がある。また、式(7)から y1 (t ) = d1 s 2 (t ) 、 y 2 (t ) = d 2 s1 (t ) であり、復元された信号の大きさは D によるスケーリングされるだけでなく、その順番が必ずしも原信号と一致しない。このことを各自実験中に観察してみること。
４. 学習アルゴリズムの設計

前節では、重み行列の候補が W = PDA -1または W = A -1 であれば、混合信号が個別に分離ができる、ということを述べた。しかし、実際には、係数 A が未知であり、また重み行列 W を構成する係数の可能な組み合わせは無数にあるため、上述の理想的な解を求めることはほとんど不可能である。本節では、近年、情報理論を駆使して考案された混合信号の分離を行うための自動的な重み W の計算手法について説明する。なお、詳細については、あらかじめ統計情報学を習得している必要があるので、ここでは計算手法の枠組みについて大まかに説明する。
まず、源信号が互いに統計的に独立（ p s (s) =nÕ p (s ) 結合分布の確率密度関数が周辺iii =1分布の確率密度関数と等しい）であるとし、独立信号ｓが A によって伝達され、ｘとして観測されるとする。基本的に、独立成分解析(ICA：Independent Component Analysis)とは、観測信号ｘのみを用いて、出力信号ｙが限りなく統計的に独立になるような、ニューラルネットワークの再構成することである。すなわち、出力信号ｙの結合分布の確率密度関数がその周辺分布の確率密度関数と等しくなったとき、つまり、
np y (y ) = Õ pi ( y i )
i =1

3-5

を満たすとき、統計的に独立した源信号ｙが推定できた、原信号を復元することになる。
この統計的な独立性を計る情報量としては次式に定義される Kullback-Leibler 情報量が用いられる。

D(y | W ) = ò p y (y ) log np y (y )

Õ p (y )idy(8)

ii =1

ここで、式(5-8)において、ｙの結合分布の確率密度関数とその周辺分布の確率密度関数の関係が p y (y ) =nÕ p ( y ) となれば、Kullback-Leibler 情報量が D(y | W ) = 0 となる。
iii =1

この情報量が最小になるように重み行列の学習アルゴリズムは¶D(y | W )= 0 で導くこ¶Wとができる。最終的に導かれた逐次型(on-line)の学習アルゴリズムはW (t + 1) = W (t ) + DW (t )(9)

= W (t ) + h (t )[I - j (y (t ))y T (t )]W (t )となる。上式において、 I は単位行列であり、 h (t ) は学習率、また j (×) は Score 関数と呼p! ( y )d
ばれる。Score 関数 j (×) の成分 j i = log pi ( y i ) = - i i は各信号の確率密度関数
dy ipi ( y i )pi (×) に依存することが判っている。
関数 j (×) の近似方法は幾つかが提案されているが、一般的には信号分布の尖（とがり）E{ y 4 }- 3を計算し、 j (×) を選択することが多い。
(E{ y 2 }) 2
例えば、音声信号( k > 0 )の場合には関数(10)j ( y ) = tanh( y )が用いられる。また、画像や通信信号( k < 0 )の場合には、次の関数を使用する。
(11)j ( y) = y 3度を表す kurtosis k ( y ) =例題 2：例題1で保存された混合信号を利用して音声分離を行う。但し、学習アルゴリズム
は式(5-9)を利用し、 j ( y ) = tanh( y ) を用いる。確認のため、源信号の尖度

k (s) =

E{s 4 }
- 3 を計算してみよう。
(E{s 2 }) 2
プログラム samp2.m% Separation of the mixed speech signals% filename: samp5_2.m% Jianting Cao, Aug. 10, 2003
% Workspace 変数をすべて削除
clear all

3-6

% 混合信号の読み込み
load x;% 学習回数、学習率、重みの初期値の設定
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,:)));

５. 画像の復元

前述したように、画像工学・医用画像処理の分野においては、複数の画像が重なったものからある特定の画像を復元することや、ノイズに埋もれた画像から原画雑に復元するという問題がよく現れる。次の例題3では、前節に述べた信号分離の方法を用いて画像の復元を試みる。
例題 3：ホームページ(http://www.sit.ac.jp/user/cao/Teaching/Teaching.htm 第６章画像データ)から「testpat1.tif」と「rice.tif」という2つの画像データをダウンロードし、ハードディスクD:にあるフォルダMYWORK に保存する。次のプログラムを十分に理解したうえで実行し、画像復元の効果を確かめよう。なお、プログラム中分からないMATLABコマンドがあれば、各自Helpで調べること。

プログラム samp3.m% Recovered the mixture of images

3-7

% filename: samp5_3.m% Jianting Cao, Aug. 10, 2003
% 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

3-8

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
上記のプログラムを実行すると、次図のようになる。

3. 3 実験課題
1. 一般化された瞬時混合モデルのブロック図は次のように示される。ここで、混合信号の数は信号源の数より多い(m>n)とする。

3-9

s1 (t )s2 (t )

x1 (t )x2 (t )!

sn (t )!

xm (t )!

yn (t )AW
y1 (t )y2 (t )

このブロック図に基づいて、1) 信号源、2) 混合行列、3) 混合信号、4) 重み行列、5)分離・復元信号、及びそれらの関係を数式で正しく記述せよ。また、特別な場合、すなわち、m=n=2で混合行列の値がA=[1 1;１１]の場合に、信号分離ができるかどうかについて論じよ。
2. 例題 1 と例題 2 において、一つのノイズ信号源(randn 関数によって生成)を加えた計 3個の混合信号を分離するようなプログラムを作成せよ。なお、レポートにプログラムと実験結果を示せ。参考のために実行結果の一例を下図に示す。時間的に余裕がある場合はプログラムを function M-ファイルで作成せよ。

3-10