# コンピュータで再現するカクテルパーティ効果

　今回の記事では数値解析用ソフトウェアとしてよく知られるMATLABを用いて，コンピュータに**カクテルパーティ効果**を再現するためのスクリプトを紹介します．正確に言うと，カクテルパーティ効果の原理の解明をモチベーションとして生まれた**独立成分分析（Independent Component Analysis: ICA）**　というアルゴリズムの紹介となっていますので，カクテルパーティ効果の考察を期待している人にとっては肩透かしな感じになってしまっていると思います（申し訳ない！）．ただ，ICA自体が「**混ざりあった波形を入れるだけで，他のパラメータは入れずに，元の分離された波形を出力してくれる**」という非常に有用な効果をもたらしてくれるので，初学者の方はその辺をフックにして読んで下さればと思います！　また，既に使ったことはあるけど原理があまり良く分からないという人にも，理解できるための一助となれると非常に嬉しいです．

## １．この記事を書いた人について
　　まずはこの場を借りて，お目汚しにならない程度で筆者について何をやっている人なのかを紹介したいと思います．
　僕はデジクリの11期生で電子工学科を卒業した後に，大学院生として生体信号を対象とした信号処理について研究しています．今回の主役となるICAについても今の研究を通して知りました．名前はすずしょーと名乗ったりとか，呼ばれたりとかしています．
　デジクリでは，有能な後輩たちに混ざってお絵描き企画やアルバム企画などに参加させてもらったり，今年度の後期限定で教卓に院生カードをかざす係をやっていたりとかしています．
　デジクリ以外では，ギターアンサンブル部でバスギターを担当した後に，今はソロギターでクラシックや日本民謡などのあまりポップでない曲を弾いたり，某技術系会社の長期インターンとしてUnityやUE4を使ったVR器材の技術デモなどを作ったりしています．
　趣味については以下の表に語れるゲームと音楽（アーティスト）を思い出せる限りリストアップしましたので，合うものがあれば気軽に話しかけてくださると嬉しいです！

| | タイトル |
| -------- | -------- |
| ゲーム | ペルソナ，ダンガンロンパ，ゼルダの伝説（時オカ，ムジュラ，トワプリ，BoW），スマブラ，バイオショック，FGO，グラブル，スカイリム，フォールアウト，シュタインズゲート，ダークソウル，デモンズソウル，Half-Life，Spec Ops: The Line，Team Fortress2，Undertale，バットマン（アーカムシリーズ），真・女神転生，Borderland2，Dishonored，紙マリオ，League Of Legend，パペッティア，Life Is Strange，Portal|
| 音楽 | 　Aphex Twin，BT, Chicane, Cosmic Baby, Cranky, David Lanz, DJ Shadow, Hybrid, Floex, London Elektricity Massive Attack, Onoken, Orbital, Paul Van Dyk, Robert Miles, Salt Tank, Sasha, Shpongle, Sven Vath, Tricky, Way Out West|

　自己紹介は以上にして，早速本題に移りたいと思います．

