■フーリエ解析（14）：　フーリエ変換を利用した音源位置の推定（CSP法）

　２次元平面上における音源の位置は、周囲に設置された３個のマイクに音が到達する時間差から推定できます。

　マイクの位置をM1、M2、M3とすると、音源位置はM1、M2を焦点とする双曲線1-2と M1、M3を焦点とする双曲線1-3の交点となります。

　２つのマイクM1とM2（あるいはM1とM3）に音が到達する時間差の算出法には幾つかの方法が知られていますが、ここではCSP法による到達時間差（Delay Of Arrival：DOA）算出法を利用して、音源位置を推定するアプレットを紹介します。

　CSP（Cross-power Spectrum Phase analysis：白色化相互相関法）では、２つのマイクの入力波形：s₁(n)、s₂(n)に対して、周波数領域でのマッチングを行い、類似度 CSP_1,2(h)が最大となる遅延量を求めてDOAとします。

具体的な計算手順は次のとおりです。

１．入力波形：s₁(n)、s₂(n) [n = 0, 1, .., N-1] を離散フーリエ変換（DFTまたはFFT）する。
　　　S₁(k) = ��_n=0^N-1 s₁ (n)W ^nk

　　　S₂(k) = ��_n=0^N-1 s₂ (n)W ^nk

　　　ここで、　W = e^-j2π/N
　　　　　　　　N　：　データ数（サンプル数）
　　　　　　　　k　：　0, 1, .., N-1

２．それらを振幅で正規化する。
　　　S₁(k) /｜S₁(k)｜ = ��_n=0^N-1 s₁ (n)W ^nk / ｜��_n=0^N-1 s₁ (n)W ^nk ｜

　　　S₂(k) /｜S₂(k)｜ = ��_n=0^N-1 s₂ (n)W ^nk / ｜��_n=0^N-1 s₂ (n)W ^nk ｜

３．両者の相互相関関数を求める。
　　　S_1,2(k) = S₁(k) S₂(k)^* /｜S₁(k)｜/｜S₂(k)｜

　　　ここで、　*　：　共役を表す

４．それを離散フーリエ逆変換（IDFTまたはIFFT）して、CSP係数を求める。
　　　CSP_1,2(h) = (1/N)��_k=0^N-1 S_1,2(k)W ^-hk

　　　ここで、　h　：　0, 1, .., N-1

５．CSP係数が最大となるサンプル差（位相差）hを求めて、δ_1,2とおく。
　　　δ_1,2 = argmax_h[CSP_1,2(h)]

６．到達時間差（DOA)を算出する。
　　　τ_1,2 = δ_1,2 / Fs

　　　ここで、　Fs　： サンプリング周波数（Hz）

・平面上（音場）にマイクが３個設置されています。
・ここで音源の位置をマウスでクリックします。
　　この点を通り、マイク位置１，２を焦点とする双曲線が青色で表示されます。
・各マイクで受音された音声データ（波形）が画面右上に表示されます。
・「計算」ボタンを押すと、この音声データからCSP法を利用して到達時間差、音源位置が算出されて画面左上に表示されます（ピンク）。
　　２本の双曲線の交点が算出（推定）された音源位置です。
・相互相関関数S_i,j(k)が画面左下に、CSP係数CSP_i,j(h)が画面右下に表示されます。



・マイク位置を変更することができます。
・「area表示」ボタンを押すと、音源位置の候補点が２つ（双曲線の交点が２個）ある領域が表示されます。
　　この領域で音が発生すると、３個のマイクではどちらが音源か判定できません。
・サンプリング周波数 Fs = 8200Hz、データ数 N = 1024としています。

　２つのマイクM1、M2に対して、音源がM1に近い場合はM2には遅れて音波が到達し、前述の到達時間差τs_1,2は正の値となります。
　逆に音源がM2に近い場合はM2には早く音波が到達しますが、前述の到達時間差τs_1,2はこの場合も正の値となります。

　音源が２つのマイクのいずれに近いかはδ_1,2の値の範囲で判断できます。
　　・δ_1,2 < N/2 の時　---> マイク１に近い
　　・δ_1,2 > N/2 の時　---> マイク２に近い

　これはサンプリング周波数 Fs（ = 8200Hz）、データ数 N（ = 1024）、音速 V（ = 340m/s）と現在の音場の広さ（マイク間の距離）の関係から導くことができます。

　位相差（サンプル差） N/2 = 512 は
　　　V x N/2/Fs = 340 x 512 / 8200 = 21.2m
の距離差に相当しますが、約20m四方の音場とマイク配置から最大距離差は約15mとなります。　従って、この条件下では距離差が20mを越えることはないわけです。

　また、CSP係数の次の性質：
　　　　CSP_2,1(N - h) = CSP_1,2(h)^*

　　　｜CSP_2,1(N - h)｜ = ｜CSP_1,2(h)｜

から、音源がマイク２に近い場合は（N - 本来の位相差）がδ_1,2として得られるので、この場合 本アプレットでは本来の位相差（N - δ_1,2）をマイナス表示にしています。

　到達時間差の算出法にはCSP法以外にも、入力波形 s₁(n)、s₂(n) [n = 0, 1, .., N-1]の相互相関関数による方法などもありますが、CSP法を用いた算出法が精度の点から好ましいと言われています。