(57)【要約】【課題】高残響下でも観測信号を有効に信号分離し、白色化歪み、スケーリング及び成分置換の問題を解決する耐高残響ブラインド信号分離装置の提供。【解決手段】周波数領域の観測信号に基づき、瞬時混合の周波数領域ＩＣＡにより成分置換のない分離行列の近似的な最適値（初期分離行列）を演算する。この初期分離行列を初期値として、最小歪み原理に基づく畳み込みＦＤＩＣＡによる周波数領域の分離行列の更新式を用いて分離行列の更新を行う。この分離行列を用いて観測信号を分離することで、高残響下でも信号を良好に分離することができる。【選択図】        図６
2006308955.tif 000002
【特許請求の範囲】
【請求項１】
複数の観測信号をフーリエ変換して周波数ごとの時系列データを算出し、当該時系列データから独立成分分析（以下「ＩＣＡ」という。）により周波数毎の分離行列を算出し、前記分離行列により分離信号を算出する耐高残響ブラインド信号分離装置であって、
周波数ごとの前記時系列データに対し、最小歪み原理に基づく畳み込み周波数領域ＩＣＡによる分離行列の更新式を用いて前記分離行列の更新を行う分離行列更新手段を有することを特徴とする耐高残響ブラインド信号分離装置。
【請求項２】
観測信号に基づき、各周波数間での行列要素の成分置換のない前記分離行列の近似的値（以下「初期分離行列」という。）を演算する分離行列初期化手段を備え、
前記分離行列更新手段は、前記初期分離行列を初期値として、前記分離行列の更新を行うことを特徴とする請求項１記載の耐高残響ブラインド信号分離装置。
【請求項３】
前記分離行列初期化手段は、前記観測信号に基づき、瞬時混合の周波数領域ＩＣＡにより成分置換のない前記初期分離行列を演算すること
を特徴とする請求項２記載の耐高残響ブラインド信号分離装置。
【請求項４】
複数の観測信号をフーリエ変換して周波数ごとの時系列データを算出し、当該時系列データから独立成分分析（以下「ＩＣＡ」という。）により周波数毎の分離行列を算出し、前記分離行列により分離信号を算出する耐高残響ブラインド信号分離方法であって、
周波数ごとの前記時系列データに対し、最小歪み原理に基づく畳み込み周波数領域ＩＣＡによる分離行列の更新式を用いて前記分離行列の更新を行う分離行列更新手順を有することを特徴とする耐高残響ブラインド信号分離方法。
【請求項５】
観測信号に基づき、各周波数間での行列要素の成分置換のない前記分離行列の近似的値（以下「初期分離行列」という。）を演算する分離行列初期化手順を備え、
前記初期分離行列を初期値として、前記分離行列更新手順において、前記分離行列の更新を行うことを特徴とする請求項４記載の耐高残響ブラインド信号分離方法。
【請求項６】
前記分離行列初期化手順において、前記観測信号に基づき、瞬時混合の周波数領域ＩＣＡにより成分置換のない前記初期分離行列を演算すること
を特徴とする請求項５記載の耐高残響ブラインド信号分離方法。
【請求項７】
コンピュータにより実行することにより、コンピュータを請求項１乃至３の何れか一に記載の耐高残響ブラインド信号分離装置として機能させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、統計的に独立な複数の原信号が線形混合した混合信号の信号分離を行う信号分離技術に関し、特に、独立成分分析（Independent Component Analysis：以下「ICA」という。）を用いて原信号の分離を行う技術に関する。
【背景技術】
【０００２】
混合信号を原信号に分離する信号処理技術はブラインド信号分離と呼ばれている。近年、ブラインド信号分離の手法としてICAが注目されており、音声信号処理、画像処理、通信等の分野において研究開発が進められている。ICAでは、例えば、複数の音源から発せられた原信号が線形混合した混合音をR個の受音器で受音して得られる複数の観測信号に対して、各原信号が統計的に独立であることを利用して観測信号を処理することで、受音器と同数のR個又はR以下の個数の原信号を分離することができる。
【０００３】
ICAは、時間領域で信号分離を行う時間領域ICA（Time Domain ICA：以下「TDICA」という。）と周波数領域で信号分離を行う周波数領域ICA（Frequency Domain ICA：以下「FDICA」という。）とに大別される。
【０００４】
2つ以上の受音器に到達する音波に時間遅れや畳み込みがなく、混合行列の要素が実定数で表される場合を「空間的混合（spatial mixture）」または「瞬時混合（instantaneous mixture）」という。瞬時混合の場合、TDICA及びFDICAの何れによっても音源を完全に分離することができる。
【０００５】
しかし、実際の環境においては、音源と受音器との距離の相違から、個々の受音器において受信される信号に時間差がある。また、反響や残響が畳み込まれて観測される。このように、空間的・時間的に原信号が混合される場合を「時空間的混合（spatio-temporal mixture）」または「畳み込み混合（convolutive mixture）」という。また、畳み込み混合の原信号の分離問題は「ブラインド・ソース・デコンボリューション（Blind Source Deconvolution：以下「BSD」という。）」と呼ばれる。ブラインド・ソース・デコンボリューションの問題は次のように定式化される。
【０００６】
〔１〕ブラインド・ソース・デコンボリューション（BSD）
音源の数をR1、受音器の数をR2とする。音源から出力される原信号をsi（t）、受音器における観測信号をxj（t）と記す。iは音源の番号を表すインデックスであり、jは受音器の番号を表すインデックスである。
【０００７】
時間領域での観測信号xj（t）を、L次のフィルタで次式のようにモデル化する。
【０００８】
【数１】
2006308955.tif 000003
【０００９】
ここで、tは時間、τは遅延時間を表す。また、各音源から各受音器までのインパルス応答をaji（τ）と記した。行列A（τ）＝（aji（τ））を「混合行列」という。一般的なインパルス応答aji（τ）の形状は、適当な時間経過後にパルス的な強い応答を持ち、時間と共に減衰していくものである。このインパルス応答長をTとする。但し、Tは未知である。
【００１０】
BSDの目的は、原信号si（t）やインパルス応答aji（τ）の情報なくして、分離のための長さLのFIR（Finite Impulse Response）フィルタの係数wij（τ）と分離信号yi（t）とを求めることにある。尚、係数の行列W（τ）＝（wij（τ））は「分離行列」という。
【００１１】
【数２】
2006308955.tif 000004
【００１２】
図12はR1＝2、R2＝2の場合のBSDの構成を表す。式（1），式（2）を周波数領域に変換すると次式のようになる。
【００１３】
【数３】
2006308955.tif 000005
【００１４】
ここで、xj（ωn，k），si（ωn，k）はそれぞれxj（t），si（t）の離散フーリエ変換（Discrete Fourier Transform：以下「DFT」という。）、aji（ωn），wij（ωn）はそれぞれaji（t），wij（t）のDFT、ωｎは正規化角周波数、kはフレーム番号を表す。
【００１５】
TDICAでは時間領域において式（2）のwij（τ）を計算することによって信号分離を行い、FDICAでは周波数領域において式（4）のwij（ωn）を計算することによって信号分離を行う。BSDにおいて、TDICAにより混合信号から元の原信号を分離しようとすると、次数（タップ長）Lが数千にも及ぶフィルタが必要とされ、アルゴリズムの安定性や計算負荷の問題が生じる。従って、一般にBSDはFDICAにより解かれることが多い。これは、時間領域での畳み込みが周波数領域では積項として表現できるため、瞬時混合の場合と同様に平易に扱うことができることによる。
【００１６】
〔２〕TDICA
次に、TDICAについての概略を簡単に説明する。TDICAでは、分離信号Y（t）＝［y1（t），…，yR1（t）］T（但し、xTはxの転置を表す。）の各成分が独立となるように分離行列W（τ）＝（wij（τ））を決定する。独立性を測る距離としては、分離信号Y（t）の各成分の同時分布と周辺分布との間のKullback-Leibler情報量（以下「KL情報量」という。）IKL（Y（t））を用いることができる。
【００１７】
【数４】
2006308955.tif 000006
【００１８】
但し、C（Y）は、Yのエントロピーを表す。分離信号Y（t）の各成分を完全に独立にするためには、IKL（Y（t））＝0となるように分離行列W（τ）を決定すればよい。KL情報量は非負なのでIKL（Y（t））のW（τ）による微分が0のとき最小となる。これを解くと、分離行列の更新式が次式のように定式化される。
【００１９】
【数５】
2006308955.tif 000007
【００２０】
ここで、IはR1＝R2のときR1×R1の単位行列、R1＜R2のときR1×R2の(I:O)なる行列、OはR1×(R2-R1)の零行列であり、ηは探索ステップ幅、δ（τ）はクロネッカーのデルタである。また、φは理想的にはyiの周辺確率密度pyi（y）を用いて次式で定義される関数である。しかし、一般には音源の正確な分布は分からないため、φはtanh(・)のような非線形関数を用いる。
【００２１】
【数６】
2006308955.tif 000008
【００２２】
ここで、まず、TDICAアルゴリズム（式（6））の非対角成分について考える。非対角成分は、次式（10）のように表される。
【００２３】
【数７】
2006308955.tif 000009
【００２４】
従って、W（τ）の収束点では、次式（11）が成り立つ。
【００２５】
【数８】
2006308955.tif 000010
【００２６】
このことは、収束後の分離信号yi（t）とyj（t）とが空間的にだけでなく時間的にも独立となることを示している。すなわち、式（6）の非対角成分は分離信号が空間的にも時間的にも独立となるように機能する。
【００２７】
一方、式（6）の対角成分について考える。対角成分は次式（12）のように表される。
【００２８】
【数９】
2006308955.tif 000011
【００２９】
従って、W（τ）の収束点では、次式（13）が成り立つ。
【００３０】
【数１０】
2006308955.tif 000012
【００３１】
このことは、分離信号yi（t）がφ（yi（t））yi（t）＝1にスケーリングされることに加えて、時間的にも独立となることを示している。すなわち、式（6）の対角成分は、分離信号をスケーリングするとともに時間的に独立化（白色化）するように働く。時間的な独立化については、非対角成分の場合、もともと異なる信号に対して独立化するように働くため問題とはならない。しかし、対角成分の場合、1つの信号を時間的に独立化することとなるため、歪みを生じる。この歪みを「白色化歪み」という。
【００３２】
図13は、TDICAからみた混合過程と分離過程を表す図である。白色信号si’（t）の分布がi.i.d（identically independent distribution）の音源と仮定し、原信号si（t）はこれを口腔内の調音器官で調音された音声信号である。原信号si（t）は白色信号si’（t）と異なり時間的な相関を持つ信号となる。従って、TDICAアルゴリズムは、分離対象の信号をi.i.dを仮定して、その過程を満たすように分離行列の更新を行う。その結果得られる分離信号yi（t）は、上述のように空間的に独立化されるのみならず時間的にも独立化（白色化）されたものとなる。そのため、音声のように時間的な相関を持つ原信号si（t）にTDICAアルゴリズムを適用した場合、分離信号yi（t）はsi（t）ではなくsi’（t）の推定値となり歪んだものとなる。この歪みは、白色雑音の歪みと似ていることから、白色化歪みと呼ばれる。
【００３３】
〔３〕最小歪みの原理
TDICAにおける白色化歪みを解決する方法としては、非特許文献3に記載の最小歪みの原理（Minimal Distortion Principle：以下「MDP」という。）による方法が公知である。以下、MDPの概略と、TDICAへの適用について簡単に説明する。
【００３４】
まず、si’（t）からxj（t）までの混合過程とxj（t）からyi（t）までの分離過程を図13のようにブロック化して考える。TDICAにより求められる分離フィルタ行列をz変換の形式で式（14）のように表す。
【００３５】
【数１１】
2006308955.tif 000013
【００３６】
TDICAにより音源si（t）を分離しようとしても、分離フィルタ行列W（z）は、混合行列の逆行列A-1（z）とはならず、式（15）のように、白色化による対角行列D（z）がA-1（z）に掛かる形で与えられる。
【００３７】
【数１２】
2006308955.tif 000014
【００３８】
ここで、D（z）はフィルタdi（z）（i＝1，…，R1）を対角要素とする対角行列である。一方、音源s（t）は白色信号s’（t）を対角成分がフィルタei（z）（i＝1，…，R1）の対角行列E（z）で変調したものと考えられる。従って、TDICAではD（z）としてE-1（z）を推定しているとも解釈することができる。
【００３９】
ここで、分離行列W（z）に対して、次式（16）のような操作を行うことにより、D（z）を相殺することができる。
【００４０】
【数１３】
2006308955.tif 000015
【００４１】
従って、式（16）の左辺を改めてW（z）とおくことによって、W（z）からD（z）の任意性を取り除くことができる。
【００４２】
【数１４】
2006308955.tif 000016
【００４３】
式（17）は式（16）においてD（z）＝diag A（z）としたものと等価である。このとき、分離信号yi（t）はaii（z）si（t）となる。従って、yi（t）のスケーリングは伝達関数で規定され理にかなったものとなり、分離信号の歪みは最小となる。
【００４４】
D（z）がdiag A（z）になるようにするには、歪みをE［‖y（t）−x（t）‖2］で評価し、この評価関数を最小とするW（z）を求めればよいことが示されている（非特許文献3参照）。これが最小歪みの原理（MDP）である。
【００４５】
次に、MDPに基づくTDICAのアルゴリズムについて説明する。MDPに従ってE［‖y（t）−x（t）‖2］を最小化するために、E［‖y（t）−x（t）‖2］の自然勾配を
【００４６】
【数１５】
2006308955.tif 000017
とすると、
【００４７】
【数１６】
2006308955.tif 000018
は式（18）で与えられる。また、そのon-line展開形は式（19）で与えられる。
【００４８】
【数１７】
2006308955.tif 000019
【００４９】
一方、TDICAの分離行列の更新式（式（14））を非対角成分ΔWoff（τ）と対角成分ΔWdiag（τ）とに分けて式（20）〜（22）のように書くことができる。
【００５０】
【数１８】
2006308955.tif 000020
【００５１】
ここで、ξはΔWoff（τ）に対するΔWdiag（τ）の比率を表し、式（6）ではξ＝1である。上述のように、非対角成分は白色化には寄与しないのでそのまま採用する。ΔWoff（τ）をon-line展開形で表すと式（23）のようになる。
【００５２】
【数１９】
2006308955.tif 000021
【００５３】
一方、対角成分ΔWdiag（τ）はスケーリング調整に寄与し白色化の原因ともなるので、MDPに基づいて導かれる
【００５４】
【数２０】
2006308955.tif 000022
と置換する。以上の結果、MDPに基づくTDICA法における分離行列の更新式のon-line展開形は、式（24）により表される。
【００５５】
【数２１】
2006308955.tif 000023
【００５６】
ここで、ψ（t）＝y（t）−x（t）、γ＝2ξである。γは非対角項に対する対角項の更新比率を表す。
【００５７】
〔４〕FDICA
次に、FDICAについての概略を簡単に説明する。FDICAでは、分離信号Y（ωn，k）＝［y1（ωn，k），…，yR1（ωn，k）］T（但し、xTはxの転置を表す。）の各成分が独立となるように分離行列W（ωn）＝（wij（ωn））を決定する。分離行列W（ωn）の更新式はTDICAの場合と同様にして求められ、次式のようになる。
【００５８】
【数２２】
2006308955.tif 000024
【００５９】
ここで、Iは式（６）に対してと同様に定義される単位行列、ηは探索ステップ幅であり、φは式（8）（9）と同様に定義される非線形関数である。
【００６０】
FDICAにおいては、成分置換（permutation）及びスケーリング（scaling）が大きな問題となる。すなわち、TDICAの場合と同様、FDICAにより得られる分離信号の周波数毎の時系列Y（ωn，k）＝［y1（ωn，k），…，yR1（ωn，k）］Tは、要素間の順序の任意性とスケーリング（大きさと位相）の任意性がある。各要素の順序が入れ替わったとしても、また、各要素のスケーリングが変化しても、各要素間の独立性は保たれるからである。式（4）より、分離信号Y（ωn，k）に関するこれらの任意性は、分離行列W（ωn）＝（wij（ωn））の任意性に対応している。従って、分離行列W（ωn）の各行成分が同じ信号に対応するように並べる操作と、各行成分のスケーリングを適切に調節する操作が必要とされる。前者が成分置換の問題であり、後者がスケーリングの問題である。
【００６１】
図14では、音源の数及び受音器の数が共に2である場合の例を示している。FDICAでは観測信号x1（t），x2（t）は時間幅Nのフレームごとにフーリエ変換されて、N個のサブバンド｛x1（ω0，k），…，x1（ωN-1，k）｝，｛x2（ω0，k），…，x2（ωN-1，k）｝に分解される。そして、各サブバンドについて、瞬時混合のFDICAが実行され、N組の分離信号｛y1（ωn，k），y2（ωn，k）｜n＝0，…，N−1｝が計算される。そして、これらの分離信号は逆フーリエ変換され、時間領域の分離信号y1（t），y2（t）として出力される。ここで、図14の例では、成分置換により分離信号y1（ωn，k）と分離信号y2（ωn，k）の成分がωn＝ω0のときとωn＝ω1のときで入れ替わっている。このような成分置換が生じると、逆フーリエ変換により時間領域の分離信号y1（t），y2（t）を再生した場合にサブバンド成分が再び混合して歪みが生じる。
【００６２】
同様に、スケーリングの任意性により、サブバンドごとに分離信号y1（ωn，k），y2（ωn，k）のゲインが異なる。従って、逆フーリエ変換により時間領域の分離信号y1（t），y2（t）を再現した際に、不自然な重み付けがなされて歪みが生じる。
【００６３】
スケーリングの問題を解決する方法としては、例えば、上述の非特許文献3に記載のMDPによる方法が公知である（特許文献4，〔0011〕参照）。
【００６４】
また、成分置換の問題を解決する方法としては、例えば、特許文献1〜3及び非特許文献1，2，5〜7等が公知である。ここでは、一例として、非特許文献1に記載の成分置換解消法について説明しておく。簡単のため、音源と受音器の数がともに2である場合について説明する。
【００６５】
まず、FDICAにより分離行列W（ωn）が求められたとし、分離信号yi（ωn,k）は観測信号xj（ωn,k）により式（4）で表されるとする。ここで、分離信号yi（ωn,k）に対して分割スペクトルvi（ωn，k）を次式のように定義する。
【００６６】
【数２３】
2006308955.tif 000025
【００６７】
この分割スペクトルは、成分置換がないときは、式（30）となり、成分置換が起きているときには式（31）となり、音源si（ωn，k）と音源−受音器間の伝達関数aji（ωn）の積で一意的に表される。
【００６８】
【数２４】
2006308955.tif 000026
【００６９】
【数２５】
2006308955.tif 000027
【００７０】
また、vi（ωn，k）の第1成分と第2成分は同じ音源si（ωn，k）に対する推定値となっているため、二つの音源s1（ωn，k），s2（ωn，k）の推定値は、v11（ωn，k）かv22（ωn，k）の何れかの分割スペクトルで代表させることができる。
【００７１】
この成分置換解消法では、音声に対する分解スペクトルのエントロピーは、雑音に対する分解スペクトルのエントロピーよりも高いことを利用する。
【００７２】
そこで、分解スペクトルのエントロピーを導入するために、まず、分割スペクトルの実部Re｛vij（ωn）｝の分布範囲をM等分したときの各区画をιm（m＝0，…，Γ−1）とし、区画ιmに入る頻度をκij（ωn，ιm）とする。また、κij（ωn，ιm）を規格化して得られる確率ρij（ωn，ιm）を下式で定義する。
【００７３】
【数２６】
2006308955.tif 000028
【００７４】
この確率ρij（ωn，ιm）を用いて、分割スペクトルの実部Re｛vij（ωn）｝の分布に対するエントロピーを次式で定義する。
【００７５】
【数２７】
2006308955.tif 000029
【００７６】
分割スペクトルvij（ωn，k）は音源si（ωn，k）と音源−受音器間の伝達関数aji（ωn）の積で一意的に表されるので、Ci1とCi2とは同じ音源に対するエントロピーを表している。そこで、両者の差（式（34））を検査して、エントロピー差ΔC（ωn）が負ならば成分置換はないと判定し、最終的な出力をz（ωn，k）＝［v11（ωn，k），v22（ωn，k）］として割り当てる、逆に、差ΔC（ωn）が正であれば、成分置換が起きていると判断し、最終的な出力をz（ωn，k）＝［v22（ωn，k），v11（ωn，k）］として割り当てる。
【００７７】
【数２８】
2006308955.tif 000030
【００７８】
尚、上記説明では、音源及び受音器の数がともに2の場合で説明したが、音源数R1、受音器数R2の場合に拡張することは容易である。以上のようにして、FDICAにおける成分置換の問題は解決することができる。
【００７９】
〔５〕畳み込みFDICA（ConvFDICA）
FDICAにおいては、残響が多くなると音源分解性能が著しく低下するという問題がある。これは、信号をDFTにより周波数領域に変換する際に、フーリエ変換の分析窓幅が狭いとインパルス応答のすべてを包含できないため、音源の直接音部分は分離できても、それに続く反射や残響部分が分離できないことに起因する。図15にこの様子を示す。図15は、インパルス応答長Tに対して窓幅Nの分析窓を使う場合を表している。この場合、Nよりも大きなインパルス応答はDFTに取り込まれない。従って、この取り込まれない残響部分が分離できないことになる。
【００８０】
そこで、分析窓幅を広くすれば残響成分も除去できるように思われる。しかしながら、同一長の観測信号について、分析窓幅を広くすると、時間分解能が粗くなり、結果的にICAに必要な十分な個数のデータが得られない。従って、分析窓幅を小さくしたとしてもICAの分解能力は低下する。また、観測信号が長ければ、分析窓幅を広くすることにより残響成分を除去することができるが、実用を考えた場合、長時間に亘って同じ環境が持続するとは考えにくい。また、移動音源に対しては短時間の観測信号で分離する必要がある。
【００８１】
この問題を解決するために、非特許文献4においては、分割畳み込みの概念を導入し、混合過程のインパルス応答を分割することにより、それまで周波数領域で瞬時混合として近似的にモデル化していたFDICAを畳み込み混合として厳密にモデル化している。そして、時間領域畳み込み混合に対するTDICAアルゴリズムを周波数領域に適用した、畳み込みFDICA（以下「ConvFDICA」という。）法が提案されている。
【００８２】
ConvFDICAでは、狭帯域の分離フィルタ行列を式（35）のように定義する。このW（ωn，z）を「狭帯域分離フィルタ行列」と呼ぶ。
【００８３】
【数２９】
2006308955.tif 000031
【００８４】
これを周波数領域での畳み込みモデルに対してTDICAの考え方を適用すると、次式（36）のような更新行列が得られる。これを、ConvFDICAと呼ぶ。
【００８５】
【数３０】
2006308955.tif 000032
【００８６】
ここで、kはフレーム番号、Mは狭帯域分離フィルタのタップ長、L=MNである。また、ΔW（ωn，p），W（ωn，p），φ（y（ωn，p）），y（ωn，p）はそれぞれ次式のように表される。
【００８７】
【数３１】
2006308955.tif 000033
【００８８】
ConvFDICAでは、分析窓幅をN、各周波数における分離フィルタのタップ長をMとして、長さMNの時間領域フィルタを周波数領域でMNと各周波数に分解する。これにより、短時間の観測信号からでも十分な統計量が得られる程度の適度な分析窓幅で高残響下でも有効に機能し収束性に優れた分離フィルタが実現される。
【００８９】
しかしながら、ConvFDICAは、TDICAの性質を継承しているため、TDICA固有の白色化歪みの問題を伴う。更に、周波数領域で分離を行うことに起因してスケーリングの不定性や成分置換の問題も生じてくる。
【特許文献１】特開2005−91560号公報
【特許文献２】特開2005−91732号公報
【特許文献３】特開2004−302122号公報
【特許文献４】特開2005−79781号公報
【非特許文献１】金田圭市，古屋武志，五反田博，「分割スペクトルのエントロピーに基づく成分置換解消法」，電子情報通信学会論文誌，Vol.J87-A, No.7, pp. 1065-1069, July 2004.
【非特許文献２】H. Sawada, R. Mukai, S. Araki, S. Makino, "A Robust and Precise Method for Solving the Permutation Problem of Frequency-Domain Blind Source Separation," IEEE Trans. Speech and Audio Processing, vol.12, no. 5, pp. 530-538, Sep. 2004.
【非特許文献３】K. Matsuoka and S. Nakashima, "Minimal Distortion Principle for Blind Source Separation". Proc. ICA 2001, pp.722-727, 2001.
【非特許文献４】C. Serviere, "Separation of speech signals with segmentation of the impulse responses under reverberant conditions   ICA, Int.Workshop on Independent Component Analysis and Signal Separation, Nara, Avril 2003
【非特許文献５】猿渡洋，栗田悟史，武田一哉，板倉文忠，鹿野清宏，「帯域分割型ICAとビームフォーミングを利用したブラインド音源分離」, 信学技報，EA2000-14，pp.1-8，2000.
【非特許文献６】石橋孝昭，金田圭市，古屋武志，五反田博，「分割スペクトルの位相に基づく音源方向の推定」, 第56回電気関係学会九州支部連合大会講演会，講演論文集(CD-ROM), 12-2A-10, (2003, H15-9, 崇城大学)
【非特許文献７】H. Gotanda, K. Nobu, T. Koya, K. Kaneda, T. Ishibashi, N. Haratani, "Permutation correction and speech extraction based on split spectrum through FastICA," Proc. ICA2003, pp.379-384, 2003.
【非特許文献８】S. Araki, S. Makino, R. Aichner, T. Nishikawa, and H. Saruwatari, "Subband based blind source separation with appropriate processing for each frequency band, " Proc. ICA2003, pp.499-504, 2003.
【発明の開示】
【発明が解決しようとする課題】
【００９０】
以上、従来の信号分離技術の問題点をまとめると次の通りである。
【００９１】
（1）TDICAでは、成分置換の問題は応用上重大な支障とはならないが、スケーリングの不定性及び白色化歪みの問題は重大な支障を及ぼす。これらの問題はMDPを適用したTDICAによって解決することができる。
【００９２】
しかしながら、反射や残響を考慮したBSDモデルにTDICAを適用しようとすると、フィルタの次数が莫大となるため実装上の問題が大きい。また、アルゴリズムの問題や計算負荷の問題も大きい。従って、BSDモデルへの適用は不適であり、残響の多い環境においてはTDICAは有効ではない。
【００９３】
（2）FDICAは、サブバンド毎に瞬時混合によって分離を行うため、BSDモデルに適用した場合でも小規模な演算回路で実現することができ、収束も速く解の安定性にも優れる。一方、FDICAでは、成分置換の問題、スケーリングの問題、及びインパルス応答が長い場合に残響を完全に分離することができないという問題がある。
【００９４】
成分置換の問題に関しては、例えば、上述の分解スペクトルのエントロピーによる成分置換解消法などによって解決されている。また、スケーリングの問題に関しては、上述のMDPに基づくTDICA法により解決されている。
【００９５】
しかしながら、観測信号のインパルス応答が長い場合に残響を完全に分離することができないという問題に関しては解決できない。
【００９６】
（3）非特許文献4に記載のConvFDICAによれば、短時間の観測信号からでも十分な統計量が得られる程度の適度な分析窓幅で、高残響下でも有効に機能し収束性にも優れた分離フィルタを実現することができる。しかしながら、ConvFDICAは、TDICAの性質を継承しているため、TDICAに固有の白色化歪みの問題を伴う。さらには、周波数領域で必然的に起こるスケーリングの不定性や成分置換の課題も残る。
【００９７】
そこで、本発明の目的は、高残響下でも観測信号を有効に分離することができ、白色化歪みの問題、スケーリングの問題及び成分置換の問題をすべて解決することが可能な信号分離技術を提供することにある。
【課題を解決するための手段】
【００９８】
本発明に係る耐高残響ブラインド信号分離装置の第1の構成は、複数の観測信号をフーリエ変換して周波数ごとの時系列データを算出し、当該時系列データから独立成分分析（以下「ICA」という。）により周波数毎の分離行列を算出し、前記分離行列により分離信号を算出する耐高残響ブラインド信号分離装置であって、周波数ごとの前記時系列データに対し、MDPに基づくConvFDICAによる分離行列の更新式を用いて前記分離行列の更新を行う分離行列更新手段を有することを特徴とする。
【００９９】
この構成によれば、MDPを用いていることにより、白色化歪みの問題を解決することができる。また、MDPを用いることにより、スケーリングの問題も解決される。さらに、周波数領域での分離行列を更新式するため、TDICAに比べて収束性が速い。また、畳み込みモデルを用いているため、残響を除去することが可能となる。
【０１００】
ここで、「MDPに基づくConvFDICAによる分離行列の更新式」とは、後述の式（85）によって与えられる更新式をいう。
【０１０１】
本発明に係る耐高残響ブラインド信号分離装置の第2の構成は、前記第1の構成において、観測信号に基づき、各周波数間での行列要素の成分置換のない前記分離行列の近似的値（以下「初期分離行列」という。）を演算する分離行列初期化手段を備え、前記分離行列更新手段は、前記初期分離行列を初期値として、前記分離行列の更新を行うことを特徴とする。
【０１０２】
この構成によれば、最初に、各周波数間での行列要素の成分置換のない初期分離行列を演算し、この初期分離行列を初期値とすることで、分離行列の最終的な収束点の近傍に初期値を設定することができる。そして、この初期値を用いて分離行列の更新を行うことにより、成分置換を生じることなく分離行列を最終的な収束点に収束させることができる。
【０１０３】
すなわち、分離係数の解空間は、模式的に示すと図3のような複雑な曲面であると考えられる。この場合、正しい収束点が独立性の尺度が最小となるが、それ以外にも成分置換を伴う解が極小値として存在する。一方、ICAは、自然勾配法によって分離係数空間での解の探索を行うため、正しい収束点の近傍に初期値を設定すれば、曲率が小さくなる方向に移動して、成分置換の収束点に飛び移ることなく正しい収束点に収束する。本発明は、この性質を利用することによって、成分置換の問題を解決したものである。
【０１０４】
ここで、「各周波数間での行列要素の成分置換のない初期分離行列を演算」する方法は、特に限定するものではなく、上述した各種公知の方法を用いることができる。
【０１０５】
本発明に係る耐高残響ブラインド信号分離装置の第3の構成は、前記第2の構成において、前記分離行列初期化手段は、前記観測信号に基づき、瞬時混合の周波数領域ICAにより成分置換のない前記初期分離行列を演算することを特徴とする。
【０１０６】
本発明に係る耐高残響ブラインド信号分離方法の第1の構成は、複数の観測信号をフーリエ変換して周波数ごとの時系列データを算出し、当該時系列データからICAにより周波数毎の分離行列を算出し、前記分離行列により分離信号を算出する耐高残響ブラインド信号分離方法であって、周波数ごとの前記時系列データに対し、MDPに基づくConvFDICAによる分離行列の更新式を用いて前記分離行列の更新を行う分離行列更新手順を有することを特徴とする。
【０１０７】
本発明に係る耐高残響ブラインド信号分離方法の第2の構成は、前記第1の構成において、観測信号に基づき、各周波数間での行列要素の成分置換のない前記分離行列の近似的値である初期分離行列を演算する分離行列初期化手順を備え、前記初期分離行列を初期値として、前記分離行列更新手順において、前記分離行列の更新を行うことを特徴とする。
【０１０８】
本発明に係る耐高残響ブラインド信号分離方法の第3の構成は、前記第2の構成において、前記分離行列初期化手順において、前記観測信号に基づき、瞬時混合の周波数領域ICAにより成分置換のない前記分離行列の近似的な最適値である初期分離行列を演算することを特徴とする。
【０１０９】
本発明に係るプログラムは、コンピュータにより実行することにより、コンピュータを前記第1乃至3の何れか一の耐高残響ブラインド信号分離装置として機能させることを特徴とする。
【発明の効果】
【０１１０】
以上のように、本発明によれば、高残響下でも観測信号を有効に分離することができ、白色化歪みの問題、スケーリングの問題及び成分置換の問題をすべて解決することが可能な信号分離技術を提供することが可能となる。
【発明を実施するための最良の形態】
【０１１１】
以下、本発明を実施するための最良の形態について、図面を参照しながら説明する。最初に、基本的な考え方についての説明を行った後に、本発明の実施例について説明する。
【０１１２】
〔１〕混合過程の厳密なモデル化
本発明の耐高残響ブラインド信号分離装置では、時間周波数領域においてICAを行う。ここで、時間周波数領域におけるICAでは、式（1）で表される観測信号をDFTによりサブバンド化し、フレーム単位で時系列化されたサブバンドの観測信号を用いてICAが行われる。しかしながら、式（1）を単純にDFTにより変換した式（3）を、再度時間領域に逆変換すると、その時間領域表現は巡回畳み込みとなる。このとき、巡回畳み込みはL＞＞Tの場合には式（1）の線形畳み込みと等しいが、L＞＞Tでないときは異なる。このことに起因して誤差を生じる。
【０１１３】
そこで、まず巡回畳み込みに起因する誤差の解消法について説明し、その後、式（1）の混合過程を周波数領域で厳密に反映させたモデルについて説明する。
【０１１４】
まず、式（1）の線形畳み込みを行列表現すると、次式（41）のように表される。
【０１１５】
【数３２】
2006308955.tif 000034
一方、式（3）は逆DFTにより時間領域に戻すと、次式（45）のように表される。
【０１１６】
【数３３】
2006308955.tif 000035
【０１１７】
ここで、式（41）と式（45）との比較から、線形畳み込みの行列表現と巡回畳み込みの行列表現とは異なることが分かる。そこで、Aji（L）の後ろにL個の0を付加したものをAji（2L）とし、長さを2Lとして式（3）の線形畳み込み表現を書き換えると、次式（49）のようになる。
【０１１８】
【数３４】
2006308955.tif 000036
【０１１９】
また、巡回畳み込み表現についても同様に書き換えると、次式（53）のようになる。
【０１２０】
【数３５】
2006308955.tif 000037
【０１２１】
ここで、式（49）と式（53）との比較から、
【０１２２】
【数３６】
2006308955.tif 000038
の後半のL個の
【０１２３】
【数３７】
2006308955.tif 000039
から
【０１２４】
【数３８】
2006308955.tif 000040
と、Xj（2L）（t）の後半のL個とは同じであるが、前半のL個は異なることが分かる。
【０１２５】
巡回畳み込みによる
【０１２６】
【数３９】
2006308955.tif 000041
と線形畳み込みによるXj（2L）（t）の間には上述のような相違があり、これに起因して誤差が生じる。しかしながら、両者は以下のような手順を踏むことによって等しくすることができる（非特許文献4参照）。
【０１２７】
まず、最初のL個の窓を式（57）により定義して、このhの各要素を対角成分とする対角行列Hを定義する。
【０１２８】
【数４０】
2006308955.tif 000042
【０１２９】
対角行列Hを式（53）の両辺に掛けることにより、次式（58）のようになり、式（59）が成立し、線形畳み込みと巡回畳み込みとは等しくなる。
【０１３０】
【数４１】
2006308955.tif 000043
【０１３１】
【数４２】
2006308955.tif 000044
【０１３２】
一方、フーリエ変換の作用素である行列Fを式（60）により定義する。
【０１３３】
【数４３】
2006308955.tif 000045
【０１３４】
これを式（58）の両辺に掛けると、式（62）となる。
【０１３５】
【数４４】
2006308955.tif 000046
【０１３６】
更に、これを変形することにより、最終的に式（63）が得られる。これは、式（58）の時間領域表現をDFTして周波数領域表現にしたものである。そこで、この式（63）を時間領域畳み込み混合過程の周波数領域における厳密なモデルと呼ぶ。
【０１３７】
【数４５】
2006308955.tif 000047
【０１３８】
〔２〕インパルス応答の分割による混合モデル
高残響下では、インパルス応答の時間Tが長いため、残響を十分に除去する為には式（1）の畳み込みのフィルタ長LをTに対して大きく採る必要がある。しかしながら、上述したように、Lを大きくすると、各サブバンドで音源を推定するのに十分な個数のデータが得られなくなり収束性に問題が生じる。一方、フィルタ長LをL＜Tとすると、インパルス応答を正確に反映したDFTとはならず、窓がインパルス応答全体を包含できないことに起因する誤差が生じる。
【０１３９】
そこで、まず適切な窓幅Nを仮定して、フィルタ長LをL＝MNとなるようにN等分する。LはTよりも十分に大きくなるように設定される。このとき、式（1）の線形畳み込みは、次式（69）のようになる。
【０１４０】
【数４６】
2006308955.tif 000048
【０１４１】
さらに、これを行列表現すると、次式（70）のように表すことができる。
【０１４２】
【数４７】
2006308955.tif 000049
【０１４３】
式（70）をそのまま周波数領域に変換しても正確に反映されないので、式（58）の関係を利用して、次式（74）のように変形する。
【０１４４】
【数４８】
2006308955.tif 000050
【０１４５】
更に、式（74）の両辺にDFTの変換行列Fを掛けて周波数領域に変換することにより、最終的に厳密なモデルが式（78）により得られる。
【０１４６】
【数４９】
2006308955.tif 000051
【０１４７】
ここで、
【０１４８】
【数５０】
2006308955.tif 000052
は式（79）を対角成分とする対角行列であり、また、それぞれの記号は、式（79）〜（83）により表される。
【０１４９】
【数５１】
2006308955.tif 000053
【０１５０】
式（78）の行列表現を各サブバンドωn（n＝0，1，…，2N−1）に分けた混合過程として記述すれば、式（84）のように、周波数領域での畳み込み表現となる。これを離散窓掛け離散フーリエ変換と呼ぶ。
【０１５１】
【数５２】
2006308955.tif 000054
【０１５２】
これを分かりやすく図で表すと図1のようになる。図1において、観測信号xj（t）は、各時間フレームkにおいて窓幅NでDFTがされるとともに、各サブバンドはフレーム0からフレームKの間で畳み込み演算がされる。これにより、すべてのインパルス応答をすべて取り込むことができると共に、適度な窓幅Nとして計算するために、各サブバンドで音源を推定するのに十分な個数のデータを短時間で得ることが可能となる。
【０１５３】
〔３〕本発明のICAアルゴリズム
本発明では、以上説明したような離散窓掛け離散フーリエ変換を用いることによって、式（1）の観測信号をサブバンド分解し、各サブバンドに対してMDPに基づくConvFDICAアルゴリズムを行う。この場合、分離行列の更新式は、下式（85）により表される。
【０１５４】
【数５３】
2006308955.tif 000055
【０１５５】
ここで、ηは探索ステップ幅、γは内分比である。また、ΔW（ωn，p）、W（ωn，p）、φ（y（ωn，k））、ψ（ωn，k），及びy（ωn，k）は、それぞれ式（86）〜（90）により表される。尚、音源の数と受音器の数はともにRとする。
【０１５６】
【数５４】
2006308955.tif 000056
【０１５７】
また、φi（ｙ（ωn，ｋ））、ψi（ωn，ｋ）は、それぞれ式（９１），（９２）により表される。
【０１５８】
【数５５】
2006308955.tif 000057
【０１５９】
このとき、ｙ（ωn，ｋ）は式（９３）のように計算される。
【０１６０】
【数５６】
2006308955.tif 000058
【０１６１】
本発明のアルゴリズムは、上述のServiereによるConvFDICAのアルゴリズム（非特許文献4参照）にMDPの概念を適用したものとなっている。以下、この本発明のアルゴリズムを「MDPに基づくConvFDICA（MDP-ConvFDICA）」と呼ぶ。このアルゴリズムを用いてICAを実行することにより、ConvFDICAにおいて問題となっていた白色化による歪みを解消することが可能となると共に、スケーリングの不定性の問題も解消することができる。
【０１６２】
〔４〕成分置換の解消法
最後に、成分置換の解消法について説明する。ICAでは、出力信号が互いに独立となるように分離フィルタの係数を更新し真の分離フィルタを求める。分離フィルタの係数をパラメータ、真の分離フィルタを収束値とすると、独立性の尺度（評価関数）が最小となるようにパラメータ空間を探索し、収束値を求める。従って、パラメータ空間での解の探索問題と考えることが可能である。
【０１６３】
タップ長1の瞬時混合の場合、図2（a）のようにパラメータ曲面は単純であり、収束性や安定性に優れる。しかしながら、残響や反射を考慮した畳み込み混合の場合、パラメータ曲面は図2（ｂ）のように複雑な曲面となって、収束の安定性が低下し、真の分離フィルタ以外の極小解（ローカルミニマム）に陥りやすくなる。
【０１６４】
しかしながら、図3のようにビームフォーミング（非特許文献5参照）や先験情報を利用した方法（非特許文献1，6，7）によって得られた結果を適切な初期値W0（z）として設定することで、安定に解が探索できると同時に成分置換も解決することができる（非特許文献8参照）。
【０１６５】
そこで、本発明では、まず、瞬時混合のFDICAにおいて成分置換とスケーリングの不定性の問題を解決し、それによって得られた分離行列を、上述の本発明のICAアルゴリズムの初期値として採択する。そして、この初期値に基づいて、本発明のICAアルゴリズムを実行すれば、成分置換のない真の分離フィルタを求めることができる。
【０１６６】
図3は適切な初期値の設定によりICAアルゴリズムの収束性を模式的に表した図である。図3のように、畳み込み混合のパラメータ曲面は複雑であり、成分置換を伴う解も極小値（ローカルミニマム）として存在する。ICAアルゴリズムは自然勾配法に基づく探索であるため、パラメータ曲面の曲率が小さくなる方向に探索点が移動する。従って、初期値を適切に設定することによってローカルミニマムに陥ることなく成分置換のない真の分離フィルタを探索することができる。
【実施例１】
【０１６７】
図4は、本発明の実施例1に係る耐高残響ブラインド信号分離装置1の全体構成を表す図である。本実施例に係る耐高残響ブラインド信号分離装置1は、2つの音源から発せられる音声信号を2つの受音器により受音して得られる観測信号x1（t），x2（t）を分離して、分離信号y1（t），y2（t）を出力するものである。
【０１６８】
耐高残響ブラインド信号分離装置1は、FFT手段2、サブバンド分離手段3_0〜3_(N-1)、及びIFFT手段4を備えている。
【０１６９】
FFT手段2は、観測信号x1（t），x2（t）のそれぞれについて、フレーム単位でFFTを行うことによってサブバンド分解を行い、N個のサブバンドの観測信号x1（ωn，k），x2（ωn，k）（n＝0，…，N-1）を出力する。ここで、ωnは角周波数、kはフレーム番号、nはサブバンド番号である。
【０１７０】
サブバンド分離手段3_n（n＝0，…，N-1）は、観測信号x1（ωn，k），x2（ωn，k）に対してICAによる分離を行い、分離信号y1（ωn，k），y2（ωn，k）を出力する。
【０１７１】
IFFT手段4は、周波数領域の分離信号y1（ωn，k），y2（ωn，k）（n＝0，…，N-1）に対してそれぞれ逆離散コサイン変換（IFFT）を行い、時間領域の分離信号y1（t），y2（t）を出力する。
【０１７２】
図5は、図4のサブバンド分離手段3_n（n＝0，…，N-1）の構成を表すブロック図である。サブバンド分離手段3_nは、分離行列記憶手段11、分離演算手段12、分離行列初期化手段13、分離行列更新手段14、分離行列合成手段15、成分置換及びスケーリング手段16、及び収束判定手段17を備えている。
【０１７３】
分離行列記憶手段11は、分離行列W（ωn，p）を記憶する。分離演算手段12は、分離行列記憶手段11に記憶された分離行列W（ωn，p）を用いて、観測信号x1（ωn，k），x2（ωn，k）から分離信号y1（ωn，k），y2（ωn，k）を演算する。
【０１７４】
分離行列初期化手段13は、観測信号x1（ωn，k），x2（ωn，k）に基づき、FDICAにより分離行列W（ωn，p）の更新値ΔW（ωn，p）を計算する。また、分離演算手段12及び分離行列合成手段15と協働することによって、各周波数間での行列要素の成分置換のない分離行列の近似値（初期分離行列）を演算する。
【０１７５】
分離行列更新手段14は、観測信号x1（ωn，k），x2（ωn，k）に基づき、式（85）を用いてMDP-ConvFDICAにより分離行列W（ωn，p）の更新値ΔW（ωn，p）を計算する。また、分離演算手段12及び分離行列合成手段15と協働することによって、分離行列の更新を行う。
【０１７６】
分離行列合成手段15は、分離行列初期化手段13又は分離行列更新手段14が出力する分離行列W（ωn，p）の更新値ΔW（ωn，p）と分離行列記憶手段11に記憶された分離行列W（ωn，p）とを加えて新たな分離行列W（ωn，p）を合成し、これにより分離行列記憶手段11の分離行列W（ωn，p）を更新する。
【０１７７】
成分置換及びスケーリング手段16は、FDICAにより得られる分離信号y1（ωn，k），y2（ωn，k）の成分置換を検出し、分離行列記憶手段11に記憶された分離行列W（ωn，p）の成分置換を行う。また、スケーリングの不定性についても解消する。
【０１７８】
収束判定手段17は、分離行列初期化手段13及び分離行列更新手段14による分離行列の更新が収束したか否かを判定し、分離行列初期化手段13及び分離行列更新手段14の動作の切り替え制御を行う。
【０１７９】
以上のように構成された本実施例に係る耐高残響ブラインド信号分離装置1について、以下その動作を説明する。
【０１８０】
図6は、本発明の実施例1に係る耐高残響ブラインド信号分離方法の流れを表すフローチャートである。図6において、ステップS1〜S6は分離行列の初期化処理であり、ステップS7〜S12はMDP-ConvFDICAによる信号分離処理である。
【０１８１】
まず、ステップS1において、時間領域の観測信号x1（t），x2（t）がFFT手段2に入力される。
【０１８２】
ステップS2において、FFT手段2は、観測信号x1（t），x2（t）のそれぞれについて、フレーム単位でFFTを行いサブバンド分解する。これにより、時間領域の観測信号x1（t），x2（t）は、N個のサブバンドの観測信号x1（ωn，k），x2（ωn，k）（n＝0，…，N-1）に変換される。各観測信号x1（ωn，k），x2（ωn，k）は、分離演算手段12及び分離行列初期化手段13に入力される。
【０１８３】
ステップS3において、各サブバンドの分離演算手段12は、分離行列記憶手段11に記憶された分離行列W（ωn，p）を読み出し、式（93）により、分離信号Y（ωn，k）＝［y1（ωn，k），y2（ωn，k）］Tを演算する。尚、初期状態においては、分離行列W（ωn，p）は単位行列に初期化されているものとする。
【０１８４】
ステップS4において、各サブバンドの分離行列初期化手段13は、FDICA（式（25））により、分離行列の更新値ΔW（ωn，p）を算出する。分離行列合成手段15は、分離行列記憶手段11に記憶された分離行列W（ωn，p）を読み出し、分離行列の更新値ΔW（ωn，p）と加えることにより、新たな分離行列Wnew（ωn，p）を生成する。そして、分離行列記憶手段11に記憶された分離行列W（ωn，p）を新たな分離行列Wnew（ωn，p）に更新する。
【０１８５】
ステップS5において、収束判定手段17は、分離行列の更新値ΔW（ωn，p）の大きさ｜ΔW（ωn，p）｜が所定の閾値Wthより小さいか否かを判定する。｜ΔW（ωn，p）｜≧Wthであれば、ステップS1に戻る。｜ΔW（ωn，p）｜＜Wthの場合には、次のステップS6に移行する。
【０１８６】
ステップS6において、成分置換及びスケーリング手段16は、成分置換に関しては、各サブバンドの分離信号y1（ωn，k），y2（ωn，k）の分割ベクトルの実部の分布に対するエントロピーH11（ωn），H22（ωn）を式（33）により計算し、その差ΔH（ωn）を式（34）により計算する。そして、ΔH（ωn）＞0であれば、分離信号y1（ωn，k），y2（ωn，k）の順序を入れ替えるとともに、分離行列記憶手段11に記憶された分離行列W（ωn，p）の行の置換を行う。次に、スケーリングに関して分離行列に対して式（１６）に準じてスケーリングの解決を行う。これにより、FDICAにおける成分置換が解消される。以上で、分離行列初期化手段13による分離行列の初期化処理が終了し、次にMDPに基づくConvFDICAによる信号分離処理（ステップS7）へ移行する。
【０１８７】
ステップS7において、時間領域の観測信号x1（t），x2（t）がFFT手段2に入力される。
【０１８８】
ステップS8において、FFT手段2は、観測信号x1（t），x2（t）のそれぞれについて、フレーム単位でFFTを行いサブバンド分解する。これにより、時間領域の観測信号x1（t），x2（t）は、N個のサブバンドの観測信号x1（ωn，k），x2（ωn，k）（n＝0，…，N-1）に変換される。各観測信号x1（ωn，k），x2（ωn，k）は、分離演算手段12及び分離行列更新手段14に入力される。
【０１８９】
ステップS9において、各サブバンドの分離演算手段12は、分離行列記憶手段11に記憶された分離行列W（ωn，p）を読み出し、式（93）により、分離信号Y（ωn，k）＝［y1（ωn，k），y2（ωn，k）］Tを演算する。
【０１９０】
ステップS10において、各サブバンドの分離行列更新手段14は、MDP-ConvFDICA（式（85））により、分離行列の更新値ΔW（ωn，p）を算出する。分離行列合成手段15は、分離行列記憶手段11に記憶された分離行列W（ωn，p）を読み出し、分離行列の更新値ΔW（ωn，p）と加えることにより、新たな分離行列Wnew（ωn，p）を生成する。そして、分離行列記憶手段11に記憶された分離行列W（ωn，p）を新たな分離行列Wnew（ωn，p）に更新する。
【０１９１】
ステップS11において、IFFT手段4は、全サブバンドの分離信号｛y1（ωn，k）｜n＝0，…，N-1｝，｛y2（ωi，k）｜n＝0，…，N-1｝について、IFFTを行い時間領域の分離信号y1（t），y2（t）に変換する。
【０１９２】
最後に、ステップS12において、IFFT手段4は、時間領域の分離信号y1（t），y2（t）を出力し、再びステップS7に戻り、以降は同じ動作を繰り返す。
【０１９３】
尚、本発明は上述の実施例に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、本発明の構成は、コンピュータによっても実現できる。この場合、各手段が有すべき機能の処理内容は、プログラムによって記述され、このプログラムをコンピュータで実行することにより、上述の処理機能をコンピュータ上で実現することができる。
【０１９４】
〔実験結果〕
最後に、本発明に係る耐高残響ブラインド信号分離方法（MDP-ConvFDICA）の実験結果について示す。ここでは、残響下で1音源2受音器で混合された音声を分離する実験の結果について説明する。
【０１９５】
混合信号は、新聞読み上げ音声コーパス（（社）日本音響学会. 新聞記事読み上げ音声コーパス．JNAS Vols. 1-16, 1997.）から選んだ8秒の音声データに、RWCP実環境音声・音響データベース（Real World Conputing Partnership, "RWCP sound scene database in real acoustic environments," http://tosa.mri.co.jp/sounddb/index.htm）から選んだ残響時間300 msのインパルス応答を畳み込むことにより作成した。ここに、インパルス応答の音源と受音器の配置は図7の通りで、受音器間距離2.83 cm、受音器の中心から音源までの距離は2 mである。音声データとしては、男女各2名による2通りの文を使用し、12通りの混合信号を作成した。尚、標本化周波数は8 kHzとした。
【０１９６】
ここでは、本発明のMDP-ConvFDICAに加えて、比較例として、自然勾配（Natural Gradient）法（S. Amari, "Natural Gradient Works Efficiently in Learning", Neural Computation, Vol.10, pp.251-276, 1998.）、FastICA（A.Hyvarinen, J.Karhunen and E. Oja, "Independent Component Analysis", Jhon Wiley & Sons, pp.165-202, 2001, ISBN 0-471-40540）、及びConvICA（非特許文献4）による分離実験も併せて行った。自然勾配法とFastICAについては、フレーム長を128，256，512，1024，2048［point］と変化させ、フレーム周期をフレーム長の1／4とし、窓関数をHamming窓としてフーリエ変換を行った。学習に用いた信号の長さは3秒、8秒とした。自然勾配法については、初期荷重は［−1，1］の乱数、更新回数は最大100回とした。FastICAについては、初期荷重は［−1，1］の乱数、更新回数は最大1000回、収束判定基準は更新前後の分離荷重の内積の絶対値が0.999999以上となることとした．
【０１９７】
ConvFDICA及び本発明のMDP-ConvFDICAについては，フレーム長を128、フレーム周期をフレーム長の1／4とし、窓関数をHamming窓としてフーリエ変換を行った。また、それぞれのフレーム長に対して、分離行列のタップ長を4，8，16と変化させた。学習に用いた信号の長さは3秒とした。ConvFDICAのタップ長4，8，16に対する更新回数はそれぞれ1000回，2000回，4000回で、探索ステップ幅は0.007, 0.002, 0.002とした。
【０１９８】
本発明のMDP-ConvFDICAのタップ長4，8，16に対する更新回数はそれぞれ最大1000回，2000回，4000回で，探索ステップ幅ηは0.02, 0.006, 0.002で、内分比γは0.1とした。
【０１９９】
尚、ConvFDICAとMDP-ConvFDICAについてはブロック化してFFTで計算することにより高速化した。ブロックサイズはタップ長と同じ長さとした。最終的な分離行列は、最も良い分離性能を与える更新回数での分離行列とした。
【０２００】
〔１〕分離性能の比較
本実験では分離性能の評価尺度として，Noise Reduction Rate（NRR；≡出力SNR［dB］−入力SNR［dB］）を用いる（非特許文献5参照）。
NRRは、目的音声と干渉音声（雑音）の比をSNRとして、次式（91）のように定義される。
【０２０１】
【数５７】
2006308955.tif 000059
【０２０２】
ここで、OSNRiとISNRiはそれぞれ出力SNRと入力SNRでi≠jである。また、gijは分離行列W(z)と混合行列A(z)の積G(z)= W(z) A(z)のi行j列要素である。図8と図9に各方法による分離性能を示す。図8では横軸はタップ長を示し、縦軸はNNRを示す。また、図9では横軸は窓幅を示し、縦軸はNNRを示す。これらの図から、FastICA、自然勾配法（NG）、MDP-ConvFDICA、ConvFDICAの順に、NNRの値は向上することが読みとれる。FastICAと自然勾配法（NG）の場合、狭帯域分離フィルタのタップ長をM=１としたことに相当する瞬時混合モデルとなっているため、白色化歪みは起きない。また、MDP-ConvFDICAの場合、これまで述べてきたように白色化歪みが起きないように工夫されている。しかし、ConvFDICAの場合、上述で再三指摘してきたように、白色化歪みが起きる。NNRは、従来、分離能力を計るための評価量として用いられてきたが、白色化歪みが大きくなると、NNRも大きくなると云う欠点がある。そのため、白色化歪みの影響も含めて分離能力を評価するSDR（Signal Distortion Ratio）が以下のように提案されている。
【０２０３】
〔２〕分離信号の音質に関する比較
本実験では、分離信号の白色化歪みに関する評価尺度として、SDR（Signal Distortion Ratio）を用いる（高谷智哉，西川剛樹，猿渡洋，鹿野清宏，「SIMOモデルを用いた高忠実度なブラインド音源分離」，信学技報，EA2002-108，pp.19-24，2003）。SDRを下式（92）のように定義する。
【０２０４】
【数５８】
2006308955.tif 000060
【０２０５】
これは、分離信号の音質を評価するものであり、この値が高い方が歪みなく信号を分離できたことになる。図10，図11に分離信号のSDRを示す。図10では横軸はタップ長を示し、縦軸はSDRを示している。図11では横軸は窓幅を示し、縦軸はSDRを示している。これらの図より、FastICAと自然勾配法（Natural Gradient：NG）では高い値を保っており、本発明のMDP-ConvFDICAでもそれほどの劣化は見られないが、ConvFDICAでは低い値を示しており、音質が劣化していることが読み取れる。このことから、MDP-ConvFDICAはConvFDICAよりSDRが約６ｄB高く原音を忠実に再現できていることが分かる。
【０２０６】
以上より、本発明に係る耐高残響ブラインド信号分離装置及び耐高残響ブラインド信号分離方法によれば、以下のような効果を得ることができる。
【０２０７】
（1）混合過程を周波数領域で畳み込みモデルで厳密に定式化した分離モデルを用いることによって、高残響下でも信号を良好に分離することができる。
（2）TDICAやConvFDICAで問題となっていた白色化による歪みが解消される。
（3）本発明における初期値設定法を用いることで、成分置換のない最適解が得られる。これにより、これまで課題となっていた成分置換の問題が解決される。
（4）ConvICAでは未解決であった成分置換やスケーリングの問題が解消される。
【０２０８】
また、本発明に係る耐高残響ブラインド信号分離方法が高残響下で有効に機能することについて、実験により以下の点が確認された。
【０２０９】
（1）インパルス応答を各周波数ビンに分割するので、残響が長くても各周波数ビンでフィルタ長は短くてすむ。これにより、TDICAやMDPに基づくTDICAで残響時間が長い、すなわちフィルタ長を長く取ると収束しないという欠点が軽減される。
（2）インパルス応答を各周波数ビンに分割するので、短い窓幅でも全ての残響を包含できる。FDICAの窓幅を大きくすると十分な統計量が得られないという欠点が軽減される。
（3）各周波数ビンにMDPを適用することにより、白色化による歪みの影響を回避できる。これにより、ConvFDICAの白色化による歪みの影響が軽減される。
（4）本発明のMDP-ConvFDICAによる分離信号には、白色化歪みによる劣化がなく、原音をほぼ忠実に再現でき、音質もConvFDICAに比べて6dB向上することが確認された。
【図面の簡単な説明】
【０２１０】
【図１】時間窓掛け離散フーリエ変換を表す図である。
【図２】瞬時混合の場合と畳み込み混合の場合のパラメータ曲面の状態を表す模式図である。
【図３】適切な初期値の設定による解の収束性について説明する図である。
【図４】本発明の実施例1に係る耐高残響ブラインド信号分離装置1の全体構成を表す図である。
【図５】図4のサブバンド分離手段3_n（n＝0，…，N−1）の構成を表すブロック図である。
【図６】本発明の実施例1に係る耐高残響ブラインド信号分離方法の流れを表すフローチャートである。
【図７】実験における受音器と音源との配置を示す図である。
【図８】FDICAのNNRの測定結果を示す図である。
【図９】ConvFDICAとMDP-ConvFDICAのNNRの測定結果を示す図である。
【図１０】FDICAのSDRの測定結果を示す図である。
【図１１】ConvFDICAとMDP-ConvFDICAのSDRの測定結果を示す図である。
【図１２】R1＝2、R2＝2の場合のBSDの構成を表す図である。
【図１３】TDICAからみた混合過程と分離過程を表す図である。
【図１４】成分置換の問題を問題について説明する図である。
【図１５】フーリエ変換の分析窓幅とインパルス応答との関係を表す図である。
【符号の説明】
【０２１１】
1  耐高残響ブラインド信号分離装置
2  FFT手段
3_0，3_2，…，3_(N-1)  サブバンド分離手段
4  IFFT手段
11  分離行列記憶手段
12  分離演算手段
13  分離行列初期化手段
14  分離行列更新手段
15  分離行列合成手段
16  成分置換及びスケーリング手段
17  収束判定手段
【図１】
2006308955.tif 000061
【図２】
2006308955.tif 000062
【図３】
2006308955.tif 000063
【図４】
2006308955.tif 000064
【図５】
2006308955.tif 000065
【図６】
2006308955.tif 000066
【図７】
2006308955.tif 000067
【図８】
2006308955.tif 000068
【図９】
2006308955.tif 000069
【図１０】
2006308955.tif 000070
【図１１】
2006308955.tif 000071
【図１２】
2006308955.tif 000072
【図１３】
2006308955.tif 000073
【図１４】
2006308955.tif 000074
【図１５】
2006308955.tif 000075