(57)【要約】【課題】混合信号について高精度な分離処理および音源方向推定処理を可能とする構成を提供する。【解決手段】複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ）を適用した学習処理により分離行列を求めて分離信号を生成するとともに、時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出する。具体的には、相互共分散行列の要素同士の位相差情報を利用して、位相角と周波数との比例定数（直線の傾き）を算出し、比例定数を利用して各分離信号に対応する音源方向を算出する。【選択図】図１２
2010020294.tif 000002
【特許請求の範囲】
【請求項１】
複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求める学習処理部と、
前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成する分離処理部と、
前記分離処理部において生成した分離信号の音源方向を算出する音源方向推定部を有し、
前記音源方向推定部は、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出する構成である信号処理装置。
【請求項２】
前記分離処理部は、
前記複数音源からの出力を混合した混合信号に対して短時間フーリエ変換（ＳＴＦＴ）を実行して時間周波数領域の観測信号スペクトログラムを生成し、
さらに、前記観測信号スペクトログラムに対して、前記分離行列を適用して時間周波数領域の分離信号スペクトログラムを生成する構成であり、
前記音源方向推定部は、
時間周波数領域において対応する時間区間の観測信号スペクトログラムと分離信号スペクトログラムとの間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出する構成である請求項１に記載の信号処理装置。
【請求項３】
前記音源方向推定部は、
時間周波数領域の観測信号の異なる時間の観測信号間の共分散行列を算出し、算出した共分散行列と、前記学習処理部の算出した分離行列に対するエルミート転置行列を適用して、観測信号と分離信号との間の相互共分散行列を算出する構成である請求項１に記載の信号処理装置。
【請求項４】
前記音源方向推定部は、
前記位相差の算出処理において、各音源（チャンネル）について、異なる位置に設定された信号入力源であるマイクペアに対応する位相差を個別に算出する構成である請求項１に記載の信号処理装置。
【請求項５】
前記音源方向推定部は、
相互共分散行列の要素同士の位相差情報を利用して、位相角と周波数との比例定数（直線の傾き）を算出する直線当てはめ処理を実行して、前記比例定数を利用して各分離信号に対応する音源方向を算出する請求項１に記載の信号処理装置。
【請求項６】
前記音源方向推定部は、
前記直線当てはめ処理として、
低い周波数から順に、位相角と周波数との対応情報のアンラップ処理を実行して前記比例定数の更新処理を行う再帰法を実行して比例定数（直線の傾き）を算出する請求項５に記載の信号処理装置。
【請求項７】
前記音源方向推定部は、
前記直線当てはめ処理として、
前記比例定数（直線の傾き）の値を予め設定した最小値から最大値まで少しずつ変化させて、前記位相差情報との誤差を計算し、誤差の少ない比例定数を選択するリニアサーチ法を実行して比例定数（直線の傾き）を算出する請求項５に記載の信号処理装置。
【請求項８】
前記音源方向推定部は、
前記直線当てはめ処理として、
前記比例定数（直線の傾き）の値を予め設定した最小値から最大値まで少しずつ変化させて、前記位相差情報との誤差を計算し、誤差の少ない比例定数を選択するリニアサーチ法を実行し、
さらに、前記リニアサーチ法によって得られた比例定数を初期値として、低い周波数から順に、位相角と周波数との対応情報のアンラップ処理を実行して前記比例定数の更新処理を行う再帰法を実行して比例定数（直線の傾き）を算出する請求項５に記載の信号処理装置。
【請求項９】
前記音源方向推定部は、
各音源（チャンネル）について、異なる位置に設定された信号入力源であるマイクペア各々に対応する複数の比例定数（直線の傾き）を算出し、算出した複数の比例定数（直線の傾き）から当てはめ誤差の少ないデータを選択し、選択データの平均値を算出して最終的な比例定数（直線の傾き）を算出する請求項５に記載の信号処理装置。
【請求項１０】
前記信号処理装置は、さらに、
前記直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して、前記直線からのずれの大きい位相角データを持つ信号分離結果の周波数成分を減衰させる周波数マスキングを実行する構成である請求項５に記載の信号処理装置。
【請求項１１】
前記信号処理装置は、さらに、
前記直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して、前記分離行列を生成し、生成した分離行列を前記学習処理部において利用する構成である請求項５に記載の信号処理装置。
【請求項１２】
前記音源方向推定部は、
各マイクペアを結ぶベクトル［ｐｉｍ］について大きさを１に正規化して横方向に並べて生成した行列［Ｐｋ］と、各マイクペアについての経路差［Ｄｉｍｋ］をマイク間隔で正規化して横方向に並べて生成したベクトル［Ｄｋ］を算出し、算出した［Ｐｋ］とベクトル［Ｄｋ］を適用して、音源方向の角度、または、マイクペアから見た音源方向を示すベクトルを算出する構成である請求項１に記載の信号処理装置。
【請求項１３】
信号処理装置において実行する信号処理方法であり、
学習処理部が、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求める学習処理ステップと、
分離処理部が、前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成する分離処理ステップと、
音源方向推定部が、前記分離処理ステップにおいて生成した分離信号の音源方向を算出する音源方向推定ステップを有し、
前記音源方向推定ステップは、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出するステップである信号処理方法。
【請求項１４】
信号処理装置において信号処理を実行させるプログラムであり、
学習処理部に、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求めさせる学習処理ステップと、
分離処理部に、前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成させる分離処理ステップと、
音源方向推定部に、前記分離処理ステップにおいて生成した分離信号の音源方向を算出させる音源方向推定ステップを有し、
前記音源方向推定ステップは、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出させるステップであるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理装置、および信号処理方法、並びにプログラムに関する。さらに、詳細には、本発明は、複数の信号が混合された信号を独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ；ＩＣＡ）を用いて分離し、さらに音源方向を推定する信号処理装置、および信号処理方法、並びにプログラムに関する。
【背景技術】
【０００２】
本発明は、音源方向（マイクから見て音がどの方向からやってくるかＤｉｒｅｃｔｉｏｎ  Ｏｆ  Ａｒｒｉｖａｌ：ＤＯＡともいう）を推定する技術に関するものである。特に独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ：ＩＣＡ）を利用して音源方向をリアルタイムに推定する技術に関する。
【０００３】
本発明は音源方向推定およびＩＣＡの両方に関わるため、従来技術として双方について以下の順番で説明していく。その後、リアルタイムのＩＣＡに音源方向推定を組み込んだ場合の問題点について説明する。説明は、
Ａ．従来技術の説明
Ｂ．従来技術の問題点
上記Ａ，Ｂの順番で説明する。
【０００４】
［Ａ．従来技術の説明］
まず、従来技術として以下のＡ１〜Ａ４について説明する。
Ａ１．マイク間の位相差を利用した音源方向推定
Ａ２．ＩＣＡの説明
Ａ３．ＩＣＡを用いた音源方向推定
Ａ４．ＩＣＡのリアルタイム化
以降では、上の４点についてそれぞれ説明する。
【０００５】
［Ａ１．マイク間の位相差を利用した音源方向推定］
音源方向を推定するための最も基本的な方式は、複数のマイクを使用して到達時間差や位相差などを利用するものである。この方式は、以下の文献でも従来技術として紹介されている。
特許文献１（特開２００５−７７２０５号公報）
特許文献２（特開２００５−４９１５３号公報）
特許文献３（特許第３０６９６６３号公報）
【０００６】
図１を用いてこの基本的な音源方向の推定方式について説明する。図１に示すように、一つの音源１１があり、また、２つのマイク１２，１３が間隔ｄで設置されている環境を考える。マイク１２，１３から音源１１までの距離がマイク間隔ｄに対して十分大きい場合、音波は平面波と近似することができ、この場合、２つのマイクには同一の角度θで音が到達すると仮定できる。図１において、角度θは、両マイクを結ぶ直線（マイクペア方向）と直交する方向を０としている。
【０００７】
音源１１からマイク１，１２までの経路と、音源１１からマイク２，１３までの経路とを比較すると、後者の方が［ｄ・ｓｉｎθ］だけ長いと近似することができるため、マイク２，１３で観測される信号はマイク１，１２の観測信号よりも［ｄ・ｓｉｎθ］に比例した分だけ位相が遅れて観測される。
【０００８】
なお、図１において、角度θは、両マイクを結ぶ直線（マイクペア方向）と直交する方向を０としているが、マイクペアの方向と音波進行方向とのなす角度［θ'］を用いた場合は、音源と２つのマイク間距離の差は［ｄ・ｃｏｓθ'］で表わすことができる。
【０００９】
マイク１とマイク２の観測信号をそれぞれｘ１（ｔ），ｘ２（ｔ）として、音源１１からマイク１，２までの経路差を［ｄ・ｓｉｎθ］とした場合の関係式は、以下に示す式［１．１］となる。
【００１０】
【数１】
2010020294.tif 000003
【００１１】
上記式において、ｔは、離散時間であり、サンプルの番号を表わす。マイク２の観測信号：ｘ２（ｔ）は音源からの距離がｘ１（ｔ）よりも長い分だけ位相が遅れるため、減衰を無視すれば、式［１．１］のように表わすことができる。
【００１２】
従って、各マイクの観測信号によって位相の遅れを検出できれば、上記式に従って、音源方向を決定する角度［θ］を算出することができる。以下では、観測信号を時間周波数領域に変換して音源方向を算出する方法について説明する。
【００１３】
複数の観測信号ｘ１，ｘ２に対して後述の短時間フーリエ変換（Ｓｈｏｒｔ−ｔｉｍｅ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ：ＳＴＦＴ）をかけた結果をそれぞれＸ１（ω，ｔ），Ｘ２（ω，ｔ）とする。ω，ｔは、それぞれ周波数ビンの番号とフレームの番号を表わす。また、変換前の観測信号［ｘ１，ｘ２］を時間領域の信号、短時間フーリエ変換後の信号［Ｘ１，Ｘ２］を時間周波数領域（またはＳＴＦＴ領域）の信号という。
【００１４】
時間領域での位相遅れは時間周波数領域での複素数倍に相当するため、時間領域における関係式［１．１］は時間周波数領域では、以下に示す式［１．２］のように表わすことができる。
【００１５】
【数２】
2010020294.tif 000004
【００１６】
音源方向を示す角度［θ］を含む項を取り出すためには、以下に示す式［１．３］のような操作を行なえばよい。
【００１７】
【数３】
2010020294.tif 000005
【００１８】
ただし、上記式において、ａｎｇｌｅ（）は複素数の偏角を−π〜＋πの範囲で求める関数であり、上線のついたＸ２は、Ｘ２の共役複素数を表わす。最終的に、以下に示す式［１．４］によって音源方向を推定することができる。
【００１９】
【数４】
2010020294.tif 000006
【００２０】
ただし、この式において、
ａｓｉｎはｓｉｎの逆関数を表わす。
また、ｈａｔ（θ（ω））は、周波数ビンωにおける角度θが観測値から推定された値であることを意味する。なお、本文中におけるｈａｔは記号（＾）を意味する。
【００２１】
上記式［１．４］は、ある特定のω（周波数ビン番号）とｔ（フレーム番号）における式だが、複数のωやｔについて音源方向を算出してから平均を取ることにより、安定したθの値を得ることができる。また、マイクをｎ個（ｎ＞２）用意し、ｎ（ｎ＋１）／２通りのペアそれぞれについて音源方向を計算することもできる。
【００２２】
以下に示す式［１．５］〜［１．７］は、複数マイク・複数フレームを用いた場合の式である。
【００２３】
【数５】
2010020294.tif 000007
【００２４】
各マイクの観測信号を要素とするベクトルを式［１．５］のように定義し、そのベクトルの平均を式［１．６］で定義する。ただし、この式のＥｔ［］は、ある区間のフレームでの平均であり、式［１．７］で定義される。式［１．４］において、Ｘ１（ω，ｔ），Ｘ２（ω，ｔ）の代わりに式［１．６］のＥｔ［Ｘｉ（ω，ｔ）］，Ｅｔ［Ｘｍ（ω，ｔ）］を用いることで、ｉ番目のマイクとｍ番目のマイクとのペアから計算される、ｔ０番目からｔ１番目までのフレームに対応した平均的な音源方向が求まる。
【００２５】
この、位相差を用いた音源方向推定方法は、後述のＩＣＡを用いた方式に比べて処理量が少ないため、リアルタイム（１フレーム未満の遅延）かつ高頻度（フレーム毎）に音源方向を算出することができる。
【００２６】
しかしその一方で、複数の音源が同時に鳴っている環境では音源方向を求めることができなくなる。また、音源が１個であっても、反射や残響が大きい環境では、方向の精度が低下するという問題がある。
【００２７】
［Ａ２．ＩＣＡの説明］
ここでは、独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ：ＩＣＡ）を用いた音源方向推定について説明する前に、まずＩＣＡそのものについて説明する。
【００２８】
ＩＣＡとは、多変量分析の一種であり、信号の統計的な性質を利用して多次元信号を分離する手法のことである。ＩＣＡ自体の詳細については、例えば非特許文献１［『入門・独立成分分析』（村田昇著、東京電機大学出版局）］などを参照されたい。
【００２９】
以下では、音信号のＩＣＡ、特に時間周波数領域のＩＣＡについて説明する。
図２に示すように、Ｎ個の音源から異なる音が鳴っていて、それらをｎ個のマイクで観測するという状況を考える。音源が発した音（原信号）がマイクに届くまでには、時間遅れや反射などがある。従って、マイクｊで観測される信号（観測信号）は以下に示す式［２．１］のように、原信号と伝達関数（transfer function）との畳み込み演算（convolution）を全音源について総和した式として表わすことができる。この混合を以下では「畳み込み混合」（convolutive mixtures）と呼ぶ。
【００３０】
【数６】
2010020294.tif 000008
【００３１】
また、全てのマイクについての観測信号を一つの式で表わすと、以下に示す式［２．２］のように表わせる。
【００３２】
【数７】
2010020294.tif 000009
【００３３】
ただし、ｘ（ｔ），ｓ（ｔ）はそれぞれｘｋ（ｔ），ｓｋ（ｔ）を要素とする列ベクトルであり、Ａ［ｌ］はａ［ｌ］ｋｊを要素とするｎ×Ｎの行列である。以降では、ｎ＝Ｎとする。
【００３４】
時間領域の畳み込み混合は、時間周波数領域では瞬時混合で表わされることが知られており、その特徴を利用したのが時間周波数領域のＩＣＡである。
【００３５】
時間周波数領域ＩＣＡ自体については、非特許文献２［『詳解  独立成分分析』の「１９．２．４．フーリエ変換法」］や、特許文献４（特開２００６−２３８４０９『音声信号分離装置・雑音除去装置および方法』）などを参照されたい。
【００３６】
以下では、主に本発明と関係ある点を説明する。
上記の式［２．２］の両辺を短時間フーリエ変換すると、以下に示す式［３．１］が得られる。
【００３７】
【数８】
2010020294.tif 000010
【００３８】
上記式［３．１］において、
ωは周波数ビンの番号、
ｔはフレームの番号、
である。
【００３９】
ωを固定すると、この式は瞬時混合（時間遅れのない混合）と見なせる。そこで、観測信号を分離するには、分離結果［Ｙ］の算出式［３．５］を用意した上で、分離結果：Ｙ（ω，ｔ）の各成分が最も独立になるように分離行列Ｗ（ω）を決める。
【００４０】
従来の時間周波数領域ＩＣＡでは、パーミュテーション問題と呼ばれる、「どの成分がどのチャンネルに分離されるか」が周波数ビンごとに異なるという問題が発生していたが、本願と同一発明者による前の特許出願である特許文献４［特開２００６−２３８４０９『音声信号分離装置・雑音除去装置および方法』］に示した構成によって、このパーミュテーション問題は、ほぼ解決することができた。本発明でもこの方法を用いるため、特許文献４［特開２００６−２３８４０９］に開示したパーミュテーション問題の解決手法ついて簡単に説明する。
【００４１】
特許文献４［特開２００６−２３８４０９］では、分離行列Ｗ（ω）を求めるために、以下に示す式［４．１］から式［４．３］までを分離行列Ｗ（ω）が収束するまで（または一定回数）繰り返し実行する。
【００４２】
【数９】
2010020294.tif 000011
【００４３】
この繰り返し実行を以降では「学習」と呼ぶ。ただし、式［４．１］〜式［４．３］は、全ての周波数ビンに対して行ない、さらに式［４．１］は、蓄積された観測信号の全てのフレームに対しても行なう。また、式［４．２］において、Ｅｔ［・］は全フレームについての平均を表わす。Ｙ（ω，ｔ）の右上についている上付きのＨはエルミート転置（ベクトルや行列の転置を取ると共に、要素を共役複素数に変換する）である。
【００４４】
分離結果Ｙ（ｔ）は式［４．４］で表わされる。また、分離結果の全チャンネル・全周波数ビンの要素を並べたベクトルである。φω（Ｙ（ｔ））は、式［４．５］で表わされるベクトルである。各要素φω（Ｙｋ（ｔ））はスコア関数と呼ばれ、Ｙｋ（ｔ）の多次元（多変量）確率密度関数（ＰＤＦ）の対数微分である（式［４．６］）。多次元ＰＤＦとして、例えば式［４．７］で表わされる関数を用いることができ、その場合、スコア関数φω（Ｙｋ（ｔ））は式［４．９］のように表わせる。ただし、‖Ｙｋ（ｔ）‖２はベクトルＹｋ（ｔ）のＬ−２ノルム（全要素の２乗和を求め、さらに平方根をとったもの）である。Ｌ−２ノルムを一般化したＬ−ｍノルムは式［４．８］で定義される。式［４．７］および式［４．９］のγは、Ｙｋ（ω，ｔ）のスケールを調整するための項であり、例えばｓｑｒｔ（Ｍ）（周波数ビン数の平方根）といった適切な正の定数を代入しておく。式［４．３］のηは学習率や学習係数と呼ばれる正の小さな値（例えば０．１程度）である。これは、式［４．２］で計算されたΔＷ（ω）を分離行列Ｗ（ω）に少しずつ反映させるために用いられる。
【００４５】
なお、式［４．１］は一つの周波数ビンにおける分離（図３（Ａ）参照）を表わしているが、全周波数ビンの分離を一つの式で表わす（図３（Ｂ）参照）ことも可能である。
【００４６】
そのためには、上述した式［４．４］で表わされる全周波数ビンの分離結果Ｙ（ｔ）および、以下に示す式［４．１１］で表わされる観測信号Ｘ（ｔ）、さらに式［４．１０］で表わされる全周波数ビン分の分離行列を用いればよく、それらのベクトルと行列を用いることで、分離は式［４．１２］のように表わすことができる。本発明は、必要に応じて式［４．１］と式［４．１１］とを使い分ける。
【００４７】
【数１０】
2010020294.tif 000012
【００４８】
なお、図３に示したＸ１〜ＸｎおよびＹ１〜Ｙｎの図はスペクトログラムと呼ばれ、短時間フーリエ変換（ＳＴＦＴ）の結果を周波数ビン方向とフレーム方向とに並べたものである。縦方向が周波数ビン、横方向がフレームである。式［４．４］や式［４．１１］では低い周波数を上に書いてあるが、スペクトログラムでは低い周波数を下に描いてある。
【００４９】
ＩＣＡはブラインド音源分離（Blind Source Separation: BSS）の一種であり、分離処理においては、信号源についての知識（音源の周波数分布や到来方向など）を必要とせず、音速・マイク間隔・マイクゲイン・伝達関数などが未知でも分離が行えるという特徴がある。さらに、信号源に対して、独立性以外の仮定（例えばスパース性など）も不要であり、複数の音源が鳴りっぱなしの状態でかつ周波数領域で混合が発生している場合（すなわちスパース性の仮定が成立していない場合）でも分離が可能である。これらの特徴は、後述の「方向を推定してから音源分離」という方式に対する利点となる。
【００５０】
［Ａ３．ＩＣＡを用いた音源方向推定］
ＩＣＡを用いて音源方向を推定する方式については、例えば特許文献５（特許第３８８１３６７号公報）に記載がある。
【００５１】
以下、上記特許文献５において開示された分離行列Ｗ（ω）から音源方向を求める方法について説明する。学習が収束した後の分離行列Ｗ（ω）は、混合行列Ａ（ω）の逆の過程を表わしている。従って、Ｗ（ω）の逆行列であるＷ（ω）−１を求めると、それは周波数ビンωに対応した伝達関数を表わしていると考えることができる。
【００５２】
以下に示す式［５．１］のように、分離行列Ｗ（ω）の逆行列をＢ（ω）とおき、さらにＢ（ω）の各要素を式［５．２］のように定める。
【００５３】
【数１１】
2010020294.tif 000013
【００５４】
ここで、分離行列Ｗ（ω）のｋ行目（横の列）および逆行列Ｂ（ω）のｋ列目（縦の列）に着目すると（図４（ａ）参照）、分離行列Ｗ（ω）のｋ行目は観測信号のベクトルＸ（ω，ｔ）からｋチャンネル目の分離結果（＝推定音源）Ｙｋ（ω，ｔ）を生成するベクトルであるのに対し、逆行列Ｂ（ω）である推定伝達関数行列Ｂ（ω）のｋ列目は推定音源Ｙｋ（ω，ｔ）から各マイクまでの伝達関数（図４（ｂ）参照）を表わしている。
【００５５】
すなわち、ｍ番目のマイクで観測される観測信号Ｘｍ（ω，ｔ）のうち、音源Ｙｋ（ω，ｔ）に由来する成分を成分信号Ｘｍｋ（ω，ｔ）と書くと、成分信号Ｘｍｋ（ω，ｔ）は以下に示す式［５．３］のように表わすことができる。
【００５６】
【数１２】
2010020294.tif 000014
【００５７】
分離行列Ｗ（ω）の逆行列Ｂ（ω）を用いると、観測信号のうち、各音源に由来する成分をそれぞれ求めることができる。従って、最初に説明した「１．マイク間の位相差を利用した音源方向推定」と同様の方法で、音源の方向を算出することができる。
【００５８】
例えば、ｉ番目のマイクとｍ番目のマイクとを用いて分離結果Ｙｋ（ω，ｔ）の音源方向であるｈａｔ（θｉｍｋ（ω））を求めるためには、先に説明した式［１．４］の２つのマイクの観測信号に対応する短時間フーリエ変換後の観測信号Ｘ１（ω，ｔ）とＸ２（ω，ｔ）とに対して、上記式［５．３］から算出される成分信号Ｘｉｋ（ω，ｔ）とＸｍｋ（ω，ｔ）とをそれぞれ代入すればよく、結果として、以下に示す式［５．４］を得る。
【００５９】
【数１３】
2010020294.tif 000015
【００６０】
ただしこの式において、ｄｉｍはｉ番目のマイクとｍ番目のマイクとの距離である。
上記式［５．４］においては、フレーム番号［ｔ］が消えていることに注意されたい。すなわち、同一の分離行列Ｗ（ω）を使って分離結果Ｙ（ω，ｔ）＝Ｗ（ω）Ｘ（ω，ｔ）を計算している限り、どのフレームｔに対しても同一の音源方向が算出される。
【００６１】
このような、ＩＣＡの分離行列の逆行列から音源方向を算出する方法は、「１．マイク間の位相差を利用した音源方向推定」と異なり、複数の音源が同時に鳴っていても、それぞれの音源方向を求めることができるという利点がある。
【００６２】
また、分離結果と音源方向とがどちらも分離行列Ｗ（ω）から計算されるため、分離結果の各チャンネルがそれぞれどの方向かという対応付けが始めからなされているのも利点である。つまり、多音源対応の方向推定手段として、ＩＣＡ以外の音源分離手段と音源方向推定手段とを組み合わせた場合は、方式の組み合わせによってはチャンネル同士の対応付けが後段で別途必要となり、そこで間違いが発生する可能性もあるが、ＩＣＡ  を用いた方式ではそのような後段での対応付けは不要である。
【００６３】
なお、式［５．４］に現れる音速Ｃやマイク間隔ｄｉｍなどは、式の上での値と実環境の値とが食い違う場合がある。そのような場合、ＩＣＡを用いた方式では、音源方向自体には推定誤差が発生するものの、方向の誤差は音源分離の精度には影響しない。また、マイクのゲインには個体差があり、音源分離や方向推定の方式によってはゲインのばらつきが分離や推定の精度に悪影響を与える場合もあるが、ＩＣＡを用いた方式ではマイクゲインのばらつきは分離性能にも方向推定にも悪影響を与えない。
【００６４】
［Ａ４．ＩＣＡのリアルタイム化］
「Ａ２．ＩＣＡの説明」で説明した式［４．１］から式［４．３］までを分離行列Ｗ（ω）が収束するまで（または一定回数）繰り返し実行する学習処理は、バッチ処理によって行われる。すなわち、観測信号全体を蓄積した後で、式［４．１］〜式［４．３］を繰り返す方式である（式［４．１］〜式［４．３］の繰り返し実行は、前述したように、学習と呼ばれる。
【００６５】
このバッチ処理は、ある工夫を行なうことで、リアルタイム（低遅延）の音源分離に適用することができる。リアルタイム処理方式を実現した音源分離処理例として、本願と同一発明者の先の特許出願である「特願２００６−３３１８２３  リアルタイム音源分離装置および方法」に開示した構成について説明する。
【００６６】
特願２００６−３３１８２３で開示した処理方式は、図５に示すように、観測信号のスペクトログラムを重複のある複数のブロック１〜Ｎに分割し、ブロックごとに学習を行なって分離行列を求めている。ブロックに重複を持たせる理由は、分離行列の精度と更新頻度とを両立させるためである。
【００６７】
なお、特願２００６−３３１８２３より以前に開示されていたリアルタイムＩＣＡ（ブロック単位のＩＣＡ）ではブロック間に重複がなかった。従って、分離行列の更新間隔を短くするにはブロック長（＝観測信号を蓄積する時間）を短くする必要があったが、ブロック長を短くすると分離精度は低下するという問題があった。
【００６８】
各ブロックから求まった分離行列は、それより以降の観測信号へ適用する（同一ブロックへは適用しない）ことで分離結果を生成する。ここでは、その方式を「ずらし適用」と呼ぶ。
【００６９】
図６は、「ずらし適用」の説明である。現時点において、ｔ番目のフレームの観測信号Ｘ（ｔ）６２が入力されたとする。この時点では、観測信号Ｘ（ｔ）を含むブロック（例えば現在時刻を含む観測信号ブロック６６）に対応した分離行列はまだ求まっていない。そこで、ブロック６６の代わりに、それより前のブロックである学習データブロック６１から学習された分離行列を観測信号Ｘ（ｔ）に乗じることで、Ｘ（ｔ）に対応した分離結果、すなわち現在時刻の分離結果Ｙ（ｔ）６４を生成する。なお、学習データブロック６１から学習された分離行列は、フレームｔの時点で既に求まっているものとする。
【００７０】
前述の通り、分離行列は混合過程の逆の過程を表わしていると考えられる。
そのため、学習データのブロッマ設定区間６１の観測信号と現在時刻の観測信号６２とで混合過程が同一であれば（例えば、音源とマイクとの位置関係が変化していなければ）、異なる区間で学習された分離行列を適用しても信号を分離することができ、そうすることで遅延の少ない分離を実現することができる。
【００７１】
重複のあるブロックから分離行列を求める処理は、例えば一般的なＰＣにおいて所定のプログラムに従った演算によって実行させることが可能である。特願２００６−３３１８２３において開示した構成では、重複のあるブロックから分離行列を求めるためスレッドと呼ばれる処理の単位を、時刻をずらして複数並列に起動する方式を提案している。図７を参照してこの並列処理方式について説明する。
【００７２】
図７には、処理単位としてのスレッド各々の時間経過に従った処理推移を示している。図７には６つのスレッド１〜６を示している。各スレッドはＡ）蓄積、Ｂ）学習、Ｃ）待機、という３つの状態を繰り返す。すなわち、スレッド長は、Ａ）蓄積、Ｂ）学習、Ｃ）待機の３つの処理のトータル時間長に対応する。図７の左から右に時間が推移する。
【００７３】
「Ａ）蓄積」は図７の濃い灰色の区間であり、スレッドがこの状態にあるときは、観測信号を蓄積する。蓄積を開始する時刻をスレッドごとにずらすことで、図５の重複つきブロックが表現できる。図７では蓄積時間の１／４だけずらしているので、１つのスレッドの蓄積時間を例えば４秒とすると、各スレッド間のずらし時間は１秒となる。
【００７４】
観測信号を一定時間（例えば４秒）蓄積すると、各スレッドは「Ｂ）学習」に状態を遷移させる。「Ｂ）学習」は図７の薄い灰色の区間であり、この状態にあるときは、蓄積された観測信号に対して先に説明した式［４．１］〜式［４．３］を繰り返し実行する。
【００７５】
学習（式［４．１］〜式［４．３］の繰り返し）によって分離行列Ｗが十分に収束したら（または単に一定回数繰り返したら）、学習を終了し、スレッドは「Ｃ）待機」状態に遷移する（図７の白い区間）。「待機」は、蓄積開始時刻および学習開始時刻をスレッド間で一定の間隔に保つためであり、その結果、学習終了時刻（＝分離行列が更新される時刻）もほぼ一定に保たれる。
【００７６】
学習で求まった分離行列Ｗは、次のスレッドの学習が終了するまでの間、分離に使用される。すなわち、図６の分離行列６３として使用される。図７の最下段に示す分離行列の時間推移に示す適用分離行列規定区間７１〜７３に使用される分離行列について説明する。
【００７７】
システムが起動してから、最初の分離行列が学習されるまでの適用分離行列規定区間７１では、図６の分離行列６３として初期値（例えば単位行列）が使用される。図７に示すスレッド１の学習が終了してからスレッド２の学習が終了するまでの区間７２では、スレッド１の観測信号蓄積区間７４に由来する分離行列が図６の分離行列６３として使用される。図７の区間７２に示す数字「１」は、この期間に使用される分離行列Ｗがスレッド１の処理によって求められたことを表わす。適用分離行列規定区間７２よりも右側の数字も同様に、分離行列が何番目のスレッドに由来するかを表わす。
【００７８】
なお、学習を開始する時点において、他のスレッドで求まった分離行列が存在する場合は、それを学習の初期値として使う。これを「分離行列の引継ぎ」と呼ぶことにする。図７に示す例では、スレッド３の最初の学習が始まるタイミングである学習開始タイミング７５では、スレッド１に由来する分離行列７２が既に求まっているので、それを学習の初期値として使う。
【００７９】
このような処理を行うことで、スレッド間でパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）の発生を防止または低減することができる。スレッド間のパーミュテーションとは、例えば、１番目のスレッドで求まった分離行列では１番目のチャンネルに音声、２番目のチャンネルに音楽が出力されるのに対し、２番目のスレッドで求まった分離行列ではそれが逆転しているなどの問題である。
【００８０】
図７を参照して説明したように、他のスレッドで求まった分離行列が存在する場合は、それを学習の初期値として使う「分離行列の引継ぎ」を行うことで、スレッド間のパーミュテーションを低減することが可能となる。また、スレッド１の学習では分離行列が十分には収束していなくても、それを次のスレッドが引き継ぐことで、収束の度合いを向上させることができる。
【００８１】
このように時間をずらしながら複数のスレッドを起動することで、分離行列は、スレッド間のずれ、すなわちブロックシフト量７６とほぼ同じ間隔で更新されていくのである。
【００８２】
［Ｂ．従来技術の問題点］
これまでに説明した「Ａ３．ＩＣＡを用いた音源方向推定」と「Ａ４．ＩＣＡのリアルタイム化」とを組み合わせると、複数音源に対応したリアルタイムの音源方向推定を行なうことができる。
【００８３】
すなわち、図７を参照して説明すると、各スレッドから求まった分離行列（図７に示す区間７２，７３・・で使用する分離行列など）それぞれに対して、先に説明した、数式［５．１］〜［５．４］を適用することで、音源方向を求めることが可能となる。
【００８４】
しかし、単に両者を組み合わせるだけでは、いくつか問題が生じる。以下、次の５つの観点から見た問題点について説明する。
Ｂ１．遅延
Ｂ２．追従遅れ
Ｂ３．分離行列の特異性やゼロ要素
Ｂ４．目的の違い
Ｂ５．パラメーター設定の柔軟性
以下、それぞれの観点について説明する。
【００８５】
［Ｂ１．遅延］
図６を参照して説明した「ずらし適用」を用いると、音源分離の遅延はほぼ無くせるが、分離行列から方向を算出する限り、音源方向推定の遅延は依然として残る。すなわち、図６において、現在時刻の観測信号Ｘ（ｔ）６２に対応した分離結果は、Ｘ（ｔ）と分離行列６３とを乗じることで求まる。すなわち、図６の現在時刻の分離結果６４として、
Ｙ（ｔ）＝ＷＸ（ｔ）
この分離結果が求まる。
【００８６】
しかし、分離行列６３から算出される音源方向は、それより前の区間である学習データブロック６１に対応したものである。言い換えると、現在時刻の観測信号６２に対応した音源方向が求まるのは、現在時刻を含むブロック６６から新たな分離行列が学習された時点である。
【００８７】
なお、学習データブロック６１の区間と現在時刻とで音源方向が同一であると仮定し、分離行列６３から算出された音源方向をＹ（ｔ）の音源方向として流用すると、見かけ上は遅延が発生しなくなる（これを「時刻の付け替え」と呼ぶ）。しかし、学習データブロック６１の区間と現在時刻との間で音源が変化した場合（音源が移動したり急に鳴り始めたりした場合）は、分離行列と観測信号とで不一致が発生するため、不正確な音源方向が付与される。しかも、分離行列だけからでは、不一致が発生しているかどうかは判定できない。
【００８８】
また、「ずらし適用」を行なうと、遅延の発生と同じ理由により、音源方向と分離結果とが厳密には対応しなくなる。例えば、分離行列（６−３）から算出される音源方向に厳密に対応した分離結果は、分離行列の学習に用いた区間の学習データブロック６１と分離行列６３とを乗じることで生成されるものである。しかし、「ずらし適用」では、分離行列６３と学習データブロック６１の観測信号との乗算は行なわれない。分離行列６３と現在時刻の観測信号との積である現在時刻の分離結果６４も、分離結果の側で学習データブロック６１と同じ開始終了時刻を持つ分離結果スペクトログラム区間６５のデータも、いずれも分離行列６３から算出される音源方向とは厳密には一致していない。
【００８９】
［Ｂ２．追従遅れ］
「ずらし適用」を用いると、分離行列の学習に用いた区間（例えば学習データブロック６１）と現在時刻の観測信号６２との間で音源が変化した場合（音源が移動したり、急に鳴り始めたりした場合）に、一時的に不整合が生じる。その後、その変化した音源に対応した分離行列が学習によって求まるため、やがて不整合は消える。ここではその現象を「追従遅れ」と呼ぶ。追従遅れは、音源分離の性能には必ずしも悪影響を与えないが、その場合でも方向推定の性能には悪影響を与えることがある。以下では、その点について説明する。
【００９０】
図８は追従遅れの概念図である。音源の一つが、ある時刻（ｔ）を境に鳴り始めたとする。（例えば、ある人が急にしゃべり始めた場合など。）なお、その鳴り始めた時刻の前後において、他の音源は鳴りっぱなしでも無音でも構わない。その鳴り始めの前後の区間について、
学習の収束の度合いを（ａ）グラフ１に示す。
観測信号と分離行列との不整合の度合いを（ｂ）グラフ２に示す。
【００９１】
最初に（ａ）グラフ１について説明する。着目している音源が非アクティブ（無音）である非アクティブ区間がしばらく続いていれば、その区間の観測信号（例えば観測区間８１ａから学習された分離行列は十分収束している。この収束度合いは、（ａ）グラフ１の収束度合いデータ８５ａに対応する。
【００９２】
時刻（ｔ）において、音源が鳴り始めても、「ずらし適用」により、ある時間は「非アクティブ」で学習された分離行列が使われる。例えば、（ａ）グラフ１のデータ区間８５ｂでも、「非アクティブ」で学習された分離行列が使われる。８５ｂの長さは、学習時間（例えば図７における７７の時間長）よりは長く、ブロックシフト幅（例えば７６の時間長）と学習時間との合計より短いか同じ値であり、学習の開始するタイミングと音が鳴り始めるタイミングとの兼ね合いで決まる。
【００９３】
その後、「音源アクティブ」を含む観測区間（例えば区間８１ｂ）から学習された分離行列が使われるようになる。例えば、（ａ）グラフ１の収束度合いデータのデータ区間８５ｃでは、観測区間８１ｂから学習された分離行列が使われる。ただし、学習のループ回数を少なく設定している場合（例えば、１００回のループで収束するところを、学習時間を短くするために２０回で打ち切っているような場合）は、分離行列はまだ収束し切っていない。従って、（ａ）グラフ１のデータ区間８５ｃの収束度は低い。
【００９４】
その後、スレッド間で「分離行列の引継ぎ」（前述の「４．ＩＣＡのリアルタイム化」を参照）を繰り返すにつれて、分離行列が収束していく。（ａ）グラフ１の収束度合いデータのデータ区間８５ｃ→８５ｄ→８５ｅ→８５ｆである。データ区間８５ｆ以降は、分離行列は十分に収束しているとする。
【００９５】
次に、観測信号と分離行列との不整合の度合いを示した（ｂ）グラフ２について説明する。このグラフは、分離行列と現在時刻の観測信号との間で、不整合がどの程度発生しているかを表わす概念的なグラフであり、以下のどちらかが発生した場合に不整合の度合いが高くなる。
（１）分離行列の学習に用いられた区間と、現在時刻との間で、音源が変化した。
（２）分離行列が完全には収束していない。
【００９６】
不整合度合いデータのデータ区間８８ｂでは、分離行列は収束しているが、（１）の理由によって不整合度が高くなる。データ区間８８ｃ〜８８ｆでは、分離行列がだんだん収束していくため、不整合度も低下していく。
【００９７】
なお、不整合が生じても、音源分離に対しては必ずしも悪影響を与えない。例えば、主要な音源が一つしか鳴っていない、つまり、音源非アクティブ区間である観測区間８１で背景ノイズ程度しか存在しないような場合は、そもそも音同士が混ざっていないため、分離行列に不整合が発生しても音源分離には影響しない。また、観測信号を目的音と妨害音とに分離して目的音のみを使用する場合でも、妨害音が連続的に鳴っている音で目的音が断続的な音（急に鳴り始める音）であれば、目的音の分離結果には影響を与えない。具体的には、連続的な音は、他の音の有無にかかわらず、どれか一つのチャンネルのみに出力され続ける。一方で、他の音源が急に鳴り始めると、８８bの段階ではその音が全ての出力チャンネルに出力されるが、不整合度が下がるにつれて、その音は連続的な音のチャンネルとは別のどれか一つのチャンネルのみに出力されるようになる。したがって、妨害音が連続的である場合は、妨害音が抑圧されている出力チャンネルが目的音の鳴る前後を通じて存在しており、そのようなチャンネルを選択すれば、妨害音の除去は成功する。
【００９８】
一方、音源方向推定に対しては、分離行列と観測信号との不整合が悪影響を与える。例えば、図８中のデータ区間８８ｂの区間では、前述した「２．追従遅れ」の「時刻の付け替え」を行なった場合に、不正確な方向（音源アクティブ区間ではなくて音源非アクティブ区間から算出される方向）が付与されてしまう。また、データ区間８８ｃ〜８８ｄでは分離行列が未収束であるため、そこから計算される分離行列も不正確な値である可能性がある。
【００９９】
なお、音源が一つの場合に限定すれば、先に説明した［１．マイク間の位相差を利用した音源方向推定」では追従遅れなしで音源が算出できる。すなわち、音源が一つしかない環境で「リアルタイムＩＣＡ＋分離行列から音源方向推定」という方式を用いることは、「１．マイク間の位相差を利用した音源方向推定」よりも問題が多い。
【０１００】
［Ｂ３．分離行列の特異性やゼロ要素］
先に、［Ａ３．ＩＣＡを用いた音源方向推定］において説明した特許文献５（特許第３８８１３６７号公報）に記載の方式では分離行列の逆行列を用いる。そのため、分離行列が特異（ｓｉｎｇｕｌａｒ）に近い場合、すなわち分離行列の行列式（ｄｅｔｅｍｉｎａｎｔ）が０に近い値の場合に、逆行列が適切には求まらず、発散する可能性がある。また、逆行列が存在しても、行列の要素のどれかが０である場合、前述の式［５．４］において関数［ａｎｇｌｅ（）］が有効な値を返せなくなる。
【０１０１】
例えば、図７において、システム起動時には分離行列として単位行列が代入されているとする（すなわち図７の区間７１で単位行列が使われているとする）。単位行列には逆行列が存在する（単位行列自身）が、その逆行列は対角要素以外の値が０である。そのため、先に説明した式［５．４］では音源方向が求まらない。従って、分離行列の逆行列から音源方向を求める方式では、分離行列として単位行列が代入されている区間７１の区間の音源方向が求まらない。
【０１０２】
開始直後以外でも、分離行列が特異に近くなったり逆行列に０の要素が含まれたりすると、正しい音源方向が求まらなくなる。また、チャンネル数（マイク数）が多くなると、逆行列の計算量も増加する。そのため、逆行列なしで同様の操作が可能であるなら、その方が望ましい。
【０１０３】
［Ｂ４．目的の違い］
先に、［Ａ３．ＩＣＡを用いた音源方向推定］において説明した特許文献５（特許第３８８１３６７号公報）では、ＩＣＡの分離行列から方向を求めているが、その主な目的は周波数ビン間のパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）問題の解消である。すなわち、周波数ビンごとに音源方向を算出し、それからほぼ同一の音源方向を持つチャンネルを周波数ビン間でグルーピングすることで、周波数ビン間のパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）を揃えている。
【０１０４】
特許文献５（特許第３８８１３６７号公報）の後で、例えば、特許文献４［特開２００６−２３８４０９］において、パーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）問題を起こしにくい方式が紹介された。しーかし、この方式と特許第３８８１３６７号とを組み合わせても、方向推定の精度が向上するわけではない。
【０１０５】
［Ｂ５．パラメーター設定の柔軟性］
ＩＣＡの分離行列から音源方向を算出する場合、方向の更新頻度や対応する区間は分離行列の更新頻度や学習データ区間と一致してしまい、別々に設定することができない。すなわち、方向の更新間隔は図７に示すブロックシフト幅７６とほぼ一致し、そうして求まった方向はブロック内の観測信号（例えば図７のデータ区間７１）の平均的な値である。
【０１０６】
ブロックのシフト幅や蓄積区間長、さらに学習のループ回数（学習時間に影響する）などは、分離の性能やＣＰＵパワーなどを考慮して決める必要があるため、分離にとっては最適でも方向推定にとっては最適とは限らない。
【０１０７】
なお、上記のＩＣＡを用いた音源方向推定は、音源分離（すなわちＩＣＡ）を行なった後で方向推定を行なう方式である。他方、従来技術の中には、上記とは逆順、すなわち音源方向推定をした後にその値を用いて音源分離を行なうものも存在する。後者の方法では、チャンネル間の対応付けは容易であり、分離自体は低遅延で行なうことができる。しかし、「方向推定の後で音源分離」という方式に共通する課題として、方向推定の精度が分離精度に影響するという点が挙げられる。さらに、たとえ方向自体は正しく推定できても、音源分離の段階において他の要因によって分離精度が低下する可能性がある。
【０１０８】
方向推定の後に周波数マスキングに基づく音源分離を行なっている例としては、特開2004-325284「音源方向を推定する方法、そのためのシステム、および音源の分離方法、そのためのシステム」などが挙げられる。この方式では、周波数ビンごとに方向を推定した後、目的音が支配的な周波数ビンの信号を残すことで目的音を抽出している。しかしこの方法では周波数領域でのスパース性（同一の周波数では、目的音と妨害音とが混合することは稀である）を仮定しているため、その仮定が成り立たない信号の混合に対しては、方向推定の精度も音源分離の精度も低下する。
【０１０９】
方向推定の後にビームフォーミングに基づく音源分離を行なっている例としては、特開2008-64892「音声認識方法およびそれを用いた音声認識装置」などが挙げられる。この方式では、MUSIC法（相関行列を固有値分解することで信号部分空間と雑音部分空間を求め、任意の音源位置ベクトルと雑音部分空間の内積の逆数を求めることにより、音源の音波到来方向や位置を調べる手法）を用いて遠距離音源の到来方向や近距離音源の位置を推定した後、分散最小ビームフォーマー（minimum variance beamformer または Frost beamformer）を用いて目的音を抽出している。分散最小ビームフォーマーでは目的音の音源から各マイクまでの伝達関数が既知である必要があり、特開2008-64892では音源の方向や位置から伝達関数を推定している。しかし、伝達関数の推定において、残響のないモデルを想定しており、また、音速やマイク位置などの値を用いているため、それらの想定や値が実際の環境と食い違っている場合に、分離の精度が低下する。
【０１１０】
以上をまとめると、従来、音源方向推定の方式はいくつか存在していたが、
「マイク間の位相差を利用した音源方向推定」が持つ以下の利点：
・追従遅れも遅延も少ない
および、「ＩＣＡを用いた音源方向推定」が持つ以下の利点：
・ブラインド分離が可能
・多音源に対応
・分離結果と音源方向との対応付けが容易
これらの両方を併せ持つ方式は存在しなかった。
【０１１１】
さらに、ＩＣＡを用いた音源方向推定において、分離結果に周波数ビン間のパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）がないことを前提とし、周波数ビン間の関係を用いて方向推定の精度を上げるような方式も存在していなかった。
【先行技術文献】
【特許文献】
【０１１２】
【特許文献１】特開２００５−７７２０５号公報
【特許文献２】特開２００５−４９１５３号公報
【特許文献３】特許第３０６９６６３号公報
【特許文献４】特開２００６−２３８４０９号公報
【特許文献５】特許第３８８１３６７号公報
【非特許文献】
【０１１３】
【非特許文献１】『入門・独立成分分析』（村田昇著、東京電機大学出版局）
【非特許文献２】『詳解  独立成分分析』の「１９．２．４．フーリエ変換法」
【発明の概要】
【発明が解決しようとする課題】
【０１１４】
本発明は、このような状況に鑑みてなされたものであり、独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ；ＩＣＡ）を用いて信号毎の高精度な分離処理とリアルタイムでの音源方向推定を実現する信号処理装置、および信号処理方法、並びにプログラムを提供することを目的とする。
【０１１５】
本発明は、例えば、時間周波数領域において、観測信号と分離結果との間で相互共分散行列を計算し、その共分散行列の要素同士の位相差から音源方向を算出する。また、位相差から音源方向を算出する際に、位相差と周波数との比例係数を精度良く推定する。本発明の処理構成は、従来法である、分離行列の逆行列から音源方向を推定する方法と比べて、以下のような利点がある。（１）分離行列法よりも短い遅延で計算可能。（２）方向は１フレーム毎の頻度で算出可能。（３）分離行列の逆行列が存在しない場合や、逆行列の要素に０が含まれる場合にも方向が算出可能。これらの利点を有する。
【課題を解決するための手段】
【０１１６】
本発明の第１の側面は、
複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求める学習処理部と、
前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成する分離処理部と、
前記分離処理部において生成した分離信号の音源方向を算出する音源方向推定部を有し、
前記音源方向推定部は、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出する構成である信号処理装置にある。
【０１１７】
さらに、本発明の信号分離装置の一実施態様において、前記分離処理部は、前記複数音源からの出力を混合した混合信号に対して短時間フーリエ変換（ＳＴＦＴ）を実行して時間周波数領域の観測信号スペクトログラムを生成し、さらに、前記観測信号スペクトログラムに対して、前記分離行列を適用して時間周波数領域の分離信号スペクトログラムを生成する構成であり、前記音源方向推定部は、時間周波数領域において対応する時間区間の観測信号スペクトログラムと分離信号スペクトログラムとの間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出する構成である。
【０１１８】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、時間周波数領域の観測信号の異なる時間の観測信号間の共分散行列を算出し、算出した共分散行列と、前記学習処理部の算出した分離行列に対するエルミート転置行列を適用して、観測信号と分離信号との間の相互共分散行列を算出する構成である。
【０１１９】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、前記位相差の算出処理において、各音源（チャンネル）について、異なる位置に設定された信号入力源であるマイクペアに対応する位相差を個別に算出する構成である。
【０１２０】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、相互共分散行列の要素同士の位相差情報を利用して、位相角と周波数との比例定数（直線の傾き）を算出する直線当てはめ処理を実行して、前記比例定数を利用して各分離信号に対応する音源方向を算出する。
【０１２１】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、前記直線当てはめ処理として、低い周波数から順に、位相角と周波数との対応情報のアンラップ処理を実行して前記比例定数の更新処理を行う再帰法を実行して比例定数（直線の傾き）を算出する。
【０１２２】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、前記直線当てはめ処理として、前記比例定数（直線の傾き）の値を予め設定した最小値から最大値まで少しずつ変化させて、前記位相差情報との誤差を計算し、誤差の少ない比例定数を選択するリニアサーチ法を実行して比例定数（直線の傾き）を算出する。
【０１２３】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、前記直線当てはめ処理として、前記比例定数（直線の傾き）の値を予め設定した最小値から最大値まで少しずつ変化させて、前記位相差情報との誤差を計算し、誤差の少ない比例定数を選択するリニアサーチ法を実行し、さらに、前記リニアサーチ法によって得られた比例定数を初期値として、低い周波数から順に、位相角と周波数との対応情報のアンラップ処理を実行して前記比例定数の更新処理を行う再帰法を実行して比例定数（直線の傾き）を算出する。
【０１２４】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、各音源（チャンネル）について、異なる位置に設定された信号入力源であるマイクペア各々に対応する複数の比例定数（直線の傾き）を算出し、算出した複数の比例定数（直線の傾き）から当てはめ誤差の少ないデータを選択し、選択データの平均値を算出して最終的な比例定数（直線の傾き）を算出する。
【０１２５】
さらに、本発明の信号分離装置の一実施態様において、前記信号処理装置は、さらに、前記直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して、前記直線からのずれの大きい位相角データを持つ信号分離結果の周波数成分を減衰させる周波数マスキングを実行する構成である。
【０１２６】
さらに、本発明の信号分離装置の一実施態様において、前記信号処理装置は、さらに、前記直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して、前記分離行列を生成し、生成した分離行列を前記学習処理部において利用する構成である。
【０１２７】
さらに、本発明の信号分離装置の一実施態様において、前記音源方向推定部は、各マイクペアを結ぶベクトル［ｐｉｍ］について大きさを１に正規化して横方向に並べて生成した行列［Ｐｋ］と、各マイクペアについての経路差［Ｄｉｍｋ］をマイク間隔で正規化して横方向に並べて生成したベクトル［Ｄｋ］を算出し、算出した［Ｐｋ］とベクトル［Ｄｋ］を適用して、音源方向の角度、または、マイクペアから見た音源方向を示すベクトルを算出する構成である。
【０１２８】
さらに、本発明の第２の側面は、
信号処理装置において実行する信号処理方法であり、
学習処理部が、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求める学習処理ステップと、
分離処理部が、前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成する分離処理ステップと、
音源方向推定部が、前記分離処理ステップにおいて生成した分離信号の音源方向を算出する音源方向推定ステップを有し、
前記音源方向推定ステップは、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出するステップである信号処理方法にある。
【０１２９】
さらに、本発明の第３の側面は、
信号処理装置において信号処理を実行させるプログラムであり、
学習処理部に、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、前記混合信号を分離する分離行列を求めさせる学習処理ステップと、
分離処理部に、前記観測信号に対して、前記分離行列を適用して前記混合信号を分離して各音源対応の分離信号を生成させる分離処理ステップと、
音源方向推定部に、前記分離処理ステップにおいて生成した分離信号の音源方向を算出させる音源方向推定ステップを有し、
前記音源方向推定ステップは、
時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出して、算出した位相差を適用して各分離信号に対応する音源方向を算出させるステップであるプログラムにある。
【０１３０】
なお、本発明のプログラムは、例えば、様々なプログラム・コードを実行可能な汎用システムに対して、コンピュータ可読な形式で提供する記憶媒体、通信媒体によって提供可能なプログラムである。このようなプログラムをコンピュータ可読な形式で提供することにより、コンピュータ・システム上でプログラムに応じた処理が実現される。
【０１３１】
本発明のさらに他の目的、特徴や利点は、後述する本発明の実施例や添付する図面に基づくより詳細な説明によって明らかになるであろう。なお、本明細書においてシステムとは、複数の装置の論理的集合構成であり、各構成の装置が同一筐体内にあるものには限らない。
【発明の効果】
【０１３２】
本発明の一実施例の構成によれば、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、混合信号を分離する分離行列を求め、観測信号に対して分離行列を適用して分離信号を生成するとともに、分離信号の音源方向を算出する音源方向推定処理を実行する。音源方向推定処理においては、時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出し、算出した位相差を適用して各分離信号に対応する音源方向を算出する。具体的には、相互共分散行列の要素同士の位相差情報を利用して位相角と周波数との比例定数（直線の傾き）を算出し、比例定数を利用して各分離信号に対応する音源方向を算出する。この処理によって、高精度な信号分離および音源方向推定処理が実現される。
【図面の簡単な説明】
【０１３３】
【図１】基本的な音源方向の推定方式について説明する図である。
【図２】Ｎ個の音源から異なる音が鳴っていて、それらをｎ個のマイクで観測するという状況について説明する図である。
【図３】周波数ビンにおける分離（図３（Ａ）参照）と、全周波数ビンの分離処理（図３（Ｂ）参照）について説明する図である。
【図４】分離行列Ｗ（ω）の逆行列としての推定伝達関数行列について説明する図である。
【図５】観測信号のスペクトログラムを重複のある複数のブロック１〜Ｎに分割し、ブロックごとに学習を行なって分離行列を求める処理例について説明する図である。
【図６】各ブロックから求まった分離行列を、それより以降の観測信号へ適用する「ずらし適用」について説明する図である。
【図７】重複のあるブロックから分離行列を求めるためスレッドと呼ばれる処理の単位を、時刻をずらして複数並列に起動する方式について説明する図である。
【図８】スレッド単位での並列処理において発生する追従遅れについて説明する図である。
【図９】現在時刻を含む区間における相互共分散について説明する図である。
【図１０】相互共分散行列のｋ列目のｉ行目の要素とｍ行目の要素との間で比をとり、さらに大きさを１に正規化したデータｒｉｍｋ（ω）の位相角ａｎｇｌｅ（ｒｉｍｋ（ω））をプロットしたデータ、および直線近似処理から得られる直線の傾き［Ｋｉｍｋ］の取得とマイクペアごとの音源方向の算出処理について説明する図である。
【図１１】ｒｉｍｋ（ω）の実部と虚部とをプロットしたデータを利用した直線当てはめ処理について説明する図である。
【図１２】本発明の一実施例に従った信号処理装置の構成例について説明する図である。
【図１３】相互共分散行列を観測信号と分離結果とから計算して音源方向の推定処理を実行する音源方向推定部の構成例について説明する図である。
【図１４】観測信号と分離行列とから音源方向の推定処理を実行する音源方向推定部の構成例について説明する図である。
【図１５】学習処理部のスレッド制御部の詳細構成例について説明する図である。
【図１６】スレッド演算部において実行する処理について説明する図である。
【図１７】学習スレッドの状態遷移について説明する図である。
【図１８】学習スレッドの状態遷移について説明する図である。
【図１９】リアルタイム分離・音源方向推定処理について説明するフローチャートを示す図である。
【図２０】短時間フーリエ変換の詳細について説明する図である。
【図２１】図１９に示すフローチャートにおけるステップＳ１０１の初期化処理の詳細について説明するフローチャートを示す図である。
【図２２】スレッド制御部１１２による複数の学習スレッド１，２に対する制御シーケンスを示す図である。
【図２３】図１９に示すフローチャートにおけるステップＳ１０５においてスレッド制御部１１２によって実行されるスレッド制御処理について説明するフローチャートを示す図である。
【図２４】図２３に示すフローチャートにおけるステップＳ２０３において実行される待機中状態における処理について説明するフローチャートを示す図である。
【図２５】図２３に示すフローチャートのステップＳ２０４において実行される蓄積中状態における処理について説明するフローチャートを示す図である。
【図２６】図２３に示すフローチャートのステップＳ２０５において実行される学習中状態における処理について説明するフローチャートを示す図である。
【図２７】図２６に示すフローチャートのステップＳ２３９において実行される、分離行列更新処理について説明するフローチャートを示す図である。
【図２８】図２６に示すフローチャートのステップＳ２４１において実行される待機時間の設定処理について説明するフローチャートを示す図である。
【図２９】学習スレッドの処理について説明するフローチャートを示す図である。
【図３０】図２９に示すフローチャートのステップＳ２９４において実行される、コマンド処理について説明するフローチャートを示す図である。
【図３１】図３０に示すフローチャートのステップＳ３０５において実行される処理の一例である分離行列の学習処理例について説明するフローチャートを示す図である。
【図３２】図１９に示すフローチャートのステップＳ１０９において実行する音源方向推定処理について説明するフローチャートを示す図である。
【図３３】図１９に示すフローチャートのステップＳ１０９において実行する音源方向推定処理について説明するフローチャートを示す図である。
【図３４】図３２のフローチャートにおけるステップＳ３２４、図３３のフローチャートにおけるステップＳ３３４の「位相差の計算」の詳細について説明するフローチャートを示す図である。
【図３５】図３４に示すフローの処理結果として生成されるデータの例について説明する図である。
【図３６】図３２のフローチャートにおけるステップＳ３２５、図３３のフローチャートにおけるステップＳ３３５の「直線当てはめ」の詳細について説明するフローチャートを示す図である。
【図３７】再帰法によって、（直線）当てはめ処理を行うシーケンスを説明するフローチャートを示す図である。
【図３８】再帰法による（直線）当てはめ処理について説明する図である。
【図３９】リニアサーチによって、（直線）当てはめ処理を行うシーケンスを説明するフローチャートを示す図である。
【図４０】リニアサーチによる（直線）当てはめ処理について説明する図である。
【図４１】図３５のプロットデータを利用して直線当てはめを行った結果の例を示す図である。
【図４２】図３２に示すフローチャートにおけるステップＳ３２６の音源方向の計算、図３３に示すフローチャートにおけるステップＳ３３６の音源方向の計算処理について説明するフローチャートを示す図である。
【図４３】マイクの配置が一直線状以外である場合のマイクペアの角度の分だけ式［８．７］の値を補正する処理例について説明する図である。
【図４４】あるチャンネルのあるマイクペアにおける直線当てはめの結果の例について説明する図である。
【図４５】位相差ψの絶対値に応じて分離結果の周波数成分を減衰させることにより、分離の精度を向上させるために適用する関数の例について説明する図である。
【図４６】マイクペアの一つと音源方向ベクトルとの関係を表わした図である。
【図４７】音源方向ベクトル推定アルゴリズムの概要を説明するフローチャートを示す図である。
【図４８】各マイクペアのベクトルｐｉｍについて、大きさを１に正規化してから横方向に並べて生成した行列［Ｐｋ］、各マイクペアについて、経路差Ｄｉｍｋをマイク間隔で正規化してから横方向に並べて生成したベクトル［Ｄｋ］、これらの生成処理の詳細について説明するフローチャートを示す図である。
【図４９】図４７に示すフローチャートにおけるステップＳ５０６の「ｑｋの推定」について説明するフローチャートを示す図である。
【図５０】音源の存在する方向がある程度限定できる場合に、それを初期値に反映する処理について説明する図である。
【発明を実施するための形態】
【０１３４】
以下、図面を参照しながら本発明の信号処理装置、および信号処理方法、並びにプログラムの詳細について説明する。
【０１３５】
本発明は、複数の信号が混合された信号を独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ；ＩＣＡ）を用いて分離し、音源方向を推定する。
【０１３６】
上述したように、従来、音源方向推定の方式はいくつか存在していたが、
「マイク間の位相差を利用した音源方向推定処理」が持つ「追従遅れも遅延も少ない」という利点と、「ＩＣＡを用いた音源方向推定処理」が持つ「ブラインド分離が可能」、「多音源に対応」、「分離結果と音源方向との対応付けが容易」という利点、これらのすべての利点を併せ持つ方式は存在しない。また、「ＩＣＡを用いた音源方向推定処理」において、分離結果に周波数ビン間のパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）がないことを前提とし、周波数ビン間の関係を用いて方向推定の精度を上げるような方式も存在していないという問題がある。
【０１３７】
上記問題点を解決するため、本発明では以下の２点を新たに導入する。
１．相互分散行列
現在時刻を含む区間において、観測信号と分離結果との間で相互共分散行列を計算し、その相互共分散行列を分離行列の代わりに用いて音源方向を計算する。
２．直線当てはめ
音源方向の計算の際に、全周波数における位相差のプロットを直線で近似し、その直線の傾きから音源方向を計算する。
まず、
１．相互分散行列
２．直線当てはめ
これらについて説明する。
【０１３８】
［１．相互共分散行列］
現在時刻を含む区間における相互共分散について、図９を用いて説明する。この図９は、先に説明した図６に本発明の新規部分を追加した図である。
【０１３９】
（ａ）観測信号と（ｂ）分離結果との双方において、現在時刻を含む区間を設定する。例えば、現在のフレーム番号を［ｔ］とし、区間の開始時刻を［ｔ−Ｌ＋１］番目のフレーム、終了時刻を［ｔ］番目のフレームとする。現在時刻を含む区間のデータは、ａ）観測信号側の観測信号データ９２、（ｂ）分離結果側の分離結果データ９４である。
【０１４０】
この区間の観測信号と分離結果との間で、周波数ビンごとに、以下に示す式［６．１］で表わされる相互共分散行列ΣＸＹ（ω）を計算する。
【０１４１】
【数１４】
2010020294.tif 000016
【０１４２】
この式において、Ｅｔ［  ］は区間内での平均を表わし、具体的には式［６．２］で計算される。μＸ（ω）およびμＹ（ω）は、それぞれ区間内におけるＸ（ω，ｔ）およびＹ（ω，ｔ）の平均であり、式［６．３］および式［６．４］で計算される。なお、μＸ（ω）およびμＹ（ω）を引く操作は、信号の平均が０に近いと期待できる場合には省略することができる。
【０１４３】
このように、本発明の信号処理装置における音源方向推定部は、時間周波数領域において対応する時間区間の観測信号スペクトログラムと分離信号スペクトログラムとの間で相互共分散行列を計算する。
【０１４４】
なお、観測信号と分離結果との相互共分散行列を直接計算する代わりに、観測信号自体の共分散行列と分離行列とから計算することも可能である。すなわち、観測信号の共分散行列を式［６．６］で計算し（必要に応じてＸ（ω，ｔ）から平均μＸ（ω）を引く）、その右側に分離行列のエルミート転置Ｗ（ω）Ｈを乗じる（式［６．５］）ことによっても、ほぼ同じ値を計算することができる。（区間内で分離行列Ｗ（ω）が一定の場合は、式［６．５］の両辺は同一の値となる。）
【０１４５】
このように計算された相互共分散行列ΣＸＹ（ω）は、分離行列Ｗ（ω）と同様の使い方をすることで音源方向を求めることができる。また、相互共分散行列ΣＸＹ（ω）から推定される音源方向は、分離行列によるものと先に説明した「Ａ１．マイク間の位相差を利用した音源方向推定」によるものとを兼ね備えた特徴を持っている。以下で、その点について説明する。
【０１４６】
先に説明した分離結果の算出式［３．５］、すなわち、
Ｙ（ω，ｔ）＝Ｗ（ω）Ｘ（ω，ｔ）
の逆演算である
Ｘ（ω，ｔ）＝Ｗ−１（ω）Ｙ（ω，ｔ）
という関係を用いると、式［６．５］の右辺は式［６．８］のように変形できる。すなわち、分離行列の逆行列と、分離結果自身の共分散行列とを乗じた形である。なお、式［６．８］は説明のための変形であり、実際の処理において逆行列を計算する必要はない。
【０１４７】
分離が完全に行なわれている場合、すなわち分離行列と観測信号との間に不一致が発生せず、さらに分離行列が十分収束している場合であれば、分離結果Ｙ（ω，ｔ）の各要素は統計的に独立であるため、分離結果の共分散行列は対角要素のみ正の値（非対角要素は０）を持つ対角行列となる。対角行列を乗じても複素数の要素比の位相とは無関係のため、分離が完全のときの式［６．８］は、音源方向については分離行列の逆行列と等価である。
【０１４８】
従って、相互共分散行列ΣＸＹ（ω）内の要素を、以下に示す式［７．１］のように表わすと、音源方向は式［７．２］で推定される。（式［５．４］と同様である。）すなわち、逆行列を用いなくても、音源方向の推定は可能である。
【０１４９】
【数１５】
2010020294.tif 000017
【０１５０】
一方、分離が不十分な場合でも、同時に鳴っている音源が一つであれば、音源方向は推定可能である。式［６．５］によれば、前述の相互共分散行列は観測信号の共分散行列と分離行列のエルミート転置との積であると考えることもできる。ここで、分離行列が単位行列である場合を考える。これは、信号が最も分離されていない状態（観測信号がそのまま出力される）に相当する。
【０１５１】
この場合、前述の相互共分散行列は観測信号の共分散行列と等価となるため、式［７．２］から算出される音源方向は、式［１．６］と式［１．４］とから算出される音源方向（「Ａ１．マイク間の位相差を利用した音源方向推定」）と同じ値である。すなわち、観測信号の共分散行列ΣＸＸ（ω）はｎ×ｎの要素を持ち、各要素は式［６．７］のように表わせる（必要に応じて、平均を引く）。ΣＸＸ（ω）はｎ個の列を持つが、どの列を用いても、式［１．６］を用いた場合とほぼ同一の値が算出される。
【０１５２】
以上、２つの場合をまとめる。分離が完全な場合は、相互共分散行列から算出される音源方向は分離行列の逆行列から算出される音源方向と同じ値をとる。もう一方の極端な場合として、分離が全く行なわれてない場合は、相互共分散行列から算出される音源方向はマイク間の位相差から算出される音源方向とほぼ同じ値をとる。つまり、両者の特徴を併せ持っているのである。
【０１５３】
以下に、従来のＩＣＡを使用した場合の音源方向推定処理と、相互共分散行列を用いた本発明の手法による音源方向推定処理の比較表を示す。
【０１５４】
【表１】
2010020294.tif 000018
【０１５５】
上記表に示すように、従来のＩＣＡを使用した場合の音源方向推定処理では、分離処理が不十分である場合は、音源方向の推定精度は低下してしまうが、相互共分散行列を用いた本発明の手法による音源方向推定処理では、単一音源の場合に良好な推定が実現される。
【０１５６】
なお、同時に鳴っている音源が複数で、かつ分離が不十分な場合は、本発明の相互共分散行列を用いた手法を適用しても正確な音源方向は推定できない。しかし、後述の「直線当てはめ」を用いることにより、複数音源が存在し、かつ分離処理が不十分な場合が発生したことを、高い精度で検出することができるため、方向推定以降の処理への影響を最小限にすることができる。
【０１５７】
さらに、相互共分散を用いた方式は、分離行列を用いた場合と比べて以下の利点を持つ。
・高頻度の更新（フレーム毎に更新可能）
・低遅延（１フレーム以内）
・追従遅れが少ない（または検出可能）
【０１５８】
共分散行列（または相互共分散行列）の計算は、フレーム番号ｔが増加するごとに上記の相互共分散行列ΣＸＹ（ω）の算出式［６．１０］または式［６．１１］を計算すればよく（必要に応じて、Ｘ（ω，ｔ）およびＹ（ω，ｔ）から平均を引く）、これは逆行列の計算と比べて計算量が少ない。また、図９に示す分離結果スペクトログラムに設定される区間の長さＬ（図９参照）を大きくしても、計算量は一定である。従って、本発明の方法では、音源方向をフレーム毎の頻度でかつ低遅延で算出することができる。
【０１５９】
また、音源方向の算出に用いる区間（図９に示す観測信号スペクトログラムに設定する区間９２、および分離結果スペクトログラムに設定する区間９４）は、分離行列の学習で用いる区間（例えば、学習データブロック区間９１）とは別個に設定することができる。従って、追従性と安定性とを兼ね備えたちょうど良い長さを設定すればよい。
【０１６０】
なお、上記の説明では（相互）共分散行列を計算していたが、共分散行列そのものの代わりに、各要素に適切な正の値を乗じた行列を用いて音源方向を計算しても構わない。（正の実数を乗じても、位相角には影響を与えないため。要素ごとに異なる値を乗じても同様である。）
【０１６１】
例えば、共分散行列を算出する式である前記した式［６．２］，［６．１０］，［６．１１］において、係数１／（Ｌ−１）を省略しても、音源方向は同一の値が算出される。また、共分散行列の代わりに相関行列（式［６．１２］、ただし、σｉｋ（ω）は共分散行列の各要素であり、式［７．１］で表わされる）を用いても、同一の音源方向が算出される。
【０１６２】
さらに、乗じる数が負の値や複素数であっても、共分散行列の各列の中で同一の値を乗じていれば、音源方向は同一の値が算出される（要素の比を計算する際に相殺されるため）。例えば、相互共分散行列ΣＸＹ（ω）の代わりに、式［６．１３］で表わされるΣ'ＸＹ（ω）を用いても、同一の音源方向が算出される。（ただし、式［６．１３］の右辺の第２項は対角要素以外が０の行列であり、また、対角要素ｈ１，...，ｈｎは任意の複素数を表わす。）
【０１６３】
この性質により、相互共分散行列の計算で用いる分離結果Ｙ（ω，ｔ）や分離行列Ｗ（ω）については、リスケーリング後のものの代わりに、リスケーリング前のものを用いても、同一の音源方向が算出される。（リスケーリングとは、分離結果の周波数ビンごとの大きさを調整する処理であり、詳細は後述する。）
【０１６４】
２．直線当てはめ
前記の式［７．２］を用いることで、周波数ビンごとかつマイクペアごとの音源方向を推定することができ、それらを平均することで音源方向を求めることができる。しかし、周波数ビンごとに算出した値を平均する方法には、以下の問題点がある。
【０１６５】
（問題１）ある原信号が周波数ビンωの成分を含まない場合（例えば原信号が帯域制限されている場合）、そのωにおける音源方向は不正確な値となり、平均を取ったときに悪影響を与える。
（問題２）周波数ビンωに対応した波長がマイク間距離の２倍よりも短くなると（すなわち、ωが、ω≧１＋Ｃ／（２ｄＦ）を満たす場合は）、音源方向が一意に決まらなくなる。具体的には、真の方向以外にも前記の式［１．３］を満たすθが存在するようになるため、式［１．３］や式［５．４］において、関数ａｓｉｎがどの値を返すべきか決定できなくなる。（空間のエリアシングと呼ばれる。）そのため、高い周波数の成分を音源方向に反映させることが難しい。
（問題３）周波数ビン間の関係などを考慮していないため、パーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）の起きにくい方式のＩＣＡを用いても、方向推定の精度が向上するわけではない。
【０１６６】
このような問題を解決するため、本発明では「直線当てはめ」という手法を用いる。これは、位相差との周波数ビン番号ωとのプロットを直線で近似し、その直線の傾きから音源方向を推定する方式である。以下、原理と方法について説明する。
【０１６７】
前記した相互共分散行列を示す式［７．１］の右辺において、相互共分散行列のｋ列目に着目する。そして、ｉ行目の要素とｍ行目の要素との間で比をとり、さらに複素数の大きさを１に正規化したものをｒｉｍｋ（ω）とする。ｒｉｍｋ（ω）は以下に示す式［８．１］によって示される。
【０１６８】
【数１６】
2010020294.tif 000019
【０１６９】
ただし、ｓｉｇｎ（ｚ）は複素数ｚの大きさを１に正規化する操作を表わし、式［８．２］で定義される。また、ｓｉｇｎ（）については、ａ／ｂとａ×ｃｏｎｊ（ｂ）とは等価である（式［８．１］参照）ため、以降の「ａとｂとの要素比」という表現は両方を表わすものとする。（複素数の掛け算・割り算とｓｉｇｎ（）との順序は、どちらを先に行なっても同じ値が算出される。例えば、式［８．１］の中央の式は、分子と分母についてｓｉｇｎをとってから比を計算しても、比をとってからｓｉｇｎを計算しても、同一の値となる。）
【０１７０】
前述した数式である数式［１．１］〜［１．７］、および数式［６．１］〜［６．１３］の考察により、［ｒｉｍｋ（ω）］は上記の式［８．３］で近似できる。すなわち、ｒｉｍｋ（ω）の位相角は［ω−１］と比例関係にあり、その比例定数Ｋｉｍｋと音源方向との間には式［８．４］の関係がある。言い換えれば、ｒｉｍｋ（ω）の位相角ａｎｇｌｅ（ｒｉｍｋ（ω））をプロットしたものを直線で近似すると、その直線の傾きはＫｉｍｋであり、Ｋｉｍｋが求まればマイクペアごとの音源方向も求まる。
【０１７１】
図１０に具体例を示す。図１０は、相互共分散行列のｋ列目のｉ行目の要素とｍ行目の要素との間で比をとり、さらに大きさを１に正規化したデータｒｉｍｋ（ω）の位相角ａｎｇｌｅ（ｒｉｍｋ（ω））をプロットしたデータ、および直線近似処理から得られる直線の傾き［Ｋｉｍｋ］の取得とマイクペアごとの音源方向の算出処理について説明する図である。
【０１７２】
図１０（ａ）は、実データから計算したｒｉｍｋ（ω）について、位相角をプロットしたものであり、縦軸が周波数ビンωを、横軸が位相角を表わす。点がほぼ直線状に並んでいることが分かる。ただし、位相角ａｎｇｌｅ（）は２πの範囲（この図では−π〜π）でしか値を特定できないため、−πとπとの間で値の「飛び」が存在する。
【０１７３】
図１０（ａ）に対して、「飛び」を考慮しながら直線当てはめを行なった結果を図１０（ｂ）に示す。直線当てはめの際には、「飛び」を最小にするため、周波数ビンごとに位相角の値に２Ｎπ（Ｎは適切な整数）を加えるアンラップ操作を実行した。図１０（ｂ）に示すアンラップ操作処理領域１０３に対してこの処理を実行した。この操作は「アンラップ（ｕｎｗｒａｐ）」と呼ぶ。このアンラップ操作の後、点に対応する最も当てはまる直線、すなわち図１０に示す当てはめ直線１０４を求める。なお、この処理の具体的なアルゴリズムについては後述する。
【０１７４】
なお、比例係数Ｋｉｍｋは、ｒｉｍｋ（ω）の位相角を用いる代わりに、複素数ｒｉｍｋ（ω）を直接用いても推定可能である。その方法について図１１を参照して説明する。なお、以降ではこの方法も「直線当てはめ」と呼ぶことにする。
【０１７５】
図１１（ａ）および（ｂ）は、それぞれｒｉｍｋ（ω）の実部と虚部とをプロットしたデータである。それぞれ、コサインとサインに近い曲線を描いているのが分かる。このプロットを上記の式［８．３］のｅｘｐ（ｊＫｉｍｌ（ω−１））で近似すると、実部・虚部共に最もよく当てはまる曲線（図１１（ｃ）に示す曲線１０５および（ｄ）に示す曲線１０６）に対応したＫｉｍｋが求める比例定数である。例えば、ｒｉｍｋ（ω）とｅｘｐ（ｊＫｉｍｌ（ω−１））との誤差を式［８．５］で定義し、この誤差を最小にするＫｉｍｋを求めればよい。（式［８．６］参照。ａｒｇ  ｍｉｎは、カッコ内を最小にするＫｉｍｋを求める操作を表わす。）
【０１７６】
このような直線当てはめを行なうことにより、前述の３つの問題は回避できる。
すなわち、
（問題１）ある原信号が周波数ビンωの成分を含まない場合（例えば原信号が帯域制限されている場合）、そのωにおける音源方向は不正確な値となり、平均を取ったときに悪影響を与える。
この問題に対して、ｒｉｍｋ（ω）の位相角が、いくつかの周波数ビンにおいて直線上から外れていても、大部分の周波数ビンにおいてほぼ直線上に位置していれば、傾きＫｉｍｋはほぼ正しく推定される。
【０１７７】
（問題２）周波数ビンωに対応した波長がマイク間距離の２倍よりも短くなると（すなわち、ωが、ω≧１＋Ｃ／（２ｄＦ）を満たす場合は）、音源方向が一意に決まらなくなる。具体的には、真の方向以外にも前記の式［１．３］を満たすθが存在するようになるため、式［１．３］や式［５．４］において、関数ａｓｉｎがどの値を返すべきか決定できなくなる。（空間のエリアシングと呼ばれる。）そのため、高い周波数の成分を音源方向に反映させることが難しい。
この問題に対して、直線当てはめの際にアンラップを行なうことで、空間のエリアシングは解消する。そのため、高い周波数に対応したビンも、音源方向に反映させることができる。
【０１７８】
（問題３）周波数ビン間の関係などを考慮していないため、パーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）の起きにくい方式のＩＣＡを用いても、方向推定の精度が向上するわけではない。
この問題に対して、「分離結果にパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）が発生していなければ、位相角と周波数ビンとは比例する」という関係を利用しているため、周波数ビンごとに方向を求める方式に比べて推定精度の向上が期待できる。言い換えると、特開２００６−２３８４０９より以前の方式では、分離結果に周波数ビン間のパーミュテーションが発生するため、図１０(a)のプロットは直線状にはならず、「直線当てはめ」を用いることができなかった。
このように、直線当てはめを行なうことにより、これらの３つの問題は回避できる。
【０１７９】
また、当てはめ誤差（式［８．５］）を用いることで、推定した音源方向がどの程度正確かを判定することができる。例えば、前述の追従遅れ（図８を参照して説明したデータ８８ｂ〜８８ｅ）が発生していて、かつその時点で鳴っている音源が複数ある場合、位相角をプロットしても直線上には並ばない。そのときに直線当てはめを強制的に行なうと、当てはめ誤差は大きな値となる。そこで、式［８．５］がある値（例えば１）よりも大きい場合は推定に失敗しているとみなし、推定した音源方向を棄却するといったことが可能となる。
【０１８０】
［実施例］
次に、上述した１．相互分散行列、２．直線当てはめを適用した音源方向推定処理を実行する本発明に従った信号処理装置の具体的な構成及び処理例について説明する。
【０１８１】
（装置の構成）
本発明の一実施例に従った信号処理装置の構成を図１２に示す。なお、図１２に示す装置構成は、は、本出願人が先に出願した特許出願２００６−３３１８２３「リアルタイム音源分離装置および方法」をベースとした構成である。この特許出願２００６−３３１８２３において示した構成に、音源方向推定部１２５を追加したものである。なお、図１２に示す信号処理装置は、具体的には例えばＰＣによって実現可能である。すなわち、分離処理部１２３、学習処理部１３０、音源方向推定部、後段処理部１２４、これらの各処理部の処理は、例えば、予め規定されたプログラムに従った処理を実行するＣＰＵによって実行可能である。
【０１８２】
図１２の中央に示す分離処理部１２３は主に観測信号の分離を行なう。図１２の右側に示す学習処理部１１０は、分離行列の学習を行なう。なお、分離処理部１２３における処理と、学習処理部１１０の処理は並列に行われる。分離処理部１２３の処理は、表（ｆｏｒｅｇｒｏｕｎｄ）の処理であり、学習処理部１１０の処理は裏（ｂａｃｋｇｒｏｕｎｄ）の処理であると言える。
【０１８３】
システム全体で見ると、分離処理部１２３は、ＡＤ変換部１２２から入力する音声信号に対して分離行列を随時適用することで分離結果を生成する一方で、学習処理部１３０の処理の進行に応じて、分離行列を最新のものに置き換えるという動作をする。
【０１８４】
音源方向推定部１２５では、分離処理部１２３から、観測信号と分離結果（または、式［６．５］の関係）を入力して、観測信号と分離行列とから音源方向を算出する。
以下では、それぞれの構成要素の処理について説明する。
【０１８５】
複数のマイクロホン１２１で収録された音はＡＤ変換部１２２でデジタルデータに変換され、分離処理部１２３のフーリエ変換部１２６へ送られる。フーリエ変換部１２６では、窓つきの短時間フーリエ変換によって周波数領域のデータへ変換される（詳細は後述する）。その際、フレームと呼ばれる一定個数のデータが生成される。以降の処理は、このフレームを単位として行なわれる。フーリエ変換されたデータは、音源方向推定部１２５、分離部１２７、および学習処理部１１０のスレッド制御部１１２にそれぞれ送られる。
【０１８６】
分離部１２７では、フーリエ変換部１２６から入力する観測信号のフーリエ変換結果の１フレーム分であるＸ（ｔ）と、そのタイミングより前に求まっている分離行列Ｗを学習処理部１１０の分離行列保持部１１３から入力して、これらの入力を乗じて、分離結果であるＹ（ｔ）を生成する。
Ｙ（ｔ）＝ＷＸ（ｔ）
である。ｔはフレーム番号である。
【０１８７】
分離結果Ｙ（ｔ）は、フーリエ逆変換部１２８および音源方向推定部１２５へ送られる。音源方向推定部１２５へは、分離結果の代わりに分離行列を送ってもよい。フーリエ逆変換部１２８へ送られた分離結果［Ｙ（ｔ）］は時間領域の信号に変換され、後段の処理を実行する後段処理部ａ１２４ａへ送られる。後段の処理とは、例えば音声認識や話者識別などである。後段の処理によっては、周波数領域のデータをそのまま使用することも可能であり、その場合、フーリエ逆変換は省略可能である。
【０１８８】
一方、学習処理部１１０のスレッド制御部１１２へ送られた観測信号（フーリエ変換部１２６から入力する観測信号のフーリエ変換結果）は、複数の学習スレッド１１１−１〜Ｎを利用した演算を行うスレッド演算部１１１へ送られる。スレッド演算部１１１は、それぞれの学習スレッドに与えられた観測信号を一定量だけ蓄積した後、ＩＣＡのバッチ処理を用いて観測信号から分離行列を求める。求まった分離行列は、スレッド制御部１１２を経由して分離処理部１２３の分離部１２７へ送られる。
【０１８９】
さらに、音源方向推定部１２５へ送られた観測信号（フーリエ変換部１２６から入力する観測信号のフーリエ変換結果）は、分離部１２７で生成される分離結果と共に（または、分離行列と共に）、相互共分散行列の計算に使用され、その後、音源方向が計算される。（詳細は後述する。）
【０１９０】
こうして求めた音源方向は、後段の処理を実行する後段処理部ｂ１２４ｂへ送られる。後段処理部ｂ１２４ｂの実行する後段の処理とは、例えばダイアリゼーション（ｄｉａｒｉｚａｔｉｏｎ）と呼ばれる処理である。すなわち、音声認識や話者識別の結果などと組み合わせることで、「どの時刻に、どの方向から（どの位置で）、誰が、何を発話したか」という情報のセットを生成することができる。
【０１９１】
次に、音源方向推定部１２５の詳細について、図１３および図１４を参照して説明する。図１３は相互共分散行列を観測信号と分離結果とから計算する（式［６．１］）場合の構成図、図１４は観測信号と分離行列とから計算する（式［６．５］）場合の構成図である。
【０１９２】
まず、図１３を参照して相互共分散行列を観測信号と分離結果とから計算して音源方向の推定処理を実行する音源方向推定部１２５の一構成例について説明する。
【０１９３】
図１３に示す音源方向推定部１２５の入力は、
図１２に示す構成におけるフーリエ変換部１２６から供給される１フレーム分の観測信号Ｘ（ｔ）、
図１２に示す構成における分離部１２７から供給される１フレーム分の分離結果Ｙ（ｔ）、
これらのデータである。
【０１９４】
フーリエ変換部１２６から供給される１フレーム分の観測信号Ｘ（ｔ）は、観測信号バッファ１３３に格納される。これは、図９に示す現在時刻を含む区間の観測信号側の観測信号データ区間９２に相当する。観測信号バッファ１３３は長さＬのリングバッファであり、Ｌ＋１番目の以降の入力はバッファ内で最も古い信号に上書きされる。例えば、Ｘ（０）からＸ（Ｌ−１）まではそのまま観測信号バッファ１３３に格納されるが、Ｘ（Ｌ）を格納する際はＸ（０）の領域に上書きされる。もちろん、必要に応じて、過去の観測信号をすべて保存しても構わない。
【０１９５】
分離部１２７から供給される１フレーム分の分離結果Ｙ（ｔ）は、分離結果バッファ１３４に格納される。これは、図９に示す現在時刻を含む区間の分離結果側の分離結果データ区間９４に相当する。分離結果バッファ１３４は観測信号バッファ１３３と同様に、長さＬのリングバッファである。
【０１９６】
演算部１３２は、Ｌフレーム分の観測信号と分離結果とから周波数ビンごとに相互共分散行列を計算し、結果を相互共分散行列バッファ１３５に格納する。先に図９を参照して説明したように、時間周波数領域において対応する時間区間の観測信号スペクトログラムと分離信号スペクトログラムとの間で相互共分散行列を計算し、結果を相互共分散行列バッファ１３５に格納する。さらに、演算部１３２は、相互共分散行列から、チャンネルごとかつマイクペアごとの要素比や位相差（式［８．１］のｒｉｍｋ（ω）や式［８．９］のαｉｍｋ（ω）など）を計算し、計算結果を位相差バッファ１３６に格納する。
【０１９７】
次に、演算部１３２は、位相差バッファ１３６に格納された位相差に対して最も当てはまる比例定数（直線の傾き）Ｋｉｍｋを求め、結果を直線傾きバッファ１３７に格納し、また、当てはめ誤差ｅｒｒｉｍｋを当てはめ誤差バッファ１３８に格納する。
【０１９８】
最後に、演算部１３２は、直線傾きバッファ１３７に格納した比例定数（直線の傾き）Ｋｉｍｋと当てはめ誤差ｅｒｒｉｍｋから、音源方向［θｋ］を計算し、計算結果としての音源方向［θｋ］を音源方向バッファ１３９に格納する。
【０１９９】
なお、必要に応じて、音源方向の角度の代わりに、方向を表わす３次元（用途によっては２次元）のベクトルや、音源の位置を表わす３次元（用途によっては２次元）の座標を計算することも可能である。
【０２００】
音源方向推定部１２５における各演算は演算部１３２において行なわれる。この図では省略されているが、演算部１３２はこの図の各バッファとデータのやり取りが行なえるように結線されている。
【０２０１】
次に、図１４を参照して、観測信号と分離行列とから音源方向の推定処理を実行する音源方向推定部１２５の一構成例について説明する。
【０２０２】
図１４に示す音源方向推定部１２５の入力は、
図１２に示す構成におけるフーリエ変換部１２６から供給される１フレーム分の観測信号Ｘ（ｔ）、
同時刻において図１２に示す分離部１２７で使用されている分離行列Ｗ、
これらのデータである。
【０２０３】
フーリエ変換部１２６から供給される１フレーム分の観測信号Ｘ（ｔ）は、観測信号バッファ１４３に格納される。これは、図９に示す現在時刻を含む区間の観測信号側の観測信号データ９２に相当する。観測信号バッファ１４３は長さＬのリングバッファである。
【０２０４】
本構成では、演算部１４２は、観測信号バッファ１４３に格納されたデータから、相互共分散行列ではなく観測信号自体の共分散行列を計算し、結果を共分散行列バッファ１４４に格納する。
【０２０５】
演算部１４２は、同時刻において図１２に示す分離部１２７で使用されている分離行列Ｗを入力し、この分離行列のエルミート転置と、共分散行列バッファ１４４に格納された共分散行列とから、前述した式［６．５］に基づいて相互共分散行列［ΣＸＹ］を計算し、相互共分散行列バッファ１４５に格納する。
【０２０６】
この後の処理は、図１３を参照して説明した処理と同様の処理となる。すなわち、演算部１４２は、相互共分散行列から、チャンネルごとかつマイクペアごとの要素比や位相差（式［８．１］のｒｉｍｋ（ω）や式［８．９］のαｉｍｋ（ω）など）を計算し、計算結果を位相差バッファ１４６に格納する。次に、演算部１４２は、位相差バッファ１４６に格納された位相差に対して最も当てはまる比例定数（直線の傾き）Ｋｉｍｋを求め、結果を直線傾きバッファ１４７に格納し、また、当てはめ誤差ｅｒｒｉｍｋを当てはめ誤差バッファ１４８に格納する。
【０２０７】
最後に、演算部１４２は、直線傾きバッファ１４７に格納した比例定数（直線の傾き）Ｋｉｍｋと当てはめ誤差ｅｒｒｉｍｋから、音源方向［θｋ］を計算し、計算結果としての音源方向［θｋ］を音源方向バッファ１４９に格納する。
【０２０８】
図１４に示す構成では、分離結果についてはＬフレーム分の信号を格納する必要が
ないため、図１３の構成と比べてバッファのメモリを節約できるという利点がある。
【０２０９】
次に、図１２に示す装置構成における学習処理部１１０のスレッド制御部１１２の詳細構成について、図１５を参照して説明する。
【０２１０】
現フレーム番号保持カウンタ１５１は、観測信号が１フレーム分供給されるごとに値が１インクリメントされ、所定の値に達すると初期値に戻るようになされている。なお、カウンタのインクリメントは、フレーム番号と同期させる代わりに、時間領域信号のサンプル番号と同期させても構わない。
【０２１１】
学習初期値保持部１５２は、それぞれのスレッドにおいて学習処理を実行する場合の分離行列Ｗの初期値を保持するものである。分離行列Ｗの初期値は、基本的には最新の分離行列と同一だが、異なる値を用いるものとしても良い。例えば、学習初期値にはリスケーリング（周波数ビン間のパワーを調整する処理。詳細は後述）を適用する前の分離行列を用い、分離行列にはスケーリング適用後のものを用いるようにする。
【０２１２】
蓄積開始予定タイミング指定情報保持部１５３は、蓄積を開始するタイミングを複数のスレッド間で一定間隔にするために用いられる情報である。なお、蓄積開始予定タイミングは、相対時刻を用いて表されていてもよいし、相対時刻の代わりにフレーム番号で管理してもよいし、時間領域信号のサンプル番号で管理しても良い。これについては他の「タイミング」を管理するための情報についても同様である。
【０２１３】
観測信号蓄積タイミング情報保持部１５４は、分離部１２７で現在使用されている分離行列Ｗが、どのタイミングで取得された観測信号を基に学習されたものであるかを示す情報、すなわち、最新の分離行列に対応した観測信号の相対時刻またはフレーム番号を保持するものである。観測信号の蓄積タイミング情報保持部１５４には、対応する観測信号の蓄積開始タイミングと蓄積終了タイミングとの両方を格納しても良いが、ブロック長、すなわち、観測信号の蓄積時間が一定ならば、いずれか一方だけを保存すれば十分である。
【０２１４】
また、スレッド制御部１１２は、それぞれのスレッドへリンクされているポインタを保持したポインタ保持部１５５を有し、これを用いて、複数のスレッド１１１−１〜Ｎの処理を制御している。
【０２１５】
次に、図１６を参照してスレッド演算部１１１において実行する処理について説明する。
【０２１６】
スレッド１１１−１〜Ｎのそれぞれは、観測信号バッファ１６１、分離結果バッファ１６２、学習演算部１６３、および、分離行列保持部１６４の各モジュールの機能を用いて、バッチ処理のＩＣＡを実行する。
【０２１７】
観測信号バッファ１６１は、スレッド制御部１１２から供給される観測信号を保持するものであり、その容量は、１ブロック長に対応する観測信号の容量と同じか、それよりも大きいものである。ただし、後述する「観測信号のフレーム間引き」を行なう場合は、間引く分だけバッファの大きさを小さくしても構わない。
【０２１８】
分離結果バッファ１６２には、学習演算部１６３により演算された、分離行列収束前の分離結果が保持される。
【０２１９】
学習演算部１６３は、観測信号バッファ１６１に蓄積されている観測信号を、分離行列保持部１６４に保持されている分離処理用の分離行列Ｗに基づいて分離して、分離結果バッファ１６２に蓄積するとともに、分離結果バッファ１６２に蓄積される分離結果を用いて、学習中の分離行列を更新する処理を実行する。
【０２２０】
また、スレッドは、状態遷移マシンであり、現在の状態は、状態格納部１６５に格納されている。そして、スレッドの状態は、カウンタ１６６のカウンタ値によって、スレッド制御部１１２により制御される。カウンタ１６６は、観測信号が１フレーム分供給されるのと同期して値が変化し、この値によって状態を切り替える。詳細は後述する。
【０２２１】
観測信号の開始・終了タイミング保持部１６７には、学習に使用されている観測信号の開始タイミングと終了タイミングを示す情報のうちの少なくともいずれか一方が保持されている。タイミングを示す情報は、上述したように、フレーム番号やサンプル番号であっても良いし、相対時刻情報であっても良い。ここでも、開始タイミングと終了タイミングとの両方を格納しても良いが、ブロック長、すなわち、観測信号の蓄積時間が一定ならば、いずれか一方だけを保存すれば十分である。
【０２２２】
学習終了フラグ１６８は、学習が終了したことをスレッド制御部１１２に通知するために用いられるフラグである。スレッドの起動時においては、学習終了フラグ１６８はＯＦＦ（フラグが立っていない）にセットされ、学習が終了した時点でＯＮにセットされる。そして、スレッド制御部１１２が、学習が終了したことを認識した後、スレッド制御部１１２の制御により、学習終了フラグ１６８は、再び、ＯＦＦにセットされる。
【０２２３】
前処理用データ保持部１６９は、前処理が施された観測信号を元に戻す際に必要となるデータを保存しておく領域である。具体的には、例えば、前処理において、観測信号の正規化（分散を１に、平均を０にそろえる）が実行される場合、前処理用データ保持部１６９には、分散（または標準偏差やその逆数）や平均などの値が保持されるので、これを用いて正規化前の信号を復元することができる。また、例えば、前処理として無相関化（ｐｒｅ−ｗｈｉｔｅｎｉｎｇとも称される）が実行される場合、前処理用データ保持部１６９には、無相関化で乗じた行列が保持される。
【０２２４】
状態格納部１６５、カウンタ１６６、観測信号の開始・終了タイミング保持部１６７に保持される値は、スレッド制御部１１２等の外部モジュールから値を書き換えることができる。例えば、このスレッドにおいて学習処理ループが回っている間であっても、スレッド制御部１１２は、カウンタ１６６の値を変更することができる。
【０２２５】
状態遷移の実装については、それぞれのスレッド１１１−１〜Ｎが自分自身のカウンタ１６６の値に基づいて状態を自発的に変化させるという仕様にしても良いし、スレッド制御部１１２がカウンタ１６６の値や学習終了フラグ１６８の値に応じて、対応するスレッドに対して「指定された状態に遷移せよ」というコマンド（以降「状態遷移コマンド」と称するものとする）を発行し、それぞれのスレッドは、供給されたコマンドに応じて、その状態を遷移させるという仕様としてもよい。
【０２２６】
次に、学習スレッド１１１−１〜Ｎの状態遷移について図１７および図１８を参照して説明する。実装としては、カウンタ１６７の値に基づいて学習スレッド自身が自分の状態を変化させるという仕様にしても良いが、スレッド制御部がカウンタの値や「学習終了フラグ」１６８の値に応じて状態遷移コマンドを発行し、学習スレッドはそのコマンドを受けて状態を変化させるという仕様でも良い。以下の実施例では、後者の仕様を採用している。
【０２２７】
図１７は、先に図７を参照して説明したスレッド中の１つを示している。スレッドの各々は、観測信号の「蓄積中」状態において、指定された時間、すなわち１ブロック長の観測信号がバッファに蓄えられる。指定された時間が経過した後、状態は、学習中に遷移する。
【０２２８】
学習中の状態において、分離行列Ｗが収束するまで（または一定回数）学習処理ループが実行され、蓄積中の状態において蓄積された観測信号に対応した分離行列が求められる。分離行列Ｗが収束した後（または一定回数の学習処理ループが実行された後）状態は、待機中に遷移する。
【０２２９】
そして、待機中の状態において、指定された時間だけ、観測信号の蓄積や学習は実行されず、待機される。待機中の状態を維持する時間は、学習にかかった時間によって決まる。すなわち、図１７に示されるように、予め、「蓄積中」状態と「学習中」状態と「待機中」状態との合計の時間幅であるスレッド長（ｔｈｒｅａｄ＿ｌｅｎ）が定められ、基本的には、「学習中」状態が終了したときからスレッド長が終了するまでの間の時間が、「待機中」状態の時間（待機時間）とされる。待機時間が過ぎた後、状態は、観測信号の「蓄積中」状態へ戻る。
【０２３０】
これらの時間は、例えば、ミリ秒などの単位で管理してもよいが、短時間フーリエ変換で生成されるフレームを単位として計測するものとしても良い。以降の説明では、フレームを単位として計測する（たとえば、カウントアップをおこなう）ものとする。
【０２３１】
図１８を参照してスレッドの状態遷移についてさらに説明する。システムの起動直後は各スレッドは「初期状態」１８１にあるが、その内の１つは「蓄積中」１８３へ、残りは「待機中」１８２へ状態を遷移させる（状態遷移コマンドを発行する）。先に説明した図７を例にとると、スレッド１が「蓄積中」へ遷移したスレッド、それ以外が「待機中」へ遷移したスレッドである。以下では、「蓄積中」へ遷移したスレッドについて先に説明する。
【０２３２】
観測信号の蓄積に要する時間をブロック長（ｂｌｏｃｋ＿ｌｅｎ）と呼ぶ（図１７参照）。また、蓄積・学習・待機の１周期に要する時間をスレッド長（ｔｈｒｅａｄ＿ｌｅｎ）と呼ぶ。これらの時間は、ミリ秒などの単位で管理してもよいが、短時間フーリエ変換で生成されるフレームを単位としても良い。以降の説明では、フレームを単位としている。
【０２３３】
「蓄積中→学習中」および「待機中→蓄積中」の状態遷移については、カウンタの値に基づいて行なう。すなわち、「蓄積中」（図１７の蓄積中状態１７１および図１８の蓄積中状態１８３）から始まったスレッドの中では、観測信号が１フレーム分供給されるごとにカウンタを＋１し、カウンタの値がブロック長（ｂｌｏｃｋ＿ｌｅｎ）と同じ値になったら、「学習中」（図１７の学習中状態１７２および図１８の学習中状態１８４）へ状態を遷移させる。学習は、分離処理と並列にバックグラウンドで行なわれるが、その間も観測信号のフレームと連動してカウンタを＋１していく。
【０２３４】
学習が終了したら、状態を「待機中」（図１７の待機中状態１７３および図１８の待機中状態１８２）へ遷移させる。待機中も学習中と同様に、観測信号のフレームと連動してカウンタを＋１していく。そしてカウンタがスレッド長（ｔｈｒｅａｄ＿ｌｅｎ）と同じ値になったら、状態を「蓄積中」（図１７の蓄積中状態１７１および図１８の蓄積中状態１８３）へ遷移させると共に、カウンタを０（または適切な初期値）に戻す。
【０２３５】
一方、「初期状態」１８１から「待機中」（図１７の待機中状態１７３および図１８の待機中状態１８２）へ遷移したスレッドについては、待機させたい時間に対応した値にカウンタをセットする。例えば、図７のスレッド２は、ブロックのシフト幅（ｂｌｏｃｋ＿ｓｈｉｆｔ）だけ待機してから「蓄積中」へ遷移する。同様に、スレッド３はブロックのシフト幅の２倍（ｂｌｏｃｋ＿ｓｈｉｆｔ×２）だけ待機している。これらを実現するためには、
スレッド２のカウンタを、
（スレッド長）−（ブロックシフト幅）：（ｔｈｒｅａｄ＿ｌｅｎ）−（ｂｌｏｃｋ＿ｓｈｉｆｔ）にセットする。
また、スレッド３のカウンタを、
（スレッド長）−（２×ブロックシフト幅）：（ｔｈｒｅａｄ＿ｌｅｎ）−（ｂｌｏｃｋ＿ｓｈｉｆｔ×２）にセットする。
このような設定にすれば、カウンタの値がスレッド長（ｔｈｒｅａｄ＿ｌｅｎ）に達した後で「蓄積中」へ遷移し、それ以降はスレッド１と同様に「蓄積中→学習中→待機中」の周期を繰り返す。
【０２３６】
学習スレッドを何個用意する必要があるかは、スレッド長とブロックのシフト幅とで決める。スレッド長を  ｔｈｒｅａｄ＿ｌｅｎ、ブロックのシフト幅をｂｌｏｃｋ＿ｓｈｉｆｔとすると、必要な個数は、
（スレット長）／（ブロックシフト幅）、すなわち、
ｔｈｒｅａｄ＿ｌｅｎ／ｂｌｏｃｋ＿ｓｈｉｆｔ
で求まる。
なお、端数は切り上げる。
【０２３７】
例えば図７では、
［スレッド長（ｔｈｒｅａｄ＿ｌｅｎ）］＝１．５×［ブロック長（ｂｌｏｃｋ＿ｌｅｎ）］、
［ブロックシフト幅（ｂｌｏｃｋ＿ｓｈｉｆｔ）］＝０．２５×ブロック長（ｂｌｏｃｋ＿ｌｅｎ）］
に設定してあるため、必要なスレッド数は１．５／０．２５＝６である。
【０２３８】
［リアルタイム分離・音源方向推定処理の処理例］
以下では、本発明の装置によるリアルタイム分離・音源方向推定処理について説明する。
図１２の装置構成において、音源方向推定部１２５の処理と分離処理部１２３の「表の処理（分離）」とについては、どちらもフレームと同期した処理が可能である。そこで、両者の処理をまとめて図１９に示すフローチャートを参照して説明する。
【０２３９】
なお、学習処理部１１０の「裏の処理（学習）」については、音源方向推定＋分離処理とは別の処理単位（別スレッド・別プロセス・別プロセッサなど）で動かすため、別のフローチャートを用いて説明する。また、両者でやりとりするコマンド等については、図２２に示すシーケンス図で説明する。
【０２４０】
最初に、図１９に示すフローチャートを参照して音源方向推定部１２５の処理と分離処理部１２３の処理について説明する。
システムが起動したら、ステップＳ１０１において、各種の初期化を行なう。初期化の詳細は後述する。ステップＳ１０３の音入力から、ステップＳ１１０の送信までの処理を、システムでの処理が終了（ステップＳ１０２でＹｅｓ）するまで繰り返す。
【０２４１】
ステップＳ１０３の音入力は、オーディオデバイス（実施形態によってはネットワークやファイルなど）から一定数のサンプルを取り込み（この処理を「キャプチャー」と呼ぶ）、バッファに蓄える処理である。これをマイクの個数分だけ行なう。以降では、キャプチャーされたデータを観測信号と呼ぶ。
【０２４２】
次に、ステップＳ１０４において、観測信号を一定長ごとに切り出して、短時間フーリエ変換（ＳＴＦＴ）を行う。短時間フーリエ変換の詳細について図２０を参照して説明する。
【０２４３】
例えば図２に示すような環境においてｋ番目のマイクによって収録された観測信号ｘｋを図２０（ａ）に示す。この観測信号ｘｋから一定長を切り出した切り出しデータであるフレーム１９１〜１９３にハニング窓やサイン窓等の窓関数を作用させる。なお、切り出した単位をフレームと呼ぶ。１フレーム分のデータに対して、離散フーリエ変換（有限区間のフーリエ変換のこと。略称DFT）または高速フーリエ変換（FFT）を施すことにより、周波数領域のデータであるスペクトルＸｋ（ｔ）を得る（ｔはフレーム番号）。
【０２４４】
切り出すフレームの間には、図に示すフレーム１９１〜１９３のように重複があってもよく、そうすることで連続するフレームのスペクトルＸｋ（ｔ−１）〜Ｘｋ（ｔ＋１）を滑らかに変化させることができる。また、スペクトルをフレーム番号に従って並べたものをスペクトログラムと呼ぶ。図２０（ｂ）がスペクトログラムの例である。
【０２４５】
なお、短時間フーリエ変換（ＳＴＦＴ）において切り出すフレーム間に重複がある場合は、逆フーリエ変換（ＦＴ）においてもフレームごとの逆変換結果（波形）波形を重複つきで重ね合わせる。これをオーバーラップ加算（ｏｖｅｒｌａｐ  ａｄｄ）という。逆変換結果は、オーバーラップ加算（ｏｖｅｒｌａｐ  ａｄｄ）の前にサイン窓等の窓関数を再び作用させても良く、これを、ｗｅｉｇｈｔｅｄ  ｏｖｅｒｌａｐ  ａｄｄ（ＷＯＬＡ）という。ＷＯＬＡにより、フレーム間の不連続性に由来するノイズを低減することができる。
【０２４６】
本発明では入力チャンネルが複数（マイクの個数分）あるため、フーリエ変換もチャンネル数だけ行なう。以降では、全チャンネル・１フレーム分のフーリエ変換結果をＸ（ｔ）というベクトルで表わす（先に説明した式［４．１１］），式［４．１１］において、ｎはチャンネル数（＝マイク数）である。Ｍは周波数ビンの総数であり、短時間フーリエ変換のポイント数をＴとすると、Ｍ＝Ｔ／２＋１である。
【０２４７】
図１９のフローに戻り、説明を続ける。ステップＳ１０４において、観測信号を一定長ごとに切り出して、短時間フーリエ変換（ＳＴＦＴ）を行った後、ステップＳ１０５において、それぞれの学習スレッドに対する制御を行なう。詳細は後述する。
【０２４８】
次に、ステップＳ１０５で生成された観測信号Ｘ（ｔ）に対して、ステップＳ１０６において分離を行なう。分離行列をＷとすると（式［４．１０］）、分離結果Ｙ（ｔ）（式［４．４］）は
Ｙ（ｔ）＝ＷＸ（ｔ）
で求まる（式［４．１２］）。
【０２４９】
次に、ステップＳ１０７において、分離結果Ｙ（ｔ）に対して逆フーリエ変換（逆ＦＴ）をかけることで、時間領域の信号に戻す。それからステップＳ１０８において、分離結果を後段の処理へ送信する。
【０２５０】
次に、ステップＳ１０９において、１フレーム分の観測信号と分離結果とを用いて（または、１フレーム分の観測信号と最新の分離行列をと用いて）、音源の方向を更新する。詳細は後述する。ステップＳ１１０では、そうして求めた音源方向を、後段の処理へ送信する。以上のステップＳ１０３〜Ｓ１１０を、終了まで繰り返す。
【０２５１】
図１９に示すフローチャートにおけるステップＳ１０１の初期化処理の詳細について、図２１に示すフローチャートを参照して説明する。
【０２５２】
ステップＳ１５１において、図１２および図１５に示すスレッド制御部１１２は、自分自身を初期化する。具体的には、スレッド制御部１１２において、現フレーム番号保持カウンタ１５１（図１５参照）が初期化されてその値が０とされ、学習初期値保持部１５２（図１５参照）に適切な初期値が代入される。例えば、初期値は単位行列でも良いし、前回のシステム終了時の分離行列Ｗが保存されている場合は、前回のシステム終了時の分離行列Ｗ、またはこの分離行列に適切な変換を作用させたものを使用しても良い。また、例えば、画像や先見知識等の情報により、音源の方向がある程度の精度で推定できるような場合には、音源方向に基づいて初期値を算出して、設定するものとしてもよい。
【０２５３】
そして、蓄積開始予定タイミング指定情報保持部１５３には、（必要スレッド数−１）×［ブロックシフト幅（ｂｌｏｃｋ＿ｓｈｉｆｔ）］の値が設定される。この値は、一番大きなスレッド番号を有するスレッドの蓄積が開始するタイミング（フレーム番号）である。そして、観測信号の蓄積タイミング情報保持部１３４には、最新の分離行列に対応した観測信号を示すタイミング情報（ブレーム番号または相対時刻情報）が保持されるので、ここでは、初期化されて、０が保持される。
【０２５４】
なお、分離行列保持部１１３（図１２参照）にも、初期化された場合の学習初期値保持部１５２と同様に、適切な初期値が保持される。すなわち、分離行列保持部１１３に保持される初期値は、単位行列でも良いし、前回のシステム終了時の分離行列が保存されている場合は、前回のシステム終了時の分離行列Ｗ、またはこの分離行列に適切な変換を作用させたものを使用しても良い。また、例えば、画像や先見知識等の情報により、音源の方向がある程度の精度で推定できるような場合には、音源方向に基づいて初期値を算出して、設定するものとしてもよい。
【０２５５】
ステップＳ１５２において、スレッド制御部１１２は、スレッド演算部１１１において実行されるスレッドを必要な数ｉだけ確保し、それらの状態を「初期化」状態とする。
【０２５６】
ここで、必要なスレッドの数ｉは、スレッド長／ブロックシフト幅（ｔｈｒｅａｄ＿ｌｅｎ／ｂｌｏｃｋ＿ｓｈｉｆｔ）の小数点以下を切り上げる（すなわち、ｔｈｒｅａｄ＿ｌｅｎ／ｂｌｏｃｋ＿ｓｈｉｆｔよりも大きく最も値の近い整数）ことにより求められる。
【０２５７】
ステップＳ１５３において、スレッド制御部１１２は、スレッドループを開始して、全てのスレッドの初期化が終了するまで、初期化未処理のスレッドを検出して、ステップＳ１５４乃至ステップＳ１５９の処理を実行する。
【０２５８】
ステップＳ１５４において、スレッド制御部１１２は、スレッド番号は１であるか否かを判断する。
【０２５９】
ステップＳ１５４において、スレッド番号は１であると判断された場合、ステップＳ１５５において、スレッド制御部１１２は、スレッド番号１のスレッド（例えば、スレッド１１１−１）を制御して、そのカウンタ１６６（図１６参照）を初期化（例えば、０にセット）する。
【０２６０】
ステップＳ１５６において、スレッド制御部１１２は、スレッド番号１のスレッド（例えば、スレッド１１１−１）に、「蓄積中」状態に状態を遷移させるための状態遷移コマンドを発行して、処理は、後述するステップＳ１５９にすすむ。
【０２６１】
ステップＳ１５４において、スレッド番号は１ではないと判断された場合、ステップＳ１５７において、スレッド制御部１１２は、対応するスレッド（スレッド１１１−２乃至スレッド１１１−ｉのうちのいずれか）のカウンタ１６６の値を、ｔｈｒｅａｄ＿ｌｅｎ−ｂｌｏｃｋ＿ｓｈｉｆｔ×（スレッド番号−１）に設定する。
【０２６２】
ステップＳ１５８において、スレッド制御部１１２は、「待機中」状態に状態を遷移させるための状態遷移コマンドを発行する。
【０２６３】
ステップＳ１５６、または、ステップＳ１５８の処理の終了後、ステップＳ１５９において、スレッド制御部１１２は、スレッド内のまだ初期化されていない情報、すなわち、状態格納部１６５（図１６参照）に格納された状態を示す情報、および、カウンタ１６６のカウンタ値以外の情報を初期化する。具体的には、例えば、スレッド制御部１１２は、学習終了フラグ１６８（図１６参照）をＯＦＦにセットし、観測信号の開始・終了タイミング保持部１６７、および、前処理用データ保持部１６９の値を初期化（例えば、０にセット）する。
【０２６４】
スレッド演算処理部１１２に確保された全てのスレッド、すなわち、スレッド１１１−１乃至スレッド１１１−ｉが初期化された場合、ステップＳ１６０において、スレッドループが終了され、初期化は終了する。
【０２６５】
このような処理により、スレッド制御部１２１は、スレッド演算部１１１に確保された複数のスレッドのすべてを初期化する。
【０２６６】
なお、ここでは、立ち上げられたスレッドの数ｉだけ、ステップＳ１５３乃至ステップＳ１６０の処理のループが繰り返されるものとして説明したが、処理ループを繰り返す代わりに、スレッドの個数ｉの並列処理を実行するものとしてもよい。これ以降の処理ループを繰り返す部分についても、同様に、処理ループを繰り返す代わりに、並列処理を実行するようにしても良い。
【０２６７】
なお、図２１のステップＳ１５４〜Ｓ１５８の処理は、図２２のシーケンス図では最初の「初期化」処理と、その直後の状態遷移コマンド送信に該当する。なお、図２２は、スレッド制御部１１２による複数の学習スレッド１，２に対する制御シーケンスを示している。各スレッドは、待機、蓄積、学習という処理を繰り返し実行する。スレッド制御部は観測信号を各スレッドに提供して、各スレッドが観測データを蓄積した後、学習処理を行い、分離行列を生成してスレッド制御部に提供する。
【０２６８】
次に、図２３のフローチャートを参照して、図１９に示すフローチャートにおけるステップＳ１０５においてスレッド制御部１１２によって実行される、スレッド制御処理について説明する。
【０２６９】
なお、このフローチャートは、スレッド制御部１１２から見たものであり、学習スレッド１１１−１−Ｎから見たものではないことに注意されたい。例えば「学習中処理」とは、学習スレッドの状態が「学習中」であるときにスレッド制御部１１２が行なう処理のことである。（学習スレッド自体の処理については、図３０を参照されたい。）
【０２７０】
ステップＳ２０１〜Ｓ２０６は学習スレッドについてのループであり、図２１に示すフローのステップＳ１５２で生成された個数だけループを回す（並列処理でも構わない）。ステップＳ２０２において、学習スレッドの現在の状態を状態保持部１６５（図１６参照）から読み込み、その値によってそれぞれ「待機中処理」・「蓄積中処理」・「学習中処理」を実行する。それぞれの処理の詳細は後述する。
【０２７１】
フローの各ステップについて説明する。ステップＳ２０１において、スレッド制御部１１２は、スレッドループを開始し、制御実行するスレッドのスレッド番号を示す変数ｓをｓ＝１として、１つのスレッドの処理が終了すると変数ｓを１インクリメントして、ｓ＝ｉとなるまで、ステップＳ２０２乃至ステップＳ２０７のスレッドループの処理を繰り返し実行する。
【０２７２】
ステップＳ２０２において、スレッド制御部１１２は、変数ｓで示されるスレッド番号のスレッドの状態格納部１６５に保持されている、そのスレッドの内部状態を示す情報を取得する。変数ｓで示されるスレッド番号のスレッドの状態として、「待機中」状態であると検出された場合、ステップＳ２０３において、スレッド制御部１１２は、図２４のフローチャートを用いて後述する待機中状態における処理を実行し、処理は、後述するステップＳ２０６に進む。
【０２７３】
ステップＳ２０２において、変数ｓで示されるスレッド番号のスレッドの状態が「蓄積中」状態であると検出された場合、ステップＳ２０４において、スレッド制御部１１２は、図２５のフローチャートを用いて後述する蓄積中状態における処理を実行し、処理は、後述するステップＳ２０６に進む。
【０２７４】
ステップＳ２０２において、変数ｓで示されるスレッド番号のスレッドの状態が「学習中」状態であると検出された場合、ステップＳ２０５において、スレッド制御部１１２は、図２６のフローチャートを用いて後述する学習中状態における処理を実行する。
【０２７５】
ステップＳ２０３、ステップＳ２０４、または、ステップＳ２０５の処理の終了後、ステップＳ２０６において、スレッド制御部１１２は、変数ｓを１インクリメントする。そして、制御実行するスレッドのスレッド番号を示す変数ｓが、ｓ＝ｉとなったとき、スレッドループを終了する。
【０２７６】
ステップＳ２０７において、スレッド制御部１１２は、現フレーム番号保持カウンタ１５１（図１５参照）に保持されているフレーム番号を１インクリメントし、スレッド制御処理を終了する。
【０２７７】
このような処理により、スレッド制御部１１２は、複数のスレッドの全てを、それらの状態に応じて制御することができる。
【０２７８】
なお、ここでは、立ち上げられたスレッドの数ｉだけ、スレッドループが繰り返されるものとして説明したが、スレッドループを繰り返す代わりに、スレッドの個数ｉの並列処理を実行するものとしてもよい。
【０２７９】
次に、図２４のフローチャートを参照して、図２３に示すフローチャートにおけるステップＳ２０３において実行される、待機中状態における処理について説明する。
【０２８０】
この待機中状態における処理は、図２３を用いて説明したスレッド制御処理における変数ｓに対応するスレッドの状態が「待機中」状態であるときに、スレッド制御部１１２において実行される処理である。
【０２８１】
ステップＳ２１１において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６（図１６参照）を、１インクリメントする。
【０２８２】
ステップＳ２１２において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６の値は、スレッド長（ｔｈｒｅａｄ＿ｌｅｎ）より小さいか否かを判断する。ステップＳ２１２において、カウンタ１６６の値は、スレッド長より小さいと判断された場合、待機中処理を終了し、図２３のステップＳ２０６に進む。
【０２８３】
ステップＳ２１２において、カウンタ１６６の値がスレッド長より小さくないと判断された場合、ステップＳ２１３において、スレッド制御部１１２は、「蓄積中」状態に状態を遷移させるための状態遷移コマンドを、対応するスレッド１１１に発行する。
【０２８４】
すなわち、スレッド制御部１１２は、図１８を用いて説明した状態遷移図において、「待機中」であるスレッドを、「蓄積中」に遷移させるための状態遷移コマンドを発行する。
【０２８５】
ステップＳ２１４において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６（図１６参照）を初期化（例えば、０にセット）し、観測信号の開始・終了タイミング保持部１６７（図１６参照）に、観測信号の蓄積開始タイミング情報、すなわち、スレッド制御部１１２の現フレーム番号保持カウンタ１５１（図１５参照）に保持されている現在のフレーム番号、または、それと同等の相対時刻情報などを設定して、待機中処理を終了し、図２３のステップＳ２０６に進む。
【０２８６】
このような処理により、スレッド制御部１１２は、「待機中」状態であるスレッドを制御し、そのカウンタ１６６の値に基づいて、「蓄積中」に状態を遷移させることができる。
【０２８７】
次に、図２５のフローチャートを参照して、図２３に示すフローチャートのステップＳ２０４において実行される蓄積中状態における処理について説明する。
【０２８８】
この蓄積中状態における処理は、図２３を用いて説明したスレッド制御処理における変数ｓに対応するスレッドの状態が「蓄積中」状態であるときに、スレッド制御部１１２において実行される処理である。
【０２８９】
ステップＳ２２１において、スレッド制御部１１２は、１フレーム分の観測信号Ｘ（ｔ）を、学習のために、対応するスレッド１１１に供給する。この処理は、図２２に示すスレッド制御部からそれぞれのスレッドへの観測信号の供給に対応する。
【０２９０】
ステップＳ２２２において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６を、１インクリメントする。
【０２９１】
ステップＳ２２３において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６の値がブロック長（ｂｌｏｃｋ＿ｌｅｎ）より小さいか否か、換言すれば、対応するスレッドの観測信号バッファ１６１（図１６参照）が満杯であるか否かを判断する。ステップＳ２２３において、カウンタ１６６の値がブロック長より小さい、換言すれば、対応するスレッドの観測信号バッファ１６１が満杯ではないと判断された場合、蓄積中処理を終了し、図２３のステップＳ２０６に進む。
【０２９２】
ステップＳ２２３において、カウンタ１６６の値がブロック長より小さくない、換言すれば、対応するスレッドの観測信号バッファ１６１が満杯であると判断された場合、ステップＳ２２４において、スレッド制御部１１２は、「学習中」状態に状態を遷移させるための状態遷移コマンドを、対応するスレッド１１１に発行して、蓄積中処理を終了し、図２３のステップＳ２０６に進む。
【０２９３】
すなわち、スレッド制御部１１２は、図１８を用いて説明した状態遷移図において、「蓄積中」であるスレッドを、「学習中」に遷移させるための状態遷移コマンドを発行する。
【０２９４】
このような処理により、スレッド制御部１１２は、「蓄積中」状態であるスレッドに観測信号を供給してその蓄積を制御し、そのカウンタ１６６の値に基づいて、「蓄積中」から「学習中」に状態を遷移させることができる。
【０２９５】
次に、図２６のフローチャートを参照して、図２３に示すフローチャートのステップＳ２０５において実行される、学習中状態における処理について説明する。
【０２９６】
この学習中状態における処理は、図２３を用いて説明したスレッド制御処理における変数ｓに対応するスレッドの状態が「学習中」状態であるときに、スレッド制御部１１２において実行される処理である。
【０２９７】
ステップＳ２３１において、スレッド制御部１１２は、対応するスレッド１１１の学習終了フラグ１６８（図１６参照）がＯＮであるか否かを判断する。ステップＳ２３１において、学習フラグがＯＮであると判断された場合、処理は、後述するステップＳ２３７に進む。
【０２９８】
ステップＳ２３１において、学習フラグがＯＮではないと判断された場合、すなわち、対応するスレッドにおいて学習処理が実行中である場合、ステップＳ２３２に進み、時刻の比較処理を行う。「時刻の比較」とは、学習スレッド１１１内に記録されている、観測信号の開始時刻１６７（図１６参照）と、スレッド制御部１１２に保存されている、現在の分離行列に対応した蓄積開始時刻１５４（図１５参照）とを比較する処理である。スレッド１１１内に記録されている観測信号の開始時刻１６７（図１６参照）が、スレッド制御部１１２に保存されている現在の分離行列に対応した蓄積開始時刻１５４よりも前である場合は、以降の処理をスキップする。
【０２９９】
一方、スレッド１１１内に記録されている観測信号の開始時刻１６７（図１６参照）が、スレッド制御部１１２に保存されている現在の分離行列に対応した蓄積開始時刻１５４よりも後または同じである場合は、ステップＳ２３３に進む。ステップＳ２３３において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６を、１インクリメントする。
【０３００】
次にステップＳ２３４において、スレッド制御部１１２は、対応するスレッド１１１のカウンタ１６６の値がスレッド長（ｔｈｒｅａｄ＿ｌｅｎ）より小さいか否かを判断する。ステップＳ２３４において、カウンタ１６６の値がスレッド長より小さいと判断された場合、学習中処理を終了し、図２３のステップＳ２０６に進む。
【０３０１】
ステップＳ２３４において、カウンタ１６６の値がスレッド長より小さくないと判断された場合、ステップＳ２３５において、スレッド制御部１１２は、カウンタ１６６の値から所定の値を減算し、学習中処理を終了し、図２３のステップＳ２０６に進む。
【０３０２】
学習中にカウンタの値がスレッド長に達した場合とは、学習にかかる時間が長くなってしまい、「待機中」状態の時間が存在しなくなった場合である。その場合、学習はまだ継続しており、観測信号バッファ１６１は利用されているため、次の蓄積を開始することができない。そこで、スレッド制御部１１２は、学習が終了するまで、次の蓄積の開始、すなわち、「蓄積中」状態へ状態を遷移させるための状態遷移コマンドの発行を延期する。そのため、スレッド制御部１１２は、カウンタ１６６の値から所定の値を減算する。減算する値は、例えば、１であっても良いが、それよりも大きな値でも良く、例えば、スレッド長の１０％などといった値であっても良い。
【０３０３】
なお、「蓄積中」状態への状態の遷移の延期を行なうと、蓄積開始時刻がスレッド間で不等間隔となり、最悪の場合、複数のスレッドでほぼ同一の区間の観測信号を蓄積してしまう可能性もある。そうなると、いくつかのスレッドが無意味になるだけでなく、例えば、ＣＰＵ２１が実行するＯＳのマルチスレッドの実装によっては、１つのＣＰＵ２１で複数の学習が同時に動くことになって、更に学習時間が増大し、間隔が一層不均等になってしまう可能性がある。
【０３０４】
そのような事態を防ぐためには、他のスレッドの待機時間を調整して蓄積開始タイミングが再び等間隔になるように調整すればよい。この処理は、ステップＳ２４１において実行される。この待機時間の調整処理の詳細については後述する。
【０３０５】
ステップＳ２３１において、学習終了フラグがＯＮであると判断された場合の処理について説明する。これは、学習スレッド内の学習ループが終了する度に一回実行される処理である。ステップＳ２３１において、学習終了フラグがＯＮであり、対応するスレッドにおいて学習処理が終了したと判断した場合、ステップＳ２３７において、スレッド制御部１１２は、対応するスレッド１１１の学習終了フラグ１６８をＯＦＦにする。この処理は、この分岐が連続実行されるのを防ぐための操作である。
【０３０６】
その後、スレッド制御部１１２はスレッドの打ち切りフラグ１７０（図１６参照）がＯＮであるかＯＦＦであるかを確認し、ＯＮである場合は、ステップＳ２３９において分離行列の更新処理を行い、ステップＳ２４１において待機時間の設定処理を行なう。一方、スレッドの打ち切りフラグ１７０（図１６参照）がＯＦＦである場合は、ステップＳ２３９の分離行列の更新処理は省略し、ステップＳ２４１において待機時間の設定処理を行なう。ステップＳ２３９の分離行列の更新処理と、ステップＳ２４の待機時間の設定処理の詳細については後述する。
【０３０７】
このような処理により、スレッド制御部１１２は、対応するスレッドの学習終了フラグ１６８を参照して、「学習中」状態のスレッドの学習が終了したか否かを判断し、学習が終了した場合、分離行列Ｗを更新し、待機時間を設定するとともに、「学習中」状態から、「待機中」または「蓄積中」に状態を遷移させることができる。
【０３０８】
次に、図２７のフローチャートを参照して、図２６に示すフローチャートのステップＳ２３９において実行される、分離行列更新処理について説明する。これは、学習で求まった分離行列を、他のモジュールに反映させる処理である。
【０３０９】
ステップＳ２５１において、スレッド制御部９１は、スレッドの、観測信号の開始・終了タイミング保持部１６７（図１６参照）に保持されている観測信号の開始タイミングは、観測信号の蓄積タイミング情報保持部１５４（図１５参照）に保持されている、現在の分離行列に対応した蓄積開始タイミングよりも早いか否かを判断する。
【０３１０】
すなわち、図７に示されるように、スレッド１の学習とスレッド２の学習とは、その一部で時間が重なっている。この図では、学習区間７７のほうが、学習区間７８より先に終了しているが、例えば、それぞれの学習にかかる時間によっては、学習区間７７よりも学習区間７８のほうが先に終了してしまう場合もあり得る。
【０３１１】
ここで、ステップＳ２５１の判断が実行されず、学習の終了が遅いものが最新の分離行列として扱われてしまった場合、スレッド２由来の分離行列Ｗ２が、より古いタイミングで取得された観測信号によって学習されて得られたスレッド１由来の分離行列Ｗ１に上書きされてしまう。そこで、新しいタイミングで取得された観測信号によって得られた分離行列が最新の分離行列として扱われるように、観測信号の開始・終了タイミング保持部１６７に保持されている観測信号の開始タイミングと、観測信号の蓄積タイミング情報保持部１５４に保持されている現在の分離行列に対応した蓄積開始タイミングとが比較される。
【０３１２】
ステップＳ２０１において、観測信号の開始タイミングは、現在の分離行列に対応した蓄積開始タイミングよりも早いと判断された場合、換言すれば、このスレッドの学習の結果得られた分離行列Ｗは、現在、観測信号の蓄積タイミング情報保持部１５４に保持されている分離行列Ｗよりも早いタイミングで観測された信号に基づいて学習されていると判断された場合、このスレッドの学習の結果得られた分離行列Ｗは利用されないので、分離行列更新処理は終了する。
【０３１３】
ステップＳ２５１において、観測信号の開始タイミングが、現在の分離行列に対応した蓄積開始タイミングよりも早くないと判断された場合、すなわち、このスレッドの学習の結果得られた分離行列Ｗは、現在、観測信号の蓄積タイミング情報保持部１５４に保持されている分離行Ｗよりも遅いタイミングで観測された信号に基づいて学習されていると判断された場合、ステップＳ２５２において、スレッド制御部１１２は、対応するスレッドの学習によって得られた分離行列Ｗを取得し、分離行列保持部１１３（図１２参照）に供給して設定する。
【０３１４】
ステップＳ２５３において、スレッド制御部１１２は、学習初期値保持部１５２に保持されるそれぞれのスレッドにおける学習の初期値を設定する。
【０３１５】
具体的には、スレッド制御部１１２は、学習初期値として、対応するスレッドの学習によって得られた分離行列Ｗを設定するものとしてもよいし、対応するスレッドの学習によって得られた分離行列Ｗを用いて演算される、分離行列Ｗとは異なる値を設定するものとしても良い。例えば、分離行列保持部１１３（図１２参照）にはリスケーリング適用後の値（後述する図２８のフローのステップＳ２７０の「後処理」より後の値）を代入し、学習初期値保持部１５２にはリスケーリング適用前の値を代入するようにする。それ以外の例については、変形例で説明する。
【０３１６】
ステップＳ２５４において、スレッド制御部１１２は、対応するスレッドの観測信号の開始・終了タイミング保持部１６７（図１６参照）に保持されているタイミング情報を、観測信号の蓄積タイミング情報保持部１５４（図１５参照）に設定し、分離行列更新処理を終了する。
【０３１７】
ステップＳ２５４の処理により、現在使用中、すなわち、分離行列保持部１１３に保持されている分離行列Ｗが、どの時間区間の観測信号から学習されたものであるかが示される。
【０３１８】
このような処理により、新しいタイミングで取得された観測信号によって得られた分離行列が最新の分離行列として設定される。
【０３１９】
次に、図２８のフローチャートを参照して、図２６に示すフローチャートのステップＳ２４１において実行される待機時間の設定処理について説明する。
【０３２０】
ステップＳ２８１において、スレッド制御部１１２は、残りの待機時間を計算する。
【０３２１】
具体的には、スレッド制御部１１２は、残り待機時間（フレーム個数）をｒｅｓｔ、蓄積開始予定タイミング指定情報保持部１５３（図１５参照）に保持されている蓄積開始予定タイミング（フレーム番号、または、対応する相対時刻）をＣｔ、現フレーム番号保持カウンタ１５１に保持されている現フレーム番号をＦｔ、ブロックのシフト幅をｂｌｏｃｋ＿ｓｈｉｆｔとして、残り待機時間ｒｅｓｔを、
ｒｅｓｔ＝Ｃｔ＋ｂｌｏｃｋ＿ｓｈｉｆｔ−Ｆｔ
として算出する。すなわち、Ｃｔ＋ｂｌｏｃｋ＿ｓｈｉｆｔが、次々回蓄積開始予定時刻を意味するため、そこからＦｔを引くことで、「次々回蓄積開始予定時刻までの残り時間」が求まるのである。
【０３２２】
ステップＳ２８２において、スレッド制御部１１２は、残りの待機時間ｒｅｓｔの計算結果は正の値であるか否かを判断する。ステップＳ２８２において、残りの待機時間ｒｅｓｔの計算結果は正の値ではない、すなわち、ゼロまたは負の値であると判断された場合、処理は、後述するステップＳ２８６に進む。
【０３２３】
ステップＳ２８２において、残りの待機時間ｒｅｓｔの計算結果は正の値であると判断された場合、ステップＳ２８３において、スレッド制御部１１２は、「待機中」状態に状態を遷移させるための状態遷移コマンドを、対応するスレッドに発行する。
【０３２４】
ステップＳ２８４において、スレッド制御部１１２は、対応するスレッドのカウンタ１６６（図１６参照）の値を、ｔｈｒｅａｄ＿ｌｅｎ−ｒｅｓｔに設定する。そうすることで、カウンタの値が、ｔｈｒｅａｄ＿ｌｅｎに達するまでの間は、「待機中」状態が継続される。
【０３２５】
ステップＳ２８５において、スレッド制御部１１２は、蓄積開始予定タイミング指定情報保持部１５３（図１５参照）に保持されている値Ｃｔに、ｂｌｏｃｋ＿ｓｈｉｆｔの値を加算する、すなわち、蓄積開始予定タイミング指定情報保持部１５３に次回の蓄積開始タイミングである、Ｃｔ＋ｂｌｏｃｋ＿ｓｈｉｆｔの値を設定し、残り待機時間の計算処理を終了する。
【０３２６】
ステップＳ２８２において、残りの待機時間ｒｅｓｔの計算結果は正の値ではない、すなわち、ゼロまたは負の値であると判断された場合、予定された蓄積開始タイミングを過ぎているのにもかかわらず蓄積が始まっていないことを意味するので、直ちに蓄積を開始する必要がある。そこで、ステップＳ２８６において、スレッド制御部１１２は、「蓄積中」状態に状態を遷移させるための状態遷移コマンドを、対応するスレッドに発行する。
【０３２７】
ステップＳ２８７において、スレッド制御部１１２は、カウンタの値を初期化（例えば０をセット）する。
【０３２８】
ステップＳ２８８において、スレッド制御部１１２は、蓄積開始予定タイミング指定情報保持部１５３に次回の蓄積開始タイミング、すなわち、現フレーム番号であるＦｔを設定し、残り待機時間の計算処理を終了する。
【０３２９】
このような処理により、それぞれのスレッドにおける「学習中」状態にかかる時間に応じて、「待機中」状態とする時間を設定することができる。
【０３３０】
本発明では、スレッド制御部と各学習スレッドとは並列で動作しており、学習スレッドはスレッド制御部とは別のフローチャートに基づいて動いている。以降では、学習スレッドの処理について、図２９のフローチャートを用いて説明する。
【０３３１】
スレッド１１１は、ステップＳ２９１において、初期設定される。そして、ステップＳ２９２において、イベントが発生するまでで待機する（処理をブロックする）。（この「待機」は、学習スレッドの状態の１つである「待機中」とは別である）イベントは、以下のどれかのアクションが行なわれたときに発生する。
・状態遷移コマンドが発行された。
・フレームデータが転送された。
・終了コマンドが発行された。
どのイベントが発生したかによって、以降の処理を分岐する（ステップＳ２９３）。すなわち、スレッド制御部１１２から入力されたイベントによって、続く処理が分岐される。
【０３３２】
ステップＳ２９３において、状態遷移コマンドが入力されたと判断された場合、ステップＳ２９４において対応するコマンド処理が実行される。
【０３３３】
ステップＳ２９３において、フレームデータの転送イベントの入力を受けたと判断された場合、ステップＳ２９５において、スレッド１１１は、フレームデータを取得する。次に、ステップＳ２９６において、スレッド１１１は、取得したフレームデータを、観測信号バッファ１６１（図１６参照）に蓄積し、ステップＳ２９２に戻り、次イベントを待機する。
【０３３４】
観測信号バッファ１６１（図１６参照）は、配列またはスタックの構造をしており、観測信号はカウンタと同じ番号の個所に格納されるものとする。
【０３３５】
ステップＳ２９３において、終了コマンドが入力されたと判断された場合、ステップＳ２９７において、スレッド１１１は、例えば、メモリの開放などの適切な終了前処理を実行し、処理が終了される。
【０３３６】
このような処理により、スレッド制御部１１２の制御に基づいて、それぞれのスレッドにおいて処理が実行される。
【０３３７】
次に、図３０のフローチャートを参照して、図２９に示すフローチャートのステップＳ２９４において実行される、コマンド処理について説明する。
【０３３８】
ステップＳ３０１において、スレッド１１１は、供給された状態遷移コマンドに応じて、それ以降の処理を分岐する。なお、以降では「○○の状態へ遷移する」というコマンドを「状態遷移コマンド『○○』」と表現する。
【０３３９】
ステップＳ３０１において、供給された状態遷移コマンドが、「待機中」状態への遷移を指令する「状態遷移コマンド『待機中』」である場合、ステップＳ３０２において、スレッド１１１は、状態格納部１６５（図１６参照）に、状態が「待機中」であることを示す情報を格納する、すなわち、状態を「待機中」に遷移して、コマンド処理を終了する。
【０３４０】
ステップＳ３０１において、供給された状態遷移コマンドが「蓄積中」状態への遷移を指令する「状態遷移コマンド『蓄積中』」である場合、ステップＳ３０３において、スレッド１１１は、状態格納部１６５に、状態が「蓄積中」であることを示す情報を格納する、すなわち、状態を「蓄積中」に遷移して、コマンド処理を終了する。
【０３４１】
ステップＳ３０１において、供給された状態遷移コマンドが「学習中」状態への遷移を指令する「状態遷移コマンド『学習中』」である場合、ステップＳ３０４において、スレッド１１１は、状態格納部１６５に、状態が「学習中」であることを示す情報を格納する、すなわち、状態を「学習中」に遷移する。
【０３４２】
さらに、ステップＳ３０５において分離行列の学習処理を実行する。この処理の詳細については後述する。
【０３４３】
ステップＳ３０６において、スレッド１１１は、学習が終了したことをスレッド制御部１１２へ通知するために、学習終了フラグ１６８をＯＮにして処理を終了する。フラグを立てることで、学習が終了した直後であることをスレッド制御部１１２へ通知するのである。
【０３４４】
このような処理により、スレッド制御部１１２から供給された状態遷移コマンドに基づいて、それぞれのスレッドの状態が遷移される。
【０３４５】
次に、図３１のフローチャートを参照して、図３０に示すフローチャートのステップＳ３０５において実行される処理の一例である分離行列の学習処理例について説明する。これは、バッチによって分離行列を求める処理であり、バッチ処理であればどのアルゴリズムも適用可能である。ただし、パーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）を起こしにくい方式を用いる必要がある。以下では、出願番号２００５−０１８８２２『音声信号分離装置・雑音除去装置および方法』に開示した構成を適用した例について説明する。
【０３４６】
ステップＳ３１１において、スレッド１１１の学習演算部１６３（図１６参照）は、観測信号バッファ１６１に蓄積された観測信号に対して、必要に応じて、前処理を実行する。
【０３４７】
具体的には、学習演算部１６３は、観測信号バッファ１６１に蓄積された観測信号に対して、正規化（normalization）や無相関化（uncorrelation または ｐｒｅ−ｗｈｉｔｅｎｉｎｇ）などの処理を行なう。例えば、正規化を行なう場合、学習演算部１６３は、観測信号バッファ１６１内の各周波数ｂｉｎについて観測信号の標準偏差を求め、標準偏差の逆数からなる対角行列をＳとして、以下に示す式［９．１］により、Ｘ'＝ＳＸを計算する。
【０３４８】
【数１７】
2010020294.tif 000020
【０３４９】
式Ｘ'＝ＳＸにおいて示されるＸは、観測信号バッファ１６１内の全フレーム分の観測信号からなる行列である。式［９．１］に示すＸ（ｔ）をｔ＝１乃至Ｔについて横方向に並べたもの、すなわち、Ｘ＝［Ｘ（１），・・・Ｘ（Ｔ）］（Ｔは総フレーム数、すなわち、ブロック長ｂｌｏｃｋ＿ｌｅｎ）と表現される。以下の式に現われるＸは、前処理を行なったＸ'も表わし得るものとする。
なお、式［９．２］は、Ｘ'と分離行列Ｗと分離結果Ｙの関係式である。
式［９．３］以下については、後述するステップＳ３２０の後処理におけるリスケーリング処理の説明中において説明する。
【０３５０】
ステップＳ３１２において、学習演算部１６３は、分離行列の初期値として、スレッド制御部１１２から、スレッド制御部１１２の学習初期値保持部１５２に保持されている学習初期値Ｗを取得する。
【０３５１】
ステップＳ３１３〜Ｓ３１９の処理は、学習のループであり、これらの処理をＷが収束するか打ち切りフラグがＯＮになるまで繰り返す。打ち切りフラグは、先に説明した図２６の学習中処理のフローのステップＳ２３６においてＯＮに設定されるフラグである。後から開始した学習がそれより前に開始した学習よりも早く終了した場合にＯＮになるものである。ステップＳ３１３において、打ち切りフラグがＯＮであると判断した場合は、処理を終了する。
【０３５２】
ステップＳ３１３において、打ち切りフラグがＯＦＦであると判断した場合は、ステップＳ３１４に進む。ステップＳ３１４において、学習演算部１６３は、分離行列Ｗの値が収束したか否かを判断する。分離行列Ｗの値が収束したか否かは、例えば、行列のノルムを用いて判定する。分離行列Ｗのノルム（全要素の２乗和）である‖Ｗ‖と、ΔＷのノルムである‖Ｗ‖とをそれぞれ計算し、両者の比である‖ΔＷ‖／‖Ｗ‖が一定の値（例えば１／１０００）よりも小さければ、Ｗが収束したと判定する。または単純に、ループが一定数（例えば５０回）回ったかどうかで判定しても構わない。
【０３５３】
ステップＳ３１４において、分離行列Ｗの値が収束したと判断された場合、処理は、後述するステップＳ３２０に進み、後処理を実行して処理を終了する。すなわち、分離行列Ｗが収束するまで、学習処理ループが実行される。
【０３５４】
ステップＳ３１４において、分離行列Ｗの値が収束していないと判断された場合（またはループ回数が所定の値に達していない場合）、ステップＳ３１５〜Ｓ３１９の学習のループの中に進む。学習は、全ての周波数ビンにおいて、先に説明した式［４．１］〜式［４．３］を繰り返す処理として行なわれる。すなわち、分離行列Ｗを求めるために、式［４．１］から式［４．３］までを分離行列Ｗが収束するまで（または一定回数）繰り返し実行する。この繰り返し実行が「学習」である。なお、分離結果Ｙ（ｔ）は式［４．４］で表わされる。
【０３５５】
ステップＳ３１６が、式［４．１］に対応する。
ステップＳ３１７が、式［４．２］に対応する。
ステップＳ３１８が、式［４．３］に対応している。
【０３５６】
式［４．１］〜式［４．３］は周波数ビンごとの式であるため、ステップＳ３１５とステップＳ３１９で周波数ビンについてのループを回すことで、全周波数ビンのΔＷを求めている。
【０３５７】
これらのループ処理の終了後に、ステップＳ３１３に戻り打ち切りフラグの判定、ステップＳ３１４における分離行列の収束判定を行う。打ち切りフラグがＯＮであれば処理を終了する。ステップＳ３１４で分離行列の収束が確認された場合（あるいは規定ループ数に達した場合）は、ステップＳ３２０に進む。
【０３５８】
ステップＳ３２０では、後処理として以下の処理を実行する。
（１）分離行列を、正規化前の観測信号に対応させる。
（２）周波数ビンの間のバランスを調整する（リスケーリング）。
【０３５９】
まず、（１）分離行列を、正規化前の観測信号に対応させる処理について説明する。
前処理として正規化が行なわれた場合、上述した処理（図３１のステップＳ３１５〜Ｓ３１９）により求められる分離行列Ｗは、正規化前の観測信号Ｘを分離して得られるものと等価ではなく、正規化後の観測信号Ｘ'を分離して得られたものである。そこで、上述した処理により求められた分離行列Ｗを補正して、正規化前の観測信号Ｘを分離して得られるものと等価なものに変換する。具体的には、正規化の際に作用させた行列Ｓを用いて、Ｗ←ＷＳという補正を行なうものとすれば良い。
【０３６０】
次に、（２）周波数ビンの間のバランスを調整する（リスケーリング）処理について説明する。
ＩＣＡのアルゴリズムによっては、分離結果Ｙの周波数ビン間のバランス（スケール）が、予想される原信号のものと異なっている場合がある。そのような場合には、後処理において周波数ビンのスケールを補正する。この補正をリスケーリング（ｒｅ−ｓｃａｌｉｎｇ）という。
【０３６１】
リスケーリングは、例えば、先に観測信号Ｘに対する正規化等の前処理（ステップＳ３１１）について説明した式［９．１］の展開式として示す式［９．３］に含まれる補正用行列Ｒを式［９．６］を用いて求め、それを分離行列Ｗに乗じることで行なう（式［９．６］のｄｉａｇ（・）は、与えられたベクトルを対角要素とする対角行列を表わす）。すなわち、
Ｗ←ＲＷを行なう。ただし、式［９．６］のλｋ（ω）は、式［９．４］右辺の［ａｒｇ  ｍｉｎ］のカッコ内を最小にする値であり、具体的には式［９．５］で計算される。（Ｅｔ［］は、ブロック内の全フレームでの平均を表わす。）
【０３６２】
以上の２つの補正をまとめると、式［９．３］のように表わせる。そのため、両方をあわせた後処理は、Ｗ←ＲＷＳという補正を行なえばよい。
以上で、分離行列学習処理は終了する（２８−１１）。
【０３６３】
なお、ステップＳ３１７においては、自然勾配法（natural gradient）に基づく式として式［４．２］を用いたが、代わりに、「独立性による等分散適応分離」（Equi-variant Adaptive Separation via Independence: EASI）と呼ばれるアルゴリズムに基づいた学習規則である式［４．１３］を用いてもよい。EASI アルゴリズムの式は、勾配法と無相関化とを合わせた形をしており、自然勾配法である式［４．２］よりも収束が速いという特徴がある。EASI アルゴリズム自体についての詳細は、『詳解  独立成分分析』（東京電機大出版局）の「12.5 独立性による等分散適応分離」の節などを参照されたい。
【０３６４】
式［４．１３］にはＥｔ［］を含む項が３つ現われているが、以下に述べる式変形によって計算量を削減することができる。１番目の項であるＥｔ［Ｙ（ω，ｔ）Ｙ（ω，ｔ）Ｈ］は式［４．１４］のように変形でき、その中のＥｔ［Ｘ（ω，ｔ）Ｘ（ω，ｔ）Ｈ］はループ内で一定なので最初に１回だけ計算し（例えばステップＳ３１１で計算する）、ループ内ではＥｔ［Ｘ（ω，ｔ）Ｘ（ω，ｔ）Ｈ］の前後にＷ（ω）およびＷ（ω）Ｈを乗じる演算のみを行えばよい。一方、３番目の項であるＥｔ［Ｙ（ω，ｔ）φ（Ｙ（ｔ））Ｈ］については式［４．１５］のように変形できる。これは、２番目の項を計算した後でエルミート転置を行なえは３番目の項も得られることを表している。
【０３６５】
また、ＥＡＳＩアルゴリズムには、収束時にＥｔ［Ｙ（ω，ｔ）Ｙ（ω，ｔ）Ｈ］＝Ｉが成立するという特徴がある。すなわち、Ｅｔ［Ｙ（ω，ｔ）Ｙ（ω，ｔ）Ｈ］の対角要素についてはＥｔ［｜Ｙｋ（ω，ｔ）｜２］＝１が成り立ち、それ以外の要素については０となる。この性質を利用し、Ｅｔ［｜Ｙｋ（ω，ｔ）｜２］＝１を満たすようにループ内でＷ（ω）のスケールを調整すると、収束速度を一層速くすることができる。具体的には、ステップＳ３１６を実行した後、Ｙ（ω，ｔ）に対しては式［４．１６］の操作を、Ｗ（ω）に対しては式［４．１７］の操作を行なう。ただし、ｄｉａｇ（）は括弧内の要素を対角成分とする対角行列を表わす。このような操作を行なうことで、少ないループ回数でＥｔ［Ｙ（ω，ｔ）Ｙ（ω，ｔ）Ｈ］＝Ｉが成立するようになる。
【０３６６】
次に、図１９に示すフローチャートのステップＳ１０９において実行する音源方向推定処理について、図３２および図３３に示すフローチャートを用いて説明する。この音源方向推定処理は、図１２に示す音源方向推定部１２５において実行される。
【０３６７】
前述したように、本発明では現在時刻を含む区間において、観測信号と分離結果との間で相互共分散行列を計算し、その相互共分散行列を分離行列の代わりに用いて音源方向を計算するという処理を行う。相互共分散行列ΣＸＹ（ω）は先に説明した式［６．１］で表される。
【０３６８】
図３２に示すフローチャートは、先に説明した図１３の処理に対応する。すなわち、相互共分散行列を観測信号Ｘと分離結果Ｙとから計算する（式［６．１］）場合の処理フローである。
図３３に示すフローチャートは、先に説明した図１４の処理に対応する。すなわち、相互共分散行列を観測信号Ｘと分離行列Ｗとから計算する（式［６．５］）場合の処理フローである。
【０３６９】
まず、図３２に示すフローチャートを参照して、相互共分散行列を観測信号と分離結果とから計算する方式を適用した場合の音源方向推定処理シーケンスについて説明する。この処理は図１３を参照して説明した構成に対応する処理シーケンスとなっている。
【０３７０】
まず、ステップＳ３２１において、観測信号と分離結果との共分散行列を更新する。具体的には、フーリエ変換部１２６（図１２参照）から受け取った観測信号Ｘ（ω，ｔ）と、分離部１２７から受け取った分離結果Ｙ（ω，ｔ）、さらに観測信号バッファ１３３（図１３参照）から取り出した観測信号Ｘ（ω，ｔ−Ｌ）、分離結果バッファ１３４から取り出した分離結果Ｙ（ω，ｔ−Ｌ）という４つのベクトルを用いて、先に示した式［６．１０］によって、相互共分散行列ΣＸＹ（ω）を算出する。この処理を全ての周波数ビン（ω＝１，...，Ｍ）について行なう。
【０３７１】
次に、ステップＳ３２２において、観測信号Ｘ（ｔ）を観測信号バッファ１３３（図１３参照）に格納する。その際、Ｘ（ｔ）はそれまでＸ（ｔ−Ｌ）が格納されていた領域に上書きされる。次に、ステップＳ３２３において、分離結果に対しても同様に、分離結果バッファ１３４に格納する。
【０３７２】
次に、
ステップＳ３２４において位相差の計算、
ステップＳ３２５において直線当てはめ処理、
ステップＳ３２６において、音源方向の計算、
これらの処理を行う。このステップＳ３２４からステップＳ３２６の処理は、フレーム毎に行なっても良いが、用途によっては数フレームに１回の頻度で行なっても良い。これらの処理をスキップする場合は音源方向推定処理を終了する。
【０３７３】
ステップＳ３２４の位相差の計算は、例えば、図１０に示すプロットを、全チャンネルおよび全マイクペアについて生成する処理である。詳細は後述する。
ステップＳ３２５における直線当てはめ処理は、例えば図１０のプロットに対して、各点の位相角をアンラップさせながら、最もフィットする直線を探し出す処理である。あるいは、図１１（ｃ），（ｄ）に示すように、プロットに対して、最もフィットするコサインとサインの曲線を探し出す処理でもある。詳細については後述する。
【０３７４】
直線当てはめの結果、例えば図３５に示す位相角プロットから、図４１に示す直線が得られる。図４１には、直線当てはめの結果得られる直線を点線で示してある。なお、図は、直線を分かりやすく示すため、プロットの図から左右または上下にはみだして直線を示している。
【０３７５】
直線当てはめの結果、直線の傾きとして、先に説明した式［８．３］のｈａｔ（Ｋｉｍｋ）が推定される。最後に、ステップＳ３２７において、ｈａｔ（Ｋｉｍｋ）から音源ごとに音源方向を計算する。これも、詳細は後述する。
【０３７６】
次に、音源方向推定処理の別の処理例について、図３３のフローチャートを参照して説明する。この図３３のフローは相互共分散行列を、観測信号の共分散行列と分離行列のエルミート転置とから計算する方式である。この処理は、先に説明した図１４の構成に対応した処理である。
【０３７７】
まず、ステップＳ３３１において、観測信号の共分散行列を更新する。具体的には、フーリエ変換部１２６（図１２参照）から受け取った観測信号Ｘ（ω，ｔ）と、観測信号バッファ１３３（図１３参照）から取り出した、時間的に前の観測信号Ｘ（ω，ｔ−Ｌ）とを用いて、先に示した式［６．１１］の演算によって、共分散行列ΣＸＸ（ω）を算出する。を行なう。この処理を全ての周波数ビン（ω＝１，...，Ｍ）について行なう。
【０３７８】
次に、ステップＳ３３２において、観測信号Ｘ（ｔ）を観測信号バッファ１４３（図１４参照）に格納する。その際、Ｘ（ｔ）はそれまでＸ（ｔ−Ｌ）が格納されていた領域に上書きされる。
【０３７９】
次に、ステップＳ３３３において、相互共分散行列を計算する。具体的には、ステップＳ３３１で更新した共分散行列ΣＸＸ（ω）と、学習処理部１１０の分離行列保持部１１３に保持された最新の分離行列Ｗ（ω）とを用いて、先に示した式［６．５］
の計算によって相互共分散行列ΣＸＹ（ω）を計算する。この処理を全ての周波数ビン（ω＝１，...，Ｍ）について行なう。
【０３８０】
次に、
ステップＳ３３４において位相差の計算、
ステップＳ３３５において直線当てはめ処理、
ステップＳ３３６において、音源方向の計算、
これらの処理を行う。このステップＳ３３４からステップＳ３３６の処理は、フレーム毎に行なっても良いが、用途によっては数フレームに１回の頻度で行なっても良い。これらの処理をスキップする場合は音源方向推定処理を終了する。
【０３８１】
次に、図３２のフローチャートにおけるステップＳ３２４、図３３のフローチャートにおけるステップＳ３３４の「位相差の計算」の詳細について、図３４のフローチャートを参照して説明する。
【０３８２】
ステップＳ３４１において、結果を格納するための領域を必要に応じて確保する。
ステップＳ３４２からステップＳ３４８までは、周波数ビンについてのループであり、周波数ビンωを１からＭまで変化させながら、ステップＳ３４３〜Ｓ３４７を繰り返し処理する。（ただし、並列処理が可能なシステムを使用している場合は、並列に処理してもかまわない。以降のチャンネルループやマイクペアループについても同様である。）
【０３８３】
ステップＳ３４３からステップＳ３４７までは、チャンネル（推定音源、すなわち分離結果）についてのループである。例えばチャンネルを４としたとき、ｃｈａｎｎｅｌ−１からｃｈａｎｎｅｌ−４までの全てについて、繰り返し処理を行なう。
【０３８４】
ステップＳ３４４〜Ｓ３４６までは、マイクのペアについてのループであり、２つのマイクの番号をｉ，ｍとすると、（ｉ，ｍ）＝（１，２）から（ｉ，ｍ）＝（ｎ−１，ｎ）までのｎ（ｎ−１）／２回のループである。このループは、図３４（Ｂ）に示すステップＳ３５１〜３５４のような２重のループで表現することも可能である。
【０３８５】
ループの内部では、相互共分散行列ΣＸＹの要素σｉｋ（ω）とσｍｋ（ω）との要素比を先に示した式［８．１］によって計算し、ステップＳ３４１で確保した領域に格納する。必要に応じて、要素比の外にその要素比の位相角も、式［８．９］に従って計算する。これらの操作を、全マイクペア・全チャンネル・全周波数ビンについて行なえば、位相差の計算は完了する。
【０３８６】
図３４のフローチャートでは、マイクペアのループは全マイクの組み合わせについて行なっていたが、代わりに、限定されたマイクペアについてのみ行なってもよい。例えば、（ｉ，ｍ）＝（１，２），（２，３），…，（ｎ−１，ｎ）のｎ−１通りのペアや、（１，２），（１，３），…，（１，ｎ）のｎ−１通りのペアに対してのみ行なってもよい。また、複数のマイクの中に、位置や間隔の既知のものと未知のものとが混在している場合は、マイク間の距離が既知であるペアについてのみループを回してもよい。以降のフローチャートに現れる「マイクペアループ」についても同様に、限定されたペアについてのみ行なってもよい。
【０３８７】
このフローで説明した位相差の計算は、例えば、図１０に示すプロットを、全チャンネルおよび全マイクペアについて生成する処理である。例えば、マイク数４、チャンネル数４とした場合、図３４に示すフローの処理結果として、図３５に示すデータが生成される。
【０３８８】
図３５に示すデータは、マイク１〜４の４つのマイクから２つを選んだ６つの組み合わせ（マイク１＆２，１＆３，１＆４，２＆３，２＆４，３＆４）における複数の推定音源である各チャンネル１〜４のプロットを設定した２４個の図を示している。これらの各図は、相互共分散行列から計算したｒｉｍｋ（ω）について、位相角をプロットしたものであり、縦軸が周波数ビンωを、横軸が位相角を表わす。点がほぼ直線状に並んでいる。なお、先に図１０を参照して説明したように、位相角ａｎｇｌｅ（）は２πの範囲（この図では−π〜π）でしか値を特定できないため、−πとπとの間で値の「飛び」が存在する。
【０３８９】
次に、図３２のフローチャートにおけるステップＳ３２５、図３３のフローチャートにおけるステップＳ３３５の「直線当てはめ」の詳細について、図３６のフローチャートを参照して説明する。
【０３９０】
なお、図３６に示すフローチャートは概略であり、プロット一枚分（例えば、図３５に示す１つのプロットデータ３５１）に対してのステップＳ３６４の当てはめ処理を、ステップＳ３６３〜Ｓ３６５で全マイクペア（図３５の縦方向に相当）に対して行い、さらに、ステップＳ３６２〜Ｓ３６６で全チャンネル（図３５の横方向に相当）について行なうことを表わしている。
【０３９１】
ステップＳ３６４の当てはめ処理を行なうアルゴリズムはさまざまなものが考えられるが、ここでは以下の３種類について説明する。
（１）再帰法
（２）リニアサーチ法
（３）ハイブリッド法
以下は、各手法の説明である。
【０３９２】
（１）再帰法
再帰法とは、低い周波数ビンから順に、先に示した式［８．９］によって表現される位相角αｉｍｋ（ω）をアンラップしながら傾きＫｉｍｋを求めていく方式である。再帰法の詳細を、図３７および図３８を参照して説明する。
【０３９３】
図３７は、再帰法によって、（直線）当てはめ処理を行うシーケンスを説明するフローチャートである。図３７に示すフローにおいて、Ｋｉｍｋは最終的な傾きだけでなく、処理途中の傾きの値を格納する変数でもある。ステップＳ３７１において、傾きＫｉｍｋに適切な初期値を代入しておく。
【０３９４】
初期値は０でもよく、これは図１０や図３５では垂直に立った直線に相当する。また、前フレームで求めた傾きが保存されている場合は、それを用いても良い。
ステップＳ３７２〜Ｓ３７６は、周波数ビンのループである。ただし、直流成分であるω＝１の位相角は０と決まっているため、ω＝２から開始する。
【０３９５】
ステップＳ３７３〜Ｓ３７５は、再帰的な処理である。すなわち、２からω−１までの周波数ビンの位相角を用いてＫｉｍｋの値（ただし途中段階）が求まっている状態において、ω番目の位相角をアンラップしてからＫｉｍｋの値を更新している。
【０３９６】
ステップＳ３７３は、２からω−１までの周波数ビンにおいて当てはめた直線（図３８の直線４０１）を、ω番目の周波数ビンへ延長する（図３８のポイントＰ４０２）ことを意味する。延長した先の点を大文字のＰ、その位相角を小文字のｐで表わす。
【０３９７】
ステップＳ３７４のアンラップとは、先に図１０を参照して説明した図１０（ａ）から図１０（ｂ）への変換処理である。すなわち、周波数ビンωの位相角αｉｍｋ（ω）に２Ｎπ（Ｎは整数）を加えて点Ｐに最も近づけることである。具体的には、以下に示す式［１０．１］を満たすＮを求めてから式［１０．２］の操作を行なえばよい。
【０３９８】
【数１８】
2010020294.tif 000021
【０３９９】
図３８では、アンラップ前の位相角αｉｍｋ（ω）を点［αｉｍｋ（ω）（旧）］４０３、アンラップ後の位相角を点［αｉｍｋ（ω）（新）］４０４として示している。
【０４００】
次に、ステップＳ３７５において、アンラップ後の位相角αｉｍｋ（ω）の値を反映させて傾きＫｉｍｋを更新する（図３８の直線４０５）。この更新のためには、例えば、先に示した式［１０．３］〜式［１０．５］のどれかを使用する。
【０４０１】
式［１０．３］は、アンラップ済みである位相角αｉｍｋ（２）〜αｉｍｋ（ω）それぞれについて傾きを計算し、平均を取ることを意味する。ただしω＝２の場合、Σの部分は０とする。
【０４０２】
式［１０．３］のΣの部分は、更新前の傾きＫｉｍｋを用いて、
（ω−２）×Ｋｉｍｋ
と近似することもでき、その結果、式［１０．４］が得られる。式［１０．４］の実装においてはシグマ（総和）のループが不要なので、式［１０．３］よりも計算量が少ない。
【０４０３】
式［１０．３］および式［１０．４］は共に、ωが小さいときにはＫｉｍｋの値が大きく変化し、ωが大きくなるにつれてＫｉｍｋが変化しにくくなるという特徴がある。従って、ステップＳ３７１で代入した初期値は最終的なＫｉｍｋの値にあまり反映されない。
【０４０４】
そこで、初期値を最終的な値に強く反映させたい場合（例えば、（３）のハイブリッド方式の一部として再帰法を用いる場合や、前フレームで求めた傾きを初期値として用いる場合など）は、式［１０．４］の代わりに式［１０．５］を用いる。式［１０．５］において、βは１より大きな値であり、βを大きくするほど傾きＫｉｍｋは初期値から変化しにくくなる。
【０４０５】
ステップＳ３７２〜Ｓ３７６をω＝Ｍ（Ｍは周波数ビンの総数）になるまで繰り返すと、位相角のアンラップが完了し、そのときの傾きＫｉｍｋの値が、当てはめ直線の傾きとなっている。なお、周波数ビンのループは途中で打ち切っても良く、その分だけ計算量が削減できる。例えば、音源が全て人間の音声であると分かっているような状況では、音声の周波数帯域に対応した周波数ビンの範囲でのみループを回せばよい。
【０４０６】
ループを抜け出したら、ステップＳ３７７において、先に示した式［８．５］に従って当てはめ誤差を計算する。この値は、先に示した図３３のフローのステップＳ３３６の音源方向の計算などで使用される。
【０４０７】
再帰法は、
１）最大でＭ回の反復で処理が完了する、
２）「低い周波数ビンの位相角には『飛び』が発生しない」という特徴を有効に利用している。
といった利点がある。しかし一方で、低い周波数ビンの位相角に誤差が多い（直線から離れている点が多い）場合や、最終的な傾きが初期値と大きく離れている場合などに、直線当てはめを失敗する（無関係の方向に直線を引いてしまう）可能性もある。
【０４０８】
（２）リニアサーチ法
次に直線当てはめのもう一つのアルゴリズムであるリニアサーチについて、図３９に示すフローチャートを参照して説明する。
【０４０９】
この処理は、傾きＫｉｍｋの値を、ステップＳ３９１で設定する最小値Ｋｍｉｎから、予め設定した最大値Ｋｍａｘまで少しずつ変化させながら、先に示した式［８．５］によって当てはめ誤差を計算し、その誤差を最小にする傾きＫｉｍｋおよび最小値を求めるというものである。
【０４１０】
傾き設定範囲の最大値であるＫｍａｘは、以下に示す式［１１．１］で求める。
【０４１１】
【数１９】
2010020294.tif 000022
【０４１２】
これは、ω=2の周波数ビンにおいて、音がマイクペアの真横（θ＝π／２）から到達した場合の位相差と等しい値である。ｄの値は、マイクｉとマイクｍとの距離ｄｉｍと一致させても良いが、マイク間距離の最大値を共通に用いても良い。また、マイク間距離そのものを用いる代わりに、若干大きめの値（例えば１．１倍した値）を用いても良い。
【０４１３】
傾き設定範囲の最小値であるＫｍｉｎについては、
Ｋｍｉｎ＝−Ｋｍａｘとすればよい。
これはθ＝−π／２の方向から音が到達した場合の位相差に対応する値である。
【０４１４】
Ｋｉｍｋの増分であるΔＫ（ステップＳ３８６参照）は、比較的小さな正の値であり、例えば、
（Ｋｍａｘ−Ｋｍｉｎ）／１０００
といった値を用いる。
【０４１５】
まず、ステップＳ３８１において、変数の初期値として、
当てはめ誤差ｅｒｒｉｍｋ＝無限大、
傾きＫ＝Ｋｍｉｎ
とおく。
【０４１６】
ステップＳ３８２において、Ｋが最大値（Ｋｍａｘ）を超えていないかを確認し、超えていた場合は処理を終了する。ステップＳ３８２において、Ｋが最大値（Ｋｍａｘ）を超えていない場合は、ステップＳ３８３において、設定した角度において検出される当てはめ誤差ｅｒｒを算出する。先に示した式［８．５］によって当てはめ誤差を計算する。
【０４１７】
ステップＳ３８４において、算出した当てはめ誤差ｅｒｒと、保持されている当てはめ誤差ｅｒｒｉｍｋを比較し、算出した当てはめ誤差ｅｒｒが保持されている当てはめ誤差ｅｒｒｉｍｋをより大きい場合は、ステップＳ３８６に進み、角度を更新する。
【０４１８】
算出した当てはめ誤差ｅｒｒが保持されている当てはめ誤差ｅｒｒｉｍｋをより大きくない場合は、ステップＳ３８５に進み、変数としての角度Ｋｉｍｋを、検証した角度Ｋによって更新し、当てはめ誤差ｅｒｒｉｍｋを算出した誤差ｅｒｒによって更新する。その後、ステップＳ３８６に進み、検証角度を更新する。さらに、更新した角度Ｋについて、ステップＳ３８２以下の処理を繰り返す。
【０４１９】
このリニアサーチ法においては、計算量と最適解の精度との間にトレードオフがあ
る。すなわち、ΔＫの値として０に近いものを用いるほど、Ｋｉｍｋの精度は良くなるが、その反面計算量は増える。
【０４２０】
（３）ハイブリッド法
次に、再帰法とリニアサーチ法とを組み合わせたハイブリッド法による直線当てはめ処理について説明する。このハイブリッド処理は、最初に比較的大きなΔＫを用いて大雑把なリニアサーチを行ない、そこで求まった比例係数を初期値として、次に再帰法を行なう。
【０４２１】
ハイブリッド法では、図３６のフローにおけるステップＳ３６４の当てはめ処理が大きく２つのステップに分かれており、最初に図３９を参照して説明したリニアサーチを行ない、次にそこで求まった傾きＫｉｍｋを初期値として、図３７のフローで示す再帰法を行なう。以下、このハイブリット法と、リニアサーチ法・再帰法をそれぞれ単独に実行する場合とに違いについて説明する。
【０４２２】
ハイブリッド方では、傾きの設定範囲の最大値Ｋｍａｘの値は、以下に示す式［１２．１］および式［１２．２］から求める。
【０４２３】
【数２０】
2010020294.tif 000023
【０４２４】
上記式［１２．１］のｆｌｏｏｒ（）は小数点以下の切捨て（ｆｌｏｏｒｉｎｇ）を表わす。１／２を加えてから切捨てを行なっているため、ｄＦ／Ｃを小数点以下で四捨五入しているのと等価である。
【０４２５】
最大マイク間隔ｄ＝０．１５［ｍ］、
サンプリング周波数Ｆ＝１６０００［Ｈｚ］、
音速Ｃ＝３４０［ｍ／ｓ］、
この設定の場合、Ｎ＝７が得られる。
【０４２６】
Ｋｍｉｎについては、先に説明した（２）リニアサーチ法と同様、
Ｋｍｉｎ＝−Ｋｍａｘ
として設定する。
【０４２７】
次に、リニアサーチ法のフローとして説明した図３９のステップＳ３９６において適用する傾きの増分［ΔＫ］については、上記の式［１２．３］で計算する。すなわち、−πＮ／（Ｍ−１），−π（Ｎ−１）／（Ｍ−１），...，０，...，π（Ｎ−１）／（Ｍ−１），πＮ／（Ｍ−１）という２Ｎ＋１種類の傾きの間で、最適値を探すのである。
【０４２８】
例えば、設定する２Ｎ＋１種類の傾きは、図４０の直線４１１〜４１７の各直線にそれぞれ対応する。傾きの増分の算出式は以下のように示される。
ΔＫ＝π／（Ｍ−１）
この値は重要であり、ΔＫがこの値よりも大きい場合は、以降の再帰法で真の最適値を見つけ損なう可能性がある。ΔＫがこの値よりも小さい場合は、リニアサーチにおいて冗長な処理が発生する。以下では、その理由について説明する。位相差プロットのアンラップは式［１０．１］および式［１０．２］に基づいて行われる。式［１０．１］のＮは整数であるため、傾きＫｉｍｋの値が±π／（ω−１）の範囲で変化してもＮとしては同一の値となり、同一のアンラップ結果が得られる。すなわち、傾きＫｉｍｋの値を、π／（Ｍ−１）よりも小さな値で変化させても、全周波数ビンにおいて変化前と同一のアンラップ結果しか得られない。言い換えると、１段階目のサーチにおけるΔＫ＝π／（Ｍ−１）という値は、２段階目のサーチ（再起法）においてアンラップ結果に違いを出すための最小の値である。
【０４２９】
上記の傾きの増分［ΔＫ＝π／（Ｍ−１）］を設定した上で、図３９を参照して説明したリニアサーチを行ない、次にそこで求まった傾きＫｉｍｋを初期値として、図３７のフローで示す再帰法を行なう。
【０４３０】
この再帰法の処理において、図３７に示すステップＳ３７１において、上記の２Ｎ＋１種類の傾きの中の最適値を初期値［Ｋｉｍｋ］設定してステップＳ３７２以下の処理を実行する。
【０４３１】
ステップＳ３７５においては、Ｋｉｍｋが初期値から大きく外れるのを防ぐため、先に示した傾き［Ｋｉｍｋ］算出式である式［１０．３］や式［１０．４］の代わりに式［１０．５］用いる。ただし、式［１０．５］のβとして、周波数ビン数Ｍと同程度の値を用いる。あるいは、先に示した式［１２．４］を用いても良い。
【０４３２】
このような再帰法とリニアサーチ法とを組み合わせたハイブリッド法による直線当てはめ処理により、効率的にかつ精度の高い直線当てはめが実現される。この直線当てはめを行った結果を図４１に示す。図４１は、先に説明した図３５のプロットデータを利用して直線当てはめを行った結果である。
【０４３３】
図４１に示すデータは、マイク１〜４の４つのマイクから２つを選んだ６つの組み合わせ（マイク１＆２，１＆３，１＆４，２＆３，２＆４，３＆４）における複数の推定音源である各チャンネル１〜４のプロットを設定した２４個の図に対する直線あてはめ結果である。プロットは、相互共分散行列から計算したｒｉｍｋ（ω）について、位相角をプロットしたものであり、縦軸が周波数ビンωを、横軸が位相角を表わす。図４１には、直線当てはめの結果得られる直線を点線で示してある。なお、図は、直線を分かりやすく示すため、プロットの図から左右または上下にはみだして直線を示している。
【０４３４】
次に、図３２に示すフローチャートにおけるステップＳ３２６の音源方向の計算、図３３に示すフローチャートにおけるステップＳ３３６の音源方向の計算、これらの音源方向の計算について、図４２に示すフローチャートを用いて説明する。
【０４３５】
図４２に示すフローチャートは、マイクペア毎に求めた傾きＫｉｍｋを総合して音源の方向θｋを算出する処理である。ただし、音源に対応していない出力チャンネルや、直線当てはめに失敗しているマイクペアについては、音源方向に値が反映されないように工夫されている。
【０４３６】
ここで、処理の説明に入る前に、「音源に対応していない出力チャンネル」や「直線当てはめに失敗しているマイクペア」について説明する。
【０４３７】
同時に鳴っている音源の数がマイク数よりも少ない場合、出力チャンネルのいくつかは、どの音源にも対応しない信号が出力される。（具体的には、残響成分のみのような音が出力される。）そのような出力チャンネルは、位相角のプロットが直線状には並ばないため、容易に判別できる。
【０４３８】
例えば、先に直線当てはめの結果として説明した図４１はマイク数＝４、チャンネル数（推定音源数）＝４として生成したプロットおよび直線当てはめ結果である。しかし、図から理解されるように、チャンネル３に示すプロットデータは、直線からはずれているものが極めて多い。このデータは、この推定音源（チャンネル３）が実際には存在していないものであることを示している。すなわち、同時音源数＝３の場合の結果であることを示している。
【０４３９】
すなわち、このｃｈａｎｎｅｌ−３の列がどの音源にも対応しない出力である。このような出力チャンネルに対して音源方向を計算しても無意味であるため、このようなチャンネルを検出して棄却する必要がある。また、特定のマイクペアについてのみ直線当てはめに失敗している場合、そのペアを除外して音源方向の平均を求めた方が、より正確な値が出るはずである。
【０４４０】
図４２に示す音源方向の計算処理フローでは、先に示した式［８．５］から求めた当てはめ誤差を用いて誤差の大きなペアは音源方向に反映されないようにすることで、音源に対応しないチャンネルや直線当てはめに失敗したマイクペアなどに由来する音源方向の値を除外するようにしている。
【０４４１】
図４２に示す音源方向の計算処理フローにおけるステップＳ４０１〜Ｓ４１０はチャンネルについてのループであり、例えば図４１の各列（チャンネル）について以降の処理を行なうことを表わす。
【０４４２】
ステップＳ４０２では、結果を格納する変数に初期値を代入している。
Ｎは、有効なペアの個数、
θｋは音源方向（角度）の総和または平均、
これらを格納する変数であり、初期値は共に０である。
【０４４３】
ステップＳ４０３〜Ｓ４０６はマイクペアについてのループであり、これは図４１の各列の６枚の位相角プロットに対して処理を行なうことを意味する。すなわち、６枚のプロットそれぞれから音源方向を計算し、それらの値を平均するのである。
【０４４４】
ステップＳ４０４の当てはめ誤差ｅｒｒｉｍｋは、先に説明した図３６のステップＳ３６４の当てはめ処理においてマイクペアｉ，ｍとチャンネルｋとの組について計算された当てはめ誤差である。この値を所定の閾値と比較し、閾値を上回っている場合には、ステップＳ４０５の処理をスキップする。
【０４４５】
当てはめ誤差ｅｒｒｉｍｋの計算に式［８．５］を用いている場合は、閾値として０から４までのどれかの値（例えば１）を設定しておく。閾値を適切に設定することにより、直線当てはめに成功している（位相差の各点がもともとほぼ直線上に並んでおり、さらに、それらの点にフィットする直線が適切に求められている）プロットのみが音源方向に反映され、それ以外のプロット（例えば図４１のチャンネル３（ｃｈａｎｎｅｌ−３）の列のプロット）は採用されないようになる。
【０４４６】
ステップＳ４０４において、当てはめ誤差ｅｒｒｉｍｋの値が閾値以下であると判定した場合は、ステップＳ４０５に進む。ステップＳ４０５においては、式［８．７］に基づいてθｉｍｋを計算し、その値をθｋに加算する。また、個数Ｎをインクリメント、すなわち＋１する。なお、式［８．７］をそのまま適用する代わりに、以下の補正を施した方が良い場合もある。その点について説明する。
【０４４７】
式［８．７］の演算において、ａｓｉｎの括弧の中の値が−１〜１の範囲を超える場合がある。その場合、ａｓｉｎを式［８．８］のように拡張し、範囲を超える入力を扱えるようにする。例えば、１を若干越える入力に対しては、π／２を若干超える値が返される。
【０４４８】
また、マイクの配置が一直線状以外である場合は、マイクペアの角度の分だけ式［８．７］の値を補正する。角度の補正について、図４３を参照して説明する。
【０４４９】
図４３において、マイク１（４２１）からマイク２（４２２）の方向を角度の基準とし、マイク１からマイク３（４２３）の方向は基準に対してγだけ回転しているとする（マイクペアの回転角４２７）。この場合、式［８．７］においてｉ＝１，ｍ＝３として求めた値θ１３ｋに対して、マイクペアの回転角γを加えることで、角度の基準を合わせる。すなわち、式［８．７］でθ１３ｋを求めた後、
θ１３ｋ←θ１３ｋ＋γ
という補正を行ない、その補正後の値をステップＳ４０５においてθｋに加算する。
【０４５０】
全てのマイクペアについてステップＳ４０４，Ｓ４０５の処理を行った後では、θｋには角度の総和が格納されている。そこで、総和を個数Ｎで割ることで平均に変換する（ステップＳ４０８）。ただし、Ｎが０の場合は、θｋには無効を表わす値を代入する。
【０４５１】
例えば図４１のチャンネル３（ｃｈａｎｎｅｌ−３）の列は、６枚のプロットが全てスキップされ、結果としてθ３には無効な値がされる。ステップＳ４０１〜Ｓ４１０を全てのチャンネルについて行なうことで、音源方向の計算は完了する。
以上で、図１９のフローチャートに示された全ての処理について、説明を完了する。
【０４５２】
［変形例］
本発明で用いている「直線当てはめ」の結果は、音源方向の他にも、各種用途に利用可能である。ここでは変形例として、次の３点について説明する。
１．周波数マスキング
２．学習初期値の再構築
３．マイク配置が直線的でない場合
以下、それぞれについて説明する。
【０４５３】
［１．周波数マスキング］
音源分離技術の一つに、周波数マスキングと呼ばれる方法がある。これは、大部分の種類の音が時間周波数領域ではスパースであるという性質を利用し、目的外の音の成分をマスクすることで、目的の音のみを残すという方法である。独立成分分析と周波数マスキングとの組み合わせも存在し、例えば特開２００７−２１９４７９に示された構成がそのような例である。
【０４５４】
周波数マスキングには、ハードマスクという方法（目的音以外の信号を０にする）とソフトマスクという方法（目的音以外の信号を、ある尺度に応じて弱くする）とがあり、特開２００７−２１９４７９に示された構成はハードマスクの例である。ソフトマスクは、ハードマスクと比べた場合に、目的音以外の音を除去する性能は劣るが、特定の周波数成分が欠落することによる問題（ミュージカルノイズの発生や音声認識への悪影響など）をある程度は防ぐことができる。
【０４５５】
先に説明した直線当てはめの結果を利用することにより、ソフトマスクの係数を決定することができ、目的音以外の音を一層除去することができる。
【０４５６】
図４４は、あるチャンネルのあるマイクペアにおける直線当てはめの結果であり、先に説明した図１０の直線当てはめ結果と同様の図である。縦軸が周波数ビンωを、横軸が位相角を表わす。点がほぼ直線状に並んでいることが分かる。この図のω番目の周波数ビンにおいて、位相角はαｉｍｋ（ω）５０１、直線上の点は（ω−１）Ｋｉｍｋ５０２である。両者の位相差をψ５０３とおく。一部の周波数ビンでのみψの絶対値が大きい場合、その周波数ビンにおいては分離が失敗している可能性が高いと考えることができる。
【０４５７】
そこで、位相差ψの絶対値に応じて分離結果の周波数成分を減衰させることにより、分離の精度を向上させることができる。
【０４５８】
具体的には、図４５に示すような、位相差ψの絶対値が大きくなるにつれて値が小さくなるような関数を用意する。例えば、以下に示す式［１３．１］のような関数である。
【０４５９】
【数２１】
2010020294.tif 000024
【０４６０】
この式において、
ｈはマスキングのかかり具合を調整する項であり、０以上１以下の値をとる（ｈ＝１の場合は、周波数マスキングを行なわないことを意味する）。ｐは０以上の値をとり、マスクの鋭さを表わす。（図４５はp=2 の場合である。）  pとして大きな値を用いるほど、鋭いマスク、すなわち位相差ψが０から離れるにつれて急激に減衰する形状のマスクを表現することができる。一方、ｐ＝０の場合は、周波数マスキングを行なわないことを意味する。
【０４６１】
直線当てはめを再帰法（またはハイブリッド法）で行なった場合、式［１３．１］の値は、アンラップ済みのαｉｍｋ（ω）と（ω−１）Ｋｉｍｋとを用いて式［１３．２］のように計算することができる。一方、直線当てはめをリニアサーチ法で行なった場合は、代わりに式［１３．３］によって計算することができる。（Ｒｅ［］は、カッコ内の複素数から実部を取り出す操作を表わす。）
【０４６２】
式［１３．１］で求めた値に対して、全マイクペアの間で平均を取る（式［１３．４］）。なお、全マイクペア間で平均を取る代わりに、直線当てはめに成功しているマイクペアについてのみ平均をとっても構わない。すなわち、図４３のフローチャートのステップＳ４０３〜Ｓ４０６の処理と同様に、当てはめ誤差ｅｒｒｉｍｋが閾値以下の場合のみ加算を行ない、ループ終了後に個数での割り算を行なう。
【０４６３】
式［１３．４］（または、当てはめ誤差ｅｒｒｉｍｋが閾値以下のマイクペアについての平均値）の値が求まったら、その値を、チャンネルｋ・周波数ビンω・フレーム番号ｔにおける分離結果Ｙｋ（ω，ｔ）に乗じ（式［１３．５］）、ソフトマスク適用結果であるＹ'ｋ（ω，ｔ）を生成する。
【０４６４】
先に説明した図１９のフローにおけるステップＳ１０７の「逆フーリエ変換」においては、分離結果Ｙｋ（ω，ｔ）の代わりにソフトマスク適用結果であるＹ'ｋ（ω，ｔ）を用いることにより、ソフトマスクが適用された分離結果波形が生成される。そして図１９のフローにおけるステップＳ１０８では、そのソフトマスクが適用された分離結果波形が送信される。
【０４６５】
なお、図４５および式［１３．１］では、マスク値が最小となるのは位相差ψ＝±πの場合のみであるが、代わりに、位相差ψが一定範囲を外れた場合は強制的に最小値となるようなマスクを用いてもよい。具体的には、０＜ψ０≦πを満たすψ０を用意し、−ψ０≦ψ≦ψ０の場合は式［１３．１］でｇｉｍｋ（ψ）を計算し、それ以外の場合はｇｉｍｋ（ψ）＝ｈとすればよい。
【０４６６】
このように、本発明の一実施例に係る信号処理装置は、直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して、直線からのずれの大きい位相角データを持つ信号分離結果の周波数成分を減衰させる周波数マスキングを実行する構成とすることが可能であり、この処理により信号分離性能をさらに向上させることが可能となる。
【０４６７】
［２．学習初期値の再構築］
先に学習処理部におけるスレッド処理の説明として説明した図７の説明では、学習開始時の初期値は、その時点で最新の分離行列（またはリスケーリング前の分離行列）をそのまま用いていた。例えば、図７に示すスレッド１の最初の学習で生成された分離行列７２は、スレッド３の学習開始時７５に初期値として使用され、同じく、スレッド３で生成された分離行列は、スレッド５の学習開始時に初期値として使用されている。複数の学習スレッドを介することによって、図８（ａ）のグラフ１のように収束の度合いを上げていくのである。
【０４６８】
このように初期値として分離行列をそのまま用いる代わりに、直線当てはめの結果から分離行列を再構築するようにすると、収束の速度を速くすることができる。その具体的な方法について、以下で説明する。
【０４６９】
分離行列Ｗ（ω）から直線当てはめの結果であるＫｉｍｋまでは、先に示した式［６．５］→式［７．１］→式［８．１］→式［８．３］という過程を経て推定していた。従って、この逆の操作を行なえば、分離行列を計算することができる。
【０４７０】
以下に示す式［１４．１］は、式［８．１］のｓｉｇｎ（σｉｋ（ω））に相当する要素を計算する式である。
【０４７１】
【数２２】
2010020294.tif 000025
【０４７２】
ただし、ｉ≧ｍの場合のＫｉｍｋの値については、式［１４．２］および式［１４．３］のように定義する。（例えば、Ｋ１１ｋ＝０、Ｋ２１ｋ＝−Ｋ１２ｋである。）
【０４７３】
なお、式［１４．１］においてはｎ個の項の間で平均を取っているが、図４３に示すフローのステップＳ４０３〜Ｓ４０６と同様に、当てはめ誤差ｅｒｒｉｍｋが閾値以下となる（ｉ，ｍ）の組み合わせについてのみ平均をとっても良い。（ただし、ｅｒｒｍｉｋ＝ｅｒｒｉｍｋとする。）あるいは、複数のｍについて平均をとる代わりに、特定のｍについての値のみを用いて計算量を削減してもよい。例えばｍ＝１のみを用いることは、１番目のマイクと他のマイクとの位相差のみを用いて分離行列を再構築することを意味する。
【０４７４】
次に式［１４．４］のように、式［１４．１］のｑｉｋ（ω）を要素とする行列Ｑ（ω）を生成する。Ｑ（ω）は、相互共分散行列ΣＸＹ（ω）に相当する行列である。Ｑ（ω）からＷ（ω）を得るためには、先に示した式［６．５］の関係を用いればよく、その結果、式［１４．５］を得る。あるいは式［６．８］の関係を用いてもよく、あるいは単にＱ（ω）の逆行列を用いてもよい。
【０４７５】
式［１４．５］の演算結果をＵ（ω）とおき、このＵ（ω）を初期値として、図１５に示すスレッド制御部１１２の学習初期値保持部１５２に保持する学習初期値とする。
【０４７６】
なお、Ｕ（ω）を直接初期値として用いる代わりに、学習初期値保持部１５２に保持する行列にスケールを近づけた行列を用いた方が、収束速度の一層の向上が期待できる。そのためには、式［１４．６］の演算を行なえばよい。ただし、Ｗｉｋ（ω）は学習初期値保持部１５２に保持する学習初期値の各要素であり、ｄｉａｇ（）はカッコ内の要素を対角要素に持つ対角行列を表わす。Ｕ（ω）の左側に対角行列を乗じているため、この式は、Ｕ（ω）の各行のスケールを学習初期値保持部１５２に保持する学習初期値の各行のスケールに一致させることを意味する。
【０４７７】
この変形例を施すことによって、収束までの速度が向上する。例えば図８の（ａ）グラフ１では、音源の状態が変化してから分離行列が再び収束するまでの間に８５ｃ〜８５ｆという移行段階が存在していたが、この段数が減る可能性がある。
【０４７８】
このように、本発明の信号処理装置の一実施例としては、直線当てはめ処理の結果として取得される比例定数（直線の傾き）を利用して分離行列を生成し、生成した分離行列を学習処理部において利用する構成とすることができる。この構成により、学習処理における収束を早めることが可能となる。
【０４７９】
［３．マイク配置が直線的でない場合］
次に、マイクの配置が同一直線上ではなく、２次元的（同一平面上）または３次元的に配置されている場合の処理例について説明する。
【０４８０】
マイクの配置が同一直線上ではなく、２次元的（同一平面上）または３次元的に配置されている場合は、音源方向として、角度θｋの代わりに音源方向のベクトルを推定することも可能である。以下では、その原理と方法について説明する。
【０４８１】
図４６は、マイクペアの一つと音源方向ベクトルとの関係を表わした図である。音源６０１はｋ番目のチャンネルに対応した音源であり、それぞれのマイクペアから見た音源ｋの方向をｑｋという単位ベクトル（大きさが１のベクトル）で表わす（音源方向ベクトル６０２）。
【０４８２】
音源方向ベクトルｑｋ６０２は、以下に示す式［１５．１］で表わされる３次元の列ベクトルである。
【０４８３】
【数２３】
2010020294.tif 000026
【０４８４】
また、図４６に現れるベクトルは全て同一の形式（３次元の列ベクトル）で表わされる。図４６に示す経路差［Ｄｉｍｋ］６０６は、音源ｋ６０１からマイクｉ６０３までの距離と、音源ｋ６０１からマイクｍ６０４までの距離との差であり、これは図１に示す経路差ｄｓｉｎθに相当する量である。
【０４８５】
図４６における経路差［Ｄｉｍｋ］６０６と方向ベクトルｑｋ６０２との関係、および、Ｄｉｍｋと直線当てはめ結果Ｋｉｍｋとの関係をそれぞれ式で表わすことで、ｑｋを推定することができる。
【０４８６】
直線当てはめ結果Ｋｉｍｋと経路差Ｄｉｍｋとの関係は、式［１５．４］のように表わせる。Ｆはサンプリング周波数を、Ｍは周波数ビン数、Ｃは音速を表わす。この式は、式［８．４］のｄｉｍ  ｓｉｎθｉｍｋという項をＤｉｍｋに置き換え、さらに"Ｄｉｍｋ＝"
の形に変形したものと等価である。
【０４８７】
一方、平面波を仮定すると、経路差Ｄｉｍｋは、マイクペア方向のベクトル（マイクｉからマイクｍ向かうベクトル）をｑｋに射影したものと考えることができる。すなわち、マイクペア方向をｐｉｍというベクトルで表わす（マイクペアベクトル６０５）と、Ｄｉｍｋ，ｐｉｍ，ｑｋの関係は式［１５．５］で表わすことができる。（上付きのＴは、ベクトルや行列の転置を表わす。ｘ，ｙを列ベクトルとすると、ｘＴｙはベクトルの内積である。さらに式［１５．５］は、内積の値をマイク間距離｜ｐｉｍ｜で正規化してある。）なお、マイクペアベクトルｐｉｍは、マイクｉとマイクｍの位置ベクトル（＝座標）をそれぞれｚｉ，ｚｍとして式［１５．３］で算出する。
【０４８８】
平面波を仮定すると、ｑｋはどのマイクペアから見ても同一であると考えることができる（例えば、他のマイクペア６１１の音源方向ベクトルｑｋは音源方向ベクトル６０２と同一方向）。そのため、音源方向ベクトルｑｋについての条件である式［１５．５］はマイクペアの個数分だけ得られる。
【０４８９】
それら複数の条件式を最もよく近似するものとするｑｋを、音源方向ベクトルの推定値ｈａｔ（ｑｋ）とすればよい。
【０４９０】
そのようなｈａｔ（ｑｋ）を求める方法として、特許第３８８１３６７号では、式［１５．６］に相当する操作が提案されている。ただし、この式の右辺の第１項についている上付きの"＋"は、疑似逆行列（厳密にはＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列）を表わす。また、右辺の第１項は全てのマイクペアについて、
ｐｉｍＴ／｜ｐｉｍ｜
を縦に列挙して生成された行列（行数はマイクペア数、列数は３）、第２項は同じく全てのマイクペアについて、
Ｄｉｍｋ／｜ｐｉｍ｜
を縦に列挙して生成された列ベクトルである。
（ただし、特許第３８８１３６７号における音源方向ベクトル算出の目的はパーミュテーション（ｐｅｒｍｕｔａｔｉｏｎ）の解消にあるため、式［１５．６］に該当する操作は周波数ビンごとに行なわれている。また、経路差Ｄｉｍｋに相当する量は、直線の傾きからではなく、周波数ビンごとに分離行列の逆行列から算出されている。）
【０４９１】
しかし、Ｍｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列を用いて導出された式［１５．６］には、ｑｋが単位ベクトルであるという制約（式［１５．２］）が反映されていない。従って、式［１５．２］の制約も考慮して推定された最適なｈａｔ（ｑｋ）とは異なる方向が算出される場合がある。
【０４９２】
そこで本発明では、不動点法に基づく反復アルゴリズムを用いて、最適な音源方向ベクトルを推定する。アルゴリズムの説明の前に、原理について説明する。
【０４９３】
マイクペアごとの音源方向ベクトルｑｋの条件式である式［１５．５］に加えて、ｑｋが単位ベクトルという条件（式［１５．２］）も一つの式で表わすと、式［１５．９］のように書くことができる。
ただし、Ｐｋ，Ｄｋは、それぞれ式［１５．７］，式［１５．８］で表わされるベクトルと行列である。
【０４９４】
Ｐｋは、各マイクペアのベクトルｐｉｍについて、大きさを１に正規化してから横方向に並べて生成した行列である。必ずしも全部のマイクペアベクトルを列挙する必要はなく、直線当てはめ誤差ｅｒｒｉｍｋの大きいマイクペアについては省略しても良い。また、マイクペアの方向と音源方向との関係によっては、直線の傾きがそれほど高精度には算出できない場合もある。そのような位置関係にあるマイクペアについても、除外してよい（詳細は後述する）。
【０４９５】
Ｄｋは、同じく各マイクペアについて、経路差Ｄｉｍｋをマイク間隔で正規化してから横方向に並べて生成したベクトルである。Ｐｋと同様にＤｋについても、当てはめ誤差ｅｒｒｉｍｋが大きなマイクペアについては省略しても良い。
【０４９６】
ＰｋおよびＤｋに採用されたマイクペアの個数をＮｋとすると、式［１５．９］の左辺の第１項は（Ｎｋ＋１）×３の行列（「縦の要素数×横の要素数」の形式）、右辺は縦方向の要素数がＮｋ＋１のベクトルである。
【０４９７】
不動点法とは、式［１５．９］を"ｑｋ＝"の形に変形した式を用意し（式［１５．１０］）、この式に基づいてｑｋの更新を繰り返すというものである。更新を繰り返すことにより、最終的にｑｋは式［１５．９］を最も良く近似する値に収束する。ただし本発明では、式［１５．２］の条件を一層強く反映させるため、式［１５．１０］と式［１５．１１］とを交互に繰り返す。
【０４９８】
図４７に示すフローチャートは、音源方向ベクトル推定アルゴリズムの概要である。
ステップＳ５０１〜Ｓ５０７はチャンネルについてのループであり、図４２のフローにおける「チャンネルループ」と同様である。
【０４９９】
ステップＳ５０２において、式［１５．７］のＰｋおよび式［１５．８］のＤｋを生成する。すなわち、
各マイクペアのベクトルｐｉｍについて、大きさを１に正規化してから横方向に並べて生成した行列［Ｐｋ］
各マイクペアについて、経路差Ｄｉｍｋをマイク間隔で正規化してから横方向に並べて生成したベクトル［Ｄｋ］
これらを生成する。詳細は後述する。
ＰｋおよびＤｋが生成された後では、Ｎｋという変数に「採用されたマイクペアの個数」が代入されている。
【０５００】
図４２のフローを参照して説明した処理と同様に、直線当てはめに失敗している場合はどのマイクペアもＰｋおよびＤｋには反映されない。また、音源方向ベクトルｑｋは３次元のベクトルであるため、音源方向ベクトルｑｋを推定するためにはＮｋが３以上である必要がある。そこで、ステップＳ５０３において、「採用されたマイクペアの個数」Ｎｋの値によって処理を３通りに分岐させる。
【０５０１】
「採用されたマイクペアの個数」Ｎｋ＝０の場合は、ステップＳ５０４に進む。この場合は方向推定に失敗しているか、またはそのチャンネルがどの音源にも対応していないと考えられるため、「無効」を表わす値を音源方向ベクトルｑｋに代入する。
【０５０２】
Ｎｋ＞＝３の場合は、ステップＳ５０６に進み、その中で不動点法に基づいて音源方向ベクトルｑｋを推定する。詳細は後述する。
【０５０３】
０＜Ｎｋ＜３の場合は、３次元ベクトルである音源方向ベクトルｑｋを推定することはできないが、音源方向の角度であるθｋを推定することはできる。そこでその場合はステップＳ５０５に進み、その中において、先に説明した図４２の処理フローに従った方法でθｋを算出する。
【０５０４】
次に、ステップＳ５０２の式［１５．７］のＰｋおよび式［１５．８］のＤｋの生成処理、すなわち、
各マイクペアのベクトルｐｉｍについて、大きさを１に正規化してから横方向に並べて生成した行列［Ｐｋ］
各マイクペアについて、経路差Ｄｉｍｋをマイク間隔で正規化してから横方向に並べて生成したベクトル［Ｄｋ］
これらの生成処理の詳細について、図４８に示すフローを参照して説明する。
【０５０５】
ステップＳ５１１において、何個のマイクペアがＰｋおよびＤｋに反映されたかを表わす変数Ｎｋを０に設定する初期値設定を行う。
【０５０６】
ステップＳ５１２〜Ｓ５１７はマイクペアについてのループであり、これは図４２のフローにおける「マイクペアループ」と同様である。ステップＳ５１３は、図４２のフローにおけるステップＳ４０４の処理と同様であり、チャンネルｋ・マイクペア（ｉ，ｍ）に対応した当てはめ誤差ｅｒｒｉｍｋが、所定の閾値を上回っているかどうかで処理を分岐させる。
【０５０７】
当てはめ誤差ｅｒｒｉｍｋが閾値を上回っている場合は、以降のステップＳ５１４〜Ｓ５１６をスキップする。当てはめ誤差ｅｒｒｉｍｋが閾値を下回っている場合はＰｋおよびＤｋにそれぞれｐｉｍおよびＤｉｍｋを追加（ステップＳ５１４,Ｓ５１５）した後、ステップＳ５１６において、Ｎｋを＋１する。
【０５０８】
ステップＳ５１４では、その時点でのＰｋに対して、ｐｉｍを新しい列としてＰｋの右側に追加するという処理である。ステップＳ５１５もほぼ同様に、その時点のＤｋに対して、Ｄｉｍｋを新しい要素としてＤｋの右側に追加する処理である。
【０５０９】
ステップＳ５１３〜Ｓ５１６をマイクペアの数だけループさせることで、最終的にＰｋとＤｋとが生成され、さらにＮｋにＰｋおよびＤｋの横の要素数が代入された状態になる。
【０５１０】
次に、図４７に示すフローチャートにおけるステップＳ５０６の「ｑｋの推定」について、図４９に示すフローを参照して説明する。
【０５１１】
ステップＳ５２１において、ｑｋに音源方向ベクトルの初期値を代入する。初期値としては、前回推定したｑｋの値を用いる方法がある。あるいは、音源の存在する方向がある程度限定できる場合は、それを初期値に反映することができる。例えば、図５０に示すように、音源の存在しそうな範囲７０１が、全マイク７０２の重心７０３から見てｙ軸方向におよそ距離Ｌ、すなわち、図に示すマイク重心音源間距離Ｌ７０４だけ離れており、さらに、マイク重心７０３から見てｚ軸方向（ｘｙ平面に垂直の方向）におよそＨだけ高い位置にあることが予め分かっている場合は、ｑｋとして［０，Ｌ，Ｈ］Ｔを代入した後、式［１５．１１］によって大きさを１に正規化すればよい。
【０５１２】
また、カメラ画像や超音波センサー赤外線センサーといった、音以外のモーダルやセンサーなどから、音源になりそうな物体の方向が予想できる場合は、それらの手段から求まった方向をｑｋの初期値としても良い。
【０５１３】
次に、ステップＳ５２３の「ｑｋの更新」処理を、ｑｋが収束するまで繰り返し、ｑｋが収束したと判定されたら（ステップＳ５２２で「はい」）、終了する。収束したかどうかを判定するためには、前回の更新値をｑ'ｋという変数に保存しておき、今回の更新結果との差分｜ｑｋ−ｑ'ｋ｜がある値よりも小さくなったかどうか（例えば０．００００１より小さいか）で判定すればよい。または単純に、５回や１０回といった一定回数だけループしても良い。
【０５１４】
ステップＳ５２３の「ｑｋの更新」とは、式［１５．１０］を実行してから式［１５．１１］を実行することである。式［１５．１０］と式［１５．１１］を交互に繰り返すことで、ｑｋは最終的に式［１５．８］を最も良く近似する値に収束するのである。
【０５１５】
なお、図４８のＳ５１３では、当てはめ誤差ｅｒｒｉｍｋの値に従ってマイクペアを選別していた。この条件に加えて、Ｋｉｍｋの値が所定の閾値を上回っていたらそのマイクペアの使用を却下するという条件分岐を加えてもよい。以降では、その理由と方法について説明する。
【０５１６】
マイクペアにとって、傾きＫｉｍｋが高い分解能で推定できる方向と低い分解能でしか推定できない方向とが存在する。図１で説明すると、角度θが０に近い場合は音源方向が高い分解能で算出可能だが、θが±πラジアン（±９０°）に近い場合は低い分解能でしか算出できない。そのため、図４６のように、音源に対して様々な位置関係のマイクペアが存在する場合は、全部のマイクペアを用いて音源方向を算出するよりも、高い分解能で音源方向を推定している思われるマイクペアを選別し、それらの間で音源方向を算出した方が、方向推定の精度が高くなる。
【０５１７】
マイクペア選別の具体的な方法は、以下の通りである。事前に、音源がどの範囲の角度ある場合にそのマイクペアを方向推定に採用するかを決めておき、上限をθｍａｘとおく。例えば、図１において音源方向の角度θが−２／３π（−６０°）から＋２／３π（＋６０°）までの間にある場合にそのマイクペアを採用するのであれば、θｍａｘ＝２／３πとする。次に角度の最大値に対応する直線の傾きを以下に示す式［１６．１］で計算する。
【０５１８】
【数２４】
2010020294.tif 000027
【０５１９】
この式は、式［８．４］の推定音源角度ｈａｔ（θｉｍｋ）をθｍａｘに置き換えたものであり、この式で算出される傾き最大値をＫｍａｘ＿ｉｍとおく。Ｋｍａｘ＿ｉｍはマイクペアごとに決まる値であり、事前に計算しておくことができる。
【０５２０】
あるマイクペアにとって音源方向が所定の範囲を超えているかどうかは、Ｋｉｍｋの絶対値がＫｍａｘ＿ｉｍを超えているかによって判定することができる（式［１６．２］）。つまり、図４３のＳ５１３において、「ｅｒｒｉｍｋ＞閾値」という判定をする代わりに、「ｅｒｒｉｍｋ＞閾値、または｜Ｋｉｍｋ｜＞Ｋｍａｘ＿ｉｍ」という判定を行ない、条件が満たされる場合は右へ分岐（そのマイクペアの使用を却下）、条件が満たされない場合は下へ分岐すればよい。
【０５２１】
以上、特定の実施例を参照しながら、本発明について詳解してきた。しかしながら、本発明の要旨を逸脱しない範囲で当業者が実施例の修正や代用を成し得ることは自明である。すなわち、例示という形態で本発明を開示してきたのであり、限定的に解釈されるべきではない。本発明の要旨を判断するためには、特許請求の範囲の欄を参酌すべきである。
【０５２２】
また、明細書中において説明した一連の処理はハードウェア、またはソフトウェア、あるいは両者の複合構成によって実行することが可能である。ソフトウェアによる処理を実行する場合は、処理シーケンスを記録したプログラムを、専用のハードウェアに組み込まれたコンピュータ内のメモリにインストールして実行させるか、あるいは、各種処理が実行可能な汎用コンピュータにプログラムをインストールして実行させることが可能である。例えば、プログラムは記録媒体に予め記録しておくことができる。記録媒体からコンピュータにインストールする他、ＬＡＮ（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ）、インターネットといったネットワークを介してプログラムを受信し、内蔵するハードディスク等の記録媒体にインストールすることができる。
【０５２３】
なお、明細書に記載された各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。また、本明細書においてシステムとは、複数の装置の論理的集合構成であり、各構成の装置が同一筐体内にあるものには限らない。
【産業上の利用可能性】
【０５２４】
以上、説明したように、本発明の一実施例の構成によれば、複数音源からの出力を混合した混合信号からなる観測信号に対して、独立成分分析（ＩＣＡ：Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ）を適用した学習処理により、混合信号を分離する分離行列を求め、観測信号に対して分離行列を適用して分離信号を生成するとともに、分離信号の音源方向を算出する音源方向推定処理を実行する。音源方向推定処理においては、時間周波数領域において対応する時間区間の観測信号と分離信号との間で相互共分散行列を計算し、相互共分散行列の要素同士の位相差を算出し、算出した位相差を適用して各分離信号に対応する音源方向を算出する。具体的には、相互共分散行列の要素同士の位相差情報を利用して位相角と周波数との比例定数（直線の傾き）を算出し、比例定数を利用して各分離信号に対応する音源方向を算出する。この処理によって、高精度な信号分離および音源方向推定処理が実現される。
【符号の説明】
【０５２５】
１１  音源
１２，１３  マイク
６１  学習データブロック
６２  現在時刻の観測信号
６３  分離行列
６４  現在時刻の分離結果
６５  分離結果スペクトログラム区間
６６  現在時刻を含む観測信号ブロック
７１〜７３  適用分離行列規定区間
７４  観測信号蓄積区間
７５  学習開始タイミング
７７，７８  学習区間
８１  観測区間
８５  収束度合いデータ
８８  不整合度合いデータ
９１  学習データブロック区間
９２  観測信号データ区間
９３  分離行列
９４  分離結果データ区間
１０３  アンラップ操作処理領域
１０４  当てはめ直線
１１０  学習処理部
１１１  スレッド演算部
１１１−１〜Ｎ  学習スレッド
１１２  スレッド制御部
１１３  分離行列保持部
１２１  マイクロホン
１２２  ＡＤ変換部
１２３  分離処理部
１２４  後段処理部
１２５  音源方向推定部
１２６  フーリエ変換部
１２７  分離部
１２８  フーリエ逆変換部
１３２  演算部
１３３  観測信号バッファ
１３４  分離結果バッファ
１３５  相互共分散行列バッファ
１３６  位相差バッファ
１３７  直線傾きバッファ
１３８  当てはめ誤差バッファ
１３９  音源方向バッファ
１４２  演算部
１４３  観測信号バッファ
１４４  共分散行列バッファ
１４５  相互共分散行列バッファ
１４６  位相差バッファ
１４７  直線傾きバッファ
１４８  当てはめ誤差バッファ
１４９  音源方向バッファ
１５１  現フレーム番号保持カウンタ
１５２  学習初期値保持部
１５３  蓄積開始予定タイミング指定情報保持部
１５４  観測信号蓄積タイミング情報保持部
１５５  ポインタ保持部
１６１  観測信号バッファ
１６２  分離結果バッファ
１６３  学習演算部
１６４  分離行列保持部
１６５  状態格納部
１６６  カウンタ
１６７  観測信号の開始・終了タイミング保持部
１６８  学習終了フラグ
１６９  前処理用データ保持部
１７０  打ち切りフラグ
１７１，１７４  蓄積中状態
１７２  学習中状態
１７３  待機中状態
１８１  初期状態
１８２  待機中状態
１８３  蓄積中状態
１８４  学習中状態
１９１〜１９３  フレーム
４０１  直線
４０５  直線
４１１〜４１７  直線
４２１〜４２３  マイク
４２４  推定音源
４２７  マイクペアの回転角
５０１  位相角
５０３  位相差
６０１  音源
６０２  音源方向ベクトル
６０３，６０４  マイク
６０５  マイクペアベクトル
６０６  経路差
６１１  他のマイクペア
７０１  音源推定範囲
７０２  マイク
７０３  全マイクの重心
７０４  マイク重心音源間距離
【図１】
2010020294.tif 000028
【図２】
2010020294.tif 000029
【図４】
2010020294.tif 000030
【図８】
2010020294.tif 000031
【図１０】
2010020294.tif 000032
【図１１】
2010020294.tif 000033
【図１２】
2010020294.tif 000034
【図１３】
2010020294.tif 000035
【図１４】
2010020294.tif 000036
【図１５】
2010020294.tif 000037
【図１６】
2010020294.tif 000038
【図１８】
2010020294.tif 000039
【図１９】
2010020294.tif 000040
【図２１】
2010020294.tif 000041
【図２２】
2010020294.tif 000042
【図２３】
2010020294.tif 000043
【図２４】
2010020294.tif 000044
【図２５】
2010020294.tif 000045
【図２６】
2010020294.tif 000046
【図２７】
2010020294.tif 000047
【図２８】
2010020294.tif 000048
【図２９】
2010020294.tif 000049
【図３０】
2010020294.tif 000050
【図３１】
2010020294.tif 000051
【図３２】
2010020294.tif 000052
【図３３】
2010020294.tif 000053
【図３４】
2010020294.tif 000054
【図３５】
2010020294.tif 000055
【図３６】
2010020294.tif 000056
【図３７】
2010020294.tif 000057
【図３８】
2010020294.tif 000058
【図３９】
2010020294.tif 000059
【図４０】
2010020294.tif 000060
【図４１】
2010020294.tif 000061
【図４２】
2010020294.tif 000062
【図４３】
2010020294.tif 000063
【図４４】
2010020294.tif 000064
【図４５】
2010020294.tif 000065
【図４６】
2010020294.tif 000066
【図４７】
2010020294.tif 000067
【図４８】
2010020294.tif 000068
【図４９】
2010020294.tif 000069
【図５０】
2010020294.tif 000070
【図３】
2010020294.tif 000071
【図５】
2010020294.tif 000072
【図６】
2010020294.tif 000073
【図７】
2010020294.tif 000074
【図９】
2010020294.tif 000075
【図１７】
2010020294.tif 000076
【図２０】
2010020294.tif 000077