## ２．カクテルパーティ効果とは？
　これについては知っている人も一定数いると思われると思いますが，聞いたことのない人の為に例を使って説明したいと思います．
　まず，自分が駅前の広場で誰かと待ち合わせをしている状況をイメージしてみてください．なるべく渋谷や新宿などの人が大勢いるところが望ましいです．多分，色んな音が不規則に入り交じる感じのガヤガヤした空間が頭の中に出来上がるのではないでしょうか．聖徳太子ならば，それぞれがどういう音なのかを聞き分けられるのかもしれませんが，大抵の人はどう頑張っても聴こえてくる音が訳の分からん雑音にしか聴こえないと思います．
ですが，そんな中でも待ち合わせの相手が自分の名前を呼ぶ声だけはハッキリと聞こえると思います．で，そういう**都合の良い音を上手く聞き分けられる――聴覚の不思議な機能のことをカクテルパーティ効果**と呼びます．（更に詳しい話は日本心理学会様の[このページ](https://psych.or.jp/interest/ff-10/)を見ると良いかもです）
　この不思議な効果に関する最新の研究がどうなっているのかは分かりませんが，今の時点ではその原理は解明されていないようです．ただ，これをモチベーションとして，信号分離や特徴抽出として非常に有用な解析アルゴリズムが開発されるようになりました．それが**独立成分分析**（Independent Component Analysis: ICA）です．
　実際にICAを使って，**カクテルパーティ効果の内，”波形の分離”という機能を再現することができる**ので，早速その効果を見てみましょう．

## ３．ICAで混ざりあった波をラクラク分解
　では，ICAの効果を検証するための前準備として**MATLABで２つの波形とそれを混ぜ合わせたものを作りたいと思います**．

%パラメータ
Fs=10000;           %サンプリング周波数
time=3;             %波形の流れる時間
T=(1:Fs*time)'/Fs;  %要素ごとの時間を記入したベクトル

%信号生成
S=zeros(length(T), 2);          %波形を格納する配列
S(:,1)=sin(2*pi*800*T)';        %1列目に800Hzのsin波
S(:,2)=sawtooth(2*pi*1680*T)';  %2列目に1680Hzのノコギリ波

　まず，上のスクリプトより，3秒の800Hzのsin波と1680Hzのノコギリ波を生成し，前者を配列Sの1列目，後者を2列目に格納しました．サンプリング周波数については，基本的にこちらが使いたい周波数の中で最も高い数値の2倍を下回っていなければOK……と言いたいところですが，実はノコギリ波の周波数特性を考えると全然足りないので，結構歪な波が出力されてしまいます．そう，例えば，できた波形をプロットした以下の図のように――

えらいこっちゃ．
ただ，今回は**２つの波の形が違ってさえいれば問題ない**ので，とりあえずこのままで行きましょう．
次にこの**２つの波を混ぜたものを2つ**出力します．

 A=rand(2, 2);   %2*2の混合行列を適当に決める
 X=S*A;         %混合波形を２つ生成する

　上のスクリプトでは最初にAというランダムに決めた行列によって，２つの波形をどのくらいの割合で混ぜるのかを決めています．その次に作った信号Sの右から先程の行列Aを乗算することによって混合波形を２つ生成しています．非常に簡単なプログラムですが，行列の計算で波形を混ぜ合わせることに違和感を感じる人が出てくるかもしれません．もしその点で気になる方が居れば『ブラインド信号源分離』について調べるか，X=S*Aを[混合波形１　混合波形２]＝[800HzのSin波　1680Hzのノコギリ波]×[a1 a2; a3 a4]と書き換えて（；は改行を意味してます），掛け算の部分を計算してみるといいでしょう．
　ともかく，これで混ざりあった波形が出力されますので，その形を見てみましょう．

　混合波形１はちょっとピンと来ないと思いますが，元々の800HzのSin波から振幅が僅かに上がったり下がったりしているので，一応はもう一方の波の影響を受けています．左の混合波形２についてはわかりやすいと思います．本当は音を聴かせて，「確かに違う音になってる！」というのをもっと直感的に示すべきなのですが，HackMDでオーディオを埋め込む方法が分からなかった（もしくは無い？）ので，もし何かいい方法がある人がいれば教えてくださると非常に助かります．

ともかく，これで準備は終わりです．では，**混ぜた波形を元の波形に戻す問題**を考えましょう．
１つ目の答えは非常に簡単です．**Aの逆行列を波形を混ぜ合わせるのに使った行列Aから直接計算**して混ぜ合わせた波形Sに右からかけてやれば瞬殺で行けます．ただ，この解き方には実用上において重大な問題があります．それは，「**Aが分からなかったらどうするのか？**」という点です．
流石に理不尽すぎる指摘なのでは？　と思う方もいるかも知れません．そこで，実際の録音を考えてみてください．教室や街の中の音をマイクで拾ってみても，色んな声が混ざりあった音が聴こえはするけど，それらが「それぞれどのくらいの距離でどのくらいの音量で鳴っているのか」「それぞれどのくらいの比率で混ざり合っているのか？」というのを断定するのはできないように思えます．つまり，**今回で言う「行列Aが分からない」というのは，「それぞれがどの程度の比率で混ざり合っているのかが断定できない」と言っているのに等しいです**．実は，先程チラリと出た『ブラインド信号源推定』というのも，そういった情報がない場合の話をするための理論となっています．
では，**逆行列を元の行列から計算する以外にどうやって混ぜ合わせた信号をバラすのか？**　その答えを簡単に与えてくれるのが**ICA**です！

%自作の（実はクソガバな）ICA
Sout=ica(X);

　はい，これで終わりです．ここで特に注目してほしいのが，**すずしょー作の（案の定，クソガバな）ica君には「こんな混合波形が２つあるよー」という情報以外，何一つ情報を与えていない**という点です．では，ica君の成果ということで，出てきた配列Soutを列ごとに見てみましょう．


あれ？　出力１は確かに振幅あまり変わらなくなったけど，出力２は――


なんか違くない？？？？？？？？？？？？　
――いや，位相が逆になっただけでちゃんとほぼほぼ再現できてますよ．だって，こうやってひっくり返して反転してみると……


　　ほら！　**ICAの出力２も元のノコギリ波とほぼほぼ一緒の形になりました**！　音として聴いてみても，やはり元の音に戻っています！

## ４．ICAの核は中心極限定理にあり
　　これを見ている人達の中にもしかしたら「ついにコンピュータも黒魔術を使えるようになったのか……」って感じで戦慄する人や嘘松認定リプを送りながらバッサリと切り捨てる仕事人が出てくるかもしれません．特に前者のような人にそのまま素敵な夢を見せ続けるのもそれはそれで悪くないとは思っているのですが，この記事のコンセプト上，魔法のままにしておくのは非常によろしくないので，早速ですがタネ明かしをしようかなと思います．
　**ICAは音声の混ざり方に関する情報を一切持たない代わりに，とある確率統計において重要な定理を手がかりにしている**のです．それが――**中心極限定理**です．凄く大雑把にまとめると「ガウス（正規）分布に従って値を取らない確率変数の列を何回か取得し，数列単位でその平均を取ると，平均を取った後の確率変数列はガウス分布に従うようになる」ということを主張している定理です（意味が分からなければ先に行っても大丈夫です！）．
そもそも**ガウス分布**ってなんだっけ？　って言う人の為に以下の図を[アタリマエ！](https://atarimae.biz/archives/9850)様よりお借りしました．

> 

　　値xを取る確率の大きさをf(x)で表した時の概形が上の図の様な形になる分布です．期待値（平均値）が山の頂点，そこから標準偏差（取る値のバラつきの大きさ）分だけ離れたところが変曲点となっています．更に**ガウス分布には重要な性質が２つ程あります**．１つ目が**尖度と呼ばれるパラメータが3**となっていること，2つ目が**同じ標準偏差を持つ確率分布の中で”エントロピー”（確率分布の無秩序性の強さ）と呼ばれる指標が最大**となるということです．
――さて，ここまでの情報で実はもうICAが何をやっているかの答えは出せるのですが，肝心の**中心極限定理**についてまだどういったものか実感が沸かない方が多いのではないでしょうか．
という訳で，中心極限定理を実感するべく，以下のスクリプトを作ってみました．．

% パラメータ
sampleNum=1;    %標本数   
plotNum=10000;  %各標本のデータの長さ

% 平均標本の作成
Output=zeros(plotNum,1);    %標本平均を格納するベクトル
for n=1:sampleNum           %標本数だけ繰り返す
    temp=rand(plotNum, 1);          %一様分布に従う標本を作成
    Output=Output+temp/sampleNum;   %平均を取る
end

　　**上記のスクリプトでは一様分布に従うplotNum(=100000)個の長さのデータをsampleNum回生成しながらそれらの標本平均を取り，その結果を変数Outputに出力しています**．今回はsampleNum=1となっていますので，一様分布に従うデータがそのまま取れると思います．では，Outputのヒストグラム（それぞれの値（横軸）を何回取ったのか（縦軸）を集計したグラフ）を見てみましょう．


　　**一様分布らしい，どの値もそれなりに同じくらい取れているデータが得られた**と思います．では，sampleNum=2，つまり**2回生成してできたデータの標本平均**を取ったらどうなるのでしょう？


　なんと，一気に**ヒストグラムの形状が変わってきました**！　真ん中の方が一気に盛り上がっているように見えます．それに続いて，**sampleNumを5,10,100と増やした**時を並べてみますと――

　**sampleNum=5の時点でかなり形状がガウス分布に近づいている**のがわかると思います．sampleNum=10でほぼほぼ形ができあがり，sampleNum=100で曲線もかなりきれいに見えるようになりました．これが，恐らく数多くの初見の人が意味不明であった「ガウス（正規）分布に従わない確率（に従う）変数の列を何回か取得し，数列単位でその平均を取ると，平均を取った後の確率変数列はガウス分布に従うようになる」という説明の正体となります．
更にこれは別に一様分布に限った話でなく，**今回使った波形**についても……


　と，**混ぜればガウス分布の形に近付こうとしているのが分かります！**　なので，中心極限定理についてはこうとも言えるでしょう．

「**どんな信号も混ぜ合わせれば絶対にガウス分布に近づく，逆に一番混ざり合っていない状態ではガウス分布から最も遠い**」

　ということは，「**混ざりあった信号に対して最もガウス分布から遠ざかる様な混合行列を与えてやると，それが必然的に元々綺麗な信号をグチャグチャにさせた混合行列の逆行列となる**」ことが言えるでしょう．じゃあ，**どうやってガウス分布から遠いのか近いのかというのを調べるか**と言うと，それが先程紹介した”**尖度**”や”**エントロピー**”です！

## ５．ICA君の中を覗いてみよう

　……という訳でついに僕が最も恐れていたクソガバICA君の中身を見てみたいと思います．

function [Output,W] = ica(X)


%白色化（無相関化）
X=X-mean(X);        %期待値を0に変換
V=cov(X);           %Xの共分散行列を生成
[v, D]=eig(V);      %Xの共分散についての固有値問題（Xの特異値問題）を解く
Z=X*v/sqrt(D)*v';  %白色化の実行
c=cov(Z);         %白色化の確認

%尖度の計算
angle=(1:360)/360;  %[0~2π]までの角度を格納するベクトル
K=zeros(360, 1);    %角度ごとの尖度を格納するベクトル
for n=1:360         %各角度で計算
    a=angle(n);                 %ラジアン単位に変換
    %二次平面上の単位ベクトルの生成
    w=[cos(a*2*pi); -sin(a*2*pi)];
    K(n)=abs(kurtosis(Z*w)-3);  %ベクトルWの向きでの尖度を計算                  
end

%混合逆行列の決定
[~,n]=max(K);   %尖度が最大値を取る時の角度を取得
a=angle(n);
W=[cos(a*2*pi) sin(a*2*pi);-sin(a*2*pi) cos(a*2*pi)];%混合逆行列の推定
Output=Z*W;     %元の信号に戻ることを願うのみ

end

　まず誤解を招かないための叫びを吐き散らします．

　**この組み方は入力する信号がないと成立しません！　普通は角度ごとに尖度を求めたりせずに，ニュートン法とかもっとスマートな方法で求めます！**　

……よし，これで誤解は解けるのかな？　実は非ガウス性の尺度についても，尖度だと外れ値に弱すぎるので，「ガウス分布のエントロピー」-「自身のエントロピー」より求められるネゲントロピー，調べたいパラメータに対する計測したデータの発生確率の関数である尤度などもよく使われたりしています（それ以外にも重要な指標がありますが，ICAについてはまだ勉強不足のところが多いので後述の参考書などを参考にすると良いと思います）．
　さて，このスクリプトが何をしているのかを簡単に説明します．「**尖度の計算**」と書いているところで，1:0から0:1の割合でなるべく細かく混ぜ合わせ方を変えていく中で，**どれが一番尖度が3から一番遠下がるのかを調べています**（ガウス分布の尖度が３な為）．こうやることで，**１つ目の信号はそれで簡単に分離**できるのですが，じゃあ，**２つ目の信号はどうするのか？**　というのを考えた時に役に立つのが最初に行った処理である「**白色化（無相関化）**」というステップです．これは**１つ目の答えが分かったら，もう一つの答えはそれに直交する混ぜ方をすれば簡単に求められることを保証するため**にやっております．
　~~正直，もうMATLAB君で作ったfigureを貼り付けるの面倒くさくなってきたので止めにしたかったのですが~~この辺りも文だけじゃ上手く説明できないので，これより**点粒祭り**を開いて，それによって説明したいと思います．
　という訳で早速また訳の分からん物が出てきましたが，これ実は（**x=800Hzのsin波，y=1680Hzのノコギリ波**）でプロットすることによって出てくる図なんです．よく考えたら当たり前なのですが，ちょっとここまで綺麗な模様になるとは思いませんでしたので最初に出たときビックリしました．そして，これを混合行列Aによって適当に混ぜ合わせた波形を作り，（**x=混成波形1，y=混成波形2**）というようにすると――

　かなり歪な形になっちゃいます．実は，**ICAのやっていることを別の表現に言い換えると，この歪になったプロットグラフを元の形・元の軸に戻す**ということをやっています．で，**形をもとに戻す為にやっているのが白色化**（無相関化や球面化とも言われたりします）となってます．実際に，歪な形になった**２つの混成波形を白色化処理**にかけると――
　なんと，**形自体はすぐにもとに戻ってくれました！**　白色化については『具体例で学ぶ数学』様の[このページ](https://mathwords.net/musoukanka)も参照すると良いでしょう．
で，最後は軸をもとに戻すような回転行列を与えてやる必要がありますが，中心極限定理を元に考えると，尖度が3から最も遠くなるような軸がガウス分布から最も遠ざかった信号――つまり，元の信号となるはずですので，それに従って，**尖度が最大となる方向に回転を与えてやる**と……
　**見事に元の形に戻ってくれました！**　この時のx軸成分とy軸成分のデータをそれぞれ見てみても，やはり元の波形の形に戻っています！　
という訳で，以上がクソガバICA君の解説でした！

## ６．もっと気になった人向けに
　本当は「ICAにはどんな欠点があるの？」とか「ICAの応用例や拡張例には何があるの？」という疑問についての答えや，カクテルパーティ効果についての考察についてもこの記事内に書きたかったのですが，先日ついに終了した研究会に関する準備~~やFGO二部三章やグラブルのイベント周回~~もあってその為の学習時間が十分とれませんでしたので，代わりに参考になるサイトや書籍をいくつか紹介して締めたいと思います．

１．[独立成分分析入門　～音の分離を題材として～](http://www.kecl.ntt.co.jp/icl/signal/sawada/mypaper/subspace2010rev.pdf)
　NTT コミュニケーション科学基礎研究所の澤田 宏先生が執筆してくださった，ICAについての**非常に分かりやすい**解説スライドです．ざっと見た感じ，今回の記事で説明した部分はほぼほぼ網羅しているように感じられたので，分かり辛いところなどがあればそちらも参照すると良いと思います．

２．[詳解独立成分分析：信号解析の新しい時代](https://www.amazon.co.jp/%E8%A9%B3%E8%A7%A3-%E7%8B%AC%E7%AB%8B%E6%88%90%E5%88%86%E5%88%86%E6%9E%90%E2%80%95%E4%BF%A1%E5%8F%B7%E8%A7%A3%E6%9E%90%E3%81%AE%E6%96%B0%E3%81%97%E3%81%84%E4%B8%96%E7%95%8C-%E3%82%A2%E3%83%BC%E3%83%9D-%E3%83%93%E3%83%90%E3%83%AA%E3%83%8D%E3%83%B3/dp/4501538600)
　**独立成分分析を勉強する際の鉄板**となる書籍です．僕もこの本を見て独立成分分析を勉強しました！（ただし，最後まで読めてない）　基本的なICAの原理については勿論，ICAの基礎的分野（確率統計・数値解析）や発展的なICAについても凄く綿密な解説がされています．数式などを使った厳密な説明がされているので結構難解であったり，発行されたのが2005年と情報が古めだったりしているののが難点ですが，それを以って有り余るほどの情報量が詰め込まれているので，これから使う方は図書館にも置いてあるものを一度読んでみることを勧めます．（特に初学者の方は1,7,8章を最初に読むのがオススメ）

　カクテルパーティ効果についての文献についても探してはみたのですが，残念ながらパッと探してみた限りですとありませんでした．ですので，もし知っている方が居れば情報提供してくださると非常に助かります！
## ７．終わりに
　今回は自分の中で語れるものとして一番アツいものとして信号処理が最初に来ており，更にその中で門外漢の方が見ても面白そうなものとしてICA――独立成分分析を紹介しました．構成に迷う場面もありましたが，最終的に結構楽しく執筆ができ，こちらとしてもICAの良い復習となりました．一方で，取っ付きやすくするために題名にカクテルパーティ効果を付けておきながら，そっちに関する記述が不十分だったのは反省しております．
　また，今回の記事は数学が苦手な方の為に，自分ができる範囲内で最大限に詳しく，また厳密性を犠牲にして直感的に書いてみるようにしましたが，正直分かりやすいかというと自信がないので「ここが分からない」とか「ここが気になる」とかあれば形は問いませんので質問をしてくださると，僕としても非常に助かります．
　最後に，ここまで読んでくださった皆さん，ここまでお付き合いありがとうございます！　企画立案者のmofusit君もこの様な機会を設けてくださり，ありがとうございます！