信号源数の推定方法、推定装置、推定プログラム及び記録媒体
【課題】 実環境において信号源の数を正しく推定する。
【解決手段】 周波数領域変換部２０が、観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換し、信号分離部３１が、この時系列データXj(f,τ)から分離信号Yi(f,τ)を生成してメモリ１０に格納する。次に、パワー算出部３２が、各分離信号Yi(f,τ)のパワー値を求めてメモリに格納し、エンベロープ相関算出部３３が、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出してメモリに格納する。そして、判定部３４が、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、残響レベル及びエンベロープ相関値のそれぞれのしきい値を示すパラメータthnoise、threv及びthcorとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する。
【解決手段】 周波数領域変換部２０が、観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換し、信号分離部３１が、この時系列データXj(f,τ)から分離信号Yi(f,τ)を生成してメモリ１０に格納する。次に、パワー算出部３２が、各分離信号Yi(f,τ)のパワー値を求めてメモリに格納し、エンベロープ相関算出部３３が、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出してメモリに格納する。そして、判定部３４が、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、残響レベル及びエンベロープ相関値のそれぞれのしきい値を示すパラメータthnoise、threv及びthcorとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、混合された未知数の信号を複数のセンサにより観測した観測値を用いて信号の数を推定する技術に関し、特に、実環境において信号源の数を正しく推定する技術に関する。
【背景技術】
【０００２】
複数の音源が混合した観測信号を短時間フーリエ変換し、各周波数ビンで空間相関行列の固有値を調べることにより、信号源の数を推定する方法が提案されている（例えば、非特許文献１参照。）。
［問題の定式化］
まず、この方法で取り扱う問題の定式化を行う。すべての信号はあるサンプリング周波数でサンプリングされ、離散的に表現されるものとする。Ｎ個の信号が混合されてＭ個のセンサで観測されたとする。以下では、信号の発生源からセンサまでに距離があり、信号が減衰・遅延し、かつ複数の経路を経てセンサに到達する状況を扱う。このような状況での混合は、信号源kからセンサjヘのインパルス応答hjk(l)による畳み込み混合
【数２】
となる。ここでPはインパルス応答hjk(l)の持続時間を、nj(t)はセンサでのノイズを表す。具体的な例としては、音信号が室内で混合される場合、音源からマイクまでの距離により音が減衰・遅延し、また壁などの反射により残響が発生し、さらにマイクに背景ノイズが付加される。
【０００３】
［固有値に基づく方法］
次に、非特許文献１で提案された信号源数の推定方法を、順を追って説明する。なお、センサの数は信号源の数と同等かそれ以上、すなわちN≦Mを仮定する。
まず、センサｊでの観測信号xj(t)にＬ点の短時間離散フーリエ変換を適用して周波数毎の時間系列
【数３】
を求める。ここでfは周波数であり、f=0, (1/L)fs,...,{(L-1)/L}fsと離散化されている（fsはサンプリング周波数)。g(l)は窓関数であり、ハニング窓g(l)=(1/2)(1+cos(2πl/L))などのg(0)にパワーの中心を持つ窓関数を用いることで、Xj(f,τ)は時刻t=τを中心とする観測信号xj(t)の周波数特性を表現する。Xj(f,τ)はLサンプルにわたる情報を含んでいるため、すべての時間τに対してXj(f,τ)を求める必要はなく、適当な間隔（例えばL/2やL/4）の時間τ毎にXj(f,τ)を求める。
【０００４】
畳み込み混合された信号には、周波数領域での操作が有効である。式(1)で示される時間領域での畳み込み混合が、周波数領域では
【数４】
と各周波数での単純混合に近似表現できるからである。ここで、Hjk(f)は信号源kからセンサｊまでの周波数応答、Sk(f,τ)やNj(f,τ)は式(2)と同様の式に従って源信号sk(t)やノイズnj(t)に短時間離散フーリエ変換を施したものである。
【０００５】
次に、X(f,τ)=[X1(f,τ),..., XM(f,τ)]Tに対して相関行列R(f)=〈X(f,τ) X(f,τ)H〉τを計算し、これをR(f)=V(f)・Λ(f)・V(f)Hのように固有値分解する。なお、V(f)=[v1(f),v2(f),...,vM(f)]であり、Λ(f)はλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列である。また、・Hは行列の共役転置を求める操作、〈・〉τは時間τに関する平均、vj(f)は固有ベクトル（M次元の縦ベクトル）、λj(f)はこれに対応する固有値であり、λ1(f)≧λ2(f)≧...≧λM(f)の順にソートされている。また、各固有値λj(f)は、[Y1(f,τ),...,YM(f,τ)]T←V(f)H・[X1(f,τ),...,XM(f,τ)]Tとしたときのｊ番目の信号Yj(f,τ)のパワー値を示す。
【０００６】
そして、分解された固有値のうち支配的な値を持つ固有値の個数Nを信号源の数と推定し、残りのM-N個の固有値の大きさをノイズのパワー値σn(f)2と推定する（λN+1(f)=…=λM(f)=σn(f)2）。
【非特許文献１】山本潔，W. F.G. van Rooijen， E. Y. Ling， 浅野太， 山田武志，北脇信彦，「ＳＶＭを用いた音源数推定法の音源分離システムヘの応用」，日本音響学会２００２年秋季研究発表会，２−５−１０，ｐｐ．５３７−５３８，２００２年９月」
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかし、従来技術の固有値に基づく方法では、実際の信号源の数を正しく推定できない場合があるという問題点がある。
例えば、上述の固有値に基づく方法を現実的な状況で用いる場合、以下に挙げる２つの問題を考慮しなければならない。
１つ目の問題は残響の影響である。一般に、残響の長さは短時間離散フーリエ変換のフレーム長Ｌよりも長いため、ある信号のある時刻の成分が複数のフレームに影響する。その結果、支配的な固有値の数が実際の信号の数よりも多く推定されることがある。
【０００８】
図１１（ａ）は、図１０に示す条件で１つの音源だけを鳴らした場合の各周波数における固有値の正規化パワー値である。この図に示すように、残響の影響により２番目に大きな値をとる固有値のパワー値が−２０ｄＢ程度になっている。上述の固有値に基づく方法の場合、所定のしきい値よりも値が大きな固有値の個数を信号源の数と判断することになるが、このしきい値が−２０ｄＢより小さかった場合、上述の２番目に値が大きな固有値も「支配的な固有値」の一つにカウントされ、音源の数が２個と推定されてしまう。すなわち、残響の影響から、このしきい値をある程度大きな値としなければ正確な音源数の推定はできない。
２つ目の問題は、各信号のパワーが固有値に適切に現れていない場合があるということである。この問題は特に位相差が小さくなる低周波数で顕著になる。
【０００９】
図１１（ｂ）は、図１０に示す条件で３音源すべてを鳴らした場合の各周波数における固有値の正規化パワー値である。この例の場合、音源数は３であるから３つの支配的な固有値が存在するはずである。しかし、この図に示すように、各音源のパワー値は同等に設定したにもかかわらず、固有値のパワー値は、２番目、３番目となるにつれ次第に小さくなっていく。この傾向は、低周波数になるほど顕著となる。そのため、この状況において多くの周波数で３音源が存在すると推定されるためには、上述のしきい値を−３０ｄＢ程度以下に設定しなければならない。しかし、しきい値を小さく設定すると、今度は残響に対応する固有値も「支配的な固有値」にカウントされ、例えば、図１１（ａ）の１音源の場合に２音源以上と推定されてしまう。
【００１０】
以上説明してきたように、従来技術の固有値に基づく方法では、残響の影響の問題と、各信号のパワーが固有値に適切に現れていない問題とにより、実際の信号源の数を正しく推定できないことがある。
本発明はこのような点に鑑みてなされたものであり、実環境でも信号源の数を正しく推定できる技術を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明では上記課題を解決するために、まず、Ｍ個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換し、この時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する。次に、上記の各分離信号Yi(f,τ)のパワー値を求めて記憶部に格納し、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープの相関値を算出して記憶部に格納する。なお、分離信号Yi(f,τ)のエンベロープとは、分離信号の絶対値の包絡線｜Yi(f,τ)｜を意味する。そして、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、記憶部に格納されている複数のパラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する。
【００１２】
ここで、源信号と残響信号とは相関を持つため、これらの間のエンベロープ相関値は高い。また、残響信号は対応する源信号よりもパワーが小さい。つまり、エンベロープ相関値が高く、パワー値が比較的小さいのが残響信号である。本発明ではこの特徴に着目し、分離信号のエンベロープ相関値やパワー値を各しきい値を示す複数のパラメータと比較して、その分離信号が源信号であるか否かを判断する。これにより、固有値のパワー値のみを指標として源信号を判別していた場合に比べ、実環境を考慮した信号源数の推定が可能となる。
【発明の効果】
【００１３】
以上のように、本発明では、分離信号のエンベロープ相関値とパワー値とを算出し、それらと複数のパラメータとを比較して源信号を判別することとしたため、実環境において信号源の数を正しく推定することが可能となる。
【発明を実施するための最良の形態】
【００１４】
以下、この発明の実施の形態を図面を参照して説明する。
〔第１の実施の形態〕
まず、本発明における第１の実施の形態について説明する。
＜全体の構成＞
図１は本形態における推定装置１の全体を示すブロック図である。
推定装置１は、例えば、ＣＰＵ（central processing unit）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、ハードディスク等がバスで接続されたノイマン型コンピュータに所定のプログラム（推定プログラム）を実行させることにより構築されるものである。
【００１５】
図１に例示するように、本形態の推定装置１は、メモリ１０、周波数領域変換部２０、信号源数推定部３０、結果統合部４０及び制御部５０を有している。ここで、信号源数推定部３０は、信号分離部３１、パワー算出部３２、エンベロープ相関算出部３３及び判定部３４を有し、メモリ１０は、観測信号領域１１、周波数毎の時系列データ領域１２、分離信号領域１３、パワー値領域１４、エンベロープ領域１５、パラメータ領域１６及び信号源数領域１７を有している。また、制御部５０はレジスタ５１を有し、推定装置１全体を制御する。また、この図における破線の矢印は理論上の情報の流れを示し、実線の矢印は現実のデータの流れ（同時に電気的或いは情報的な接続関係も）を示している。ただし、制御部５０における入出力データの表記は省略してある。
【００１６】
＜処理の概要＞
本形態では、源信号が混合された混合信号をＭ個のセンサで観測した観測信号x1(t),...,xM(t)から源信号の数を推定する。
本形態では、まず前処理として複数のパラメータ（ノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcor）を特定するデータをメモリ１０に格納しておく。入力された時間領域の観測信号xj(t)(j={1,...,M})は、それぞれ、周波数領域変換部２０で周波数毎の時系列データXj(f,τ)に変換され、信号源数推定部３０の信号分離部３１に送られる。信号分離部３１は、周波数ｆ毎にこの時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成する。そして、パワー算出部３２が、各分離信号Yi(f,τ)のパワー値を算出し、エンベロープ相関算出部３３が、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出する。これらが算出されると、判定部３４は、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値とメモリ内の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断して、周波数ｆに対する信号源の数EN(f)を推定する。
【００１７】
その後、音声などの広帯域信号に対しては、最後に結果統合部４０において周波数毎の推定値が統合され、全体としての信号源数の推定値enを得る。一方、通信分野などで用いられる狭帯域信号に対しては、周波数毎の推定値を統合する必要はなく、着目する周波数fでの推定値EN(f)を得れば良い。
＜本形態の詳細＞
図２（ａ）は図１に例示した信号分離部３１の機能構成を、図２（ｂ）はパワー算出部３２の機能構成を、図３（ａ）はエンベロープ相関算出部３３の機能構成を、図３（ｂ）は判定部の機能構成を、それぞれ例示したブロック図である。また、図４及び図５は、本形態における信号源数の推定方法を説明するためのフローチャートである。
【００１８】
以下、図１〜図５を用い、本形態における構成・処理の詳細について説明する。
［前処理］
まず、前処理としてノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorを特定するデータを、メモリ１０（「記憶部」に相当）のパラメータ領域１６に格納する。なおパラメータとしては、例えば、thnoise=0.01、threv =0.2、thcor=0.5を例示できる。ただし、実際の測定時において、鳴っている音源数が分かるサンプルがあれば、その観測データをもとに各パラメータを調整していってもよい。具体的には、例えば、第１パラメータthnoiseを、ノイズ信号の正規化パワー値よりも大きく残響信号の正規化パワー値よりも小さくなるように調整し、第２パラメータthrevを、源信号の正規化パワー値よりも小さく残響信号の正規化パワー値よりも大きくなるように調整し、第３パラメータthcorを、源信号と残響信号とのエンベロープ相関値より小さくなるように調整する。なお、正規化パワー値やエンベロープ相関値の意味については後述する。
また、信号源数の推定対象となる時間領域の観測信号xj(t)(j={1,...,M})をメモリ１０の観測信号領域１１に書き込む。なお、この観測信号xj(t)はＭ個のセンサ（マイクロホン等）での観測信号であり、下付添字のjは、その観測信号xj(t)がj番目のセンサで観測されたことを示す。
【００１９】
［周波数領域への変換］
まず、制御部５０（図１）が変数jに１を代入し、それをレジスタ５１に格納する（ステップＳ１）。次に、周波数領域変換部２０が、メモリ１０の観測信号領域１１にアクセスし、観測信号xj(t)を読み込む（ステップＳ２）。観測信号xj(t)を読み込んだ周波数領域変換部２０は、それを周波数毎の時間系列データXj(f,τ)に変換してメモリ１０の周波数毎の時系列データ領域１２に格納する（ステップＳ３）。なお、この例では、サンプリング周波数fs、Ｌ点の短時間離散フーリエ変換を利用してこの変換を行う。
次に制御部５０は、レジスタ５１に格納された変数jがＭか否かを判断する（ステップＳ４）。ここでj＝Ｍでないと判断された場合、制御部５０がjに１を加算した値を新たなjとし（ステップＳ５）、それをレジスタ５１に格納してステップＳ２の処理に戻る。一方、j＝Ｍであると判断された場合、以下の信号源推定処理に移る。
【００２０】
［信号源推定処理］
まず、制御部５０（図１）が、変数ｆに０を代入してレジスタ５１に格納する（ステップＳ６）。
独立成分分析（ＩＣＡ：Independent Component Analysis）処理：
次に、独立成分分析（ＩＣＡ）部３１ａ（図２（ａ））が、メモリ１０の周波数毎の時系列データ領域１２から時系列データXj(f,τ)を抽出し、独立成分分析（ＩＣＡ）を用い、X(f,τ)=[X1(f,τ),...,XM(f,τ)]Tから、Ｍ×Ｍ行列の分離行列Ｗ(f)とＩＣＡ分離信号Z(f,τ)=[Z1(f,τ),...,ZM(f,τ)]Tとを生成してレジスタ３１ｂ（「記憶部」に相当）に格納する（ステップＳ７）。ここでＩＣＡによる信号分離は、ＩＣＡ分離信号Z(f,τ)の各要素が互いに独立になるようにZ(f,τ)=Ｗ(f)・X(f,τ)となるＷ(f) を算出する手法である。また、ＩＣＡのアルゴリズムは、A. Hyvarinen and J. Karhunen and E. Oja, "Independent Component Analysis," John Wiley & Sons, 2001, ISBN 0-471-40540,などに様々なものが示されている。
【００２１】
なおＩＣＡの解にはスケーリングの任意性がある。Z(f,τ)のある要素にあるスカラ値を掛けても、要素間の独立性は変化しないからである。従って、この段階では、センサで観測された源信号のパワーがＩＣＡ分離信号に正しく反映されていない可能性が高い。また、源信号の数Ｎがセンサ数Ｍより少なければ、ＩＣＡ分離信号Z(f,τ)のＮ個の要素は源信号に対応し、残りのM-N個の要素はノイズや残響成分に対応するが、この段階のノイズや残響に対応する要素の大きさは一般に増幅されている。そこで、次にスケーリング部３１ｃ（図２（ａ））において、このスケーリングの任意性の問題を解決する。
【００２２】
スケーリング問題解決処理：
スケーリング部３１ｃでは、スケーリングの任意性の問題を解決するため、以下に示す操作を行う。まず、対角行列生成部３１ｃａが、レジスタ３１ｂから分離行列Ｗ(f)を読み出し、この分離行列Ｗ(f)からスケーリング問題を解決するための対角行列Λ(f)を生成する（ステップＳ８）。この対角行列Λ(f)としては、例えば、
Λ(f)=sqrt(diag[(W(f)・W(f)H)-1]) …(4)
が例示できる。ここで、・-1は逆行列、・Hは共役転置行列、diagは対角成分以外を０にする操作、sqrtは各要素の平方根を計算する操作である。
【００２３】
生成された対角行列Λ(f)は、積演算部３１ｃｂに送られ、積演算部３１ｃｂは、これとレジスタ３１ｂから読み出したＩＣＡ分離信号Z(f,τ)とを用い、[Y1(f,τ),...,YM(f,τ)]T←Λ(f)・[Z1(f,τ),...,ZM(f,τ)]Tの演算によって、スケーリング問題を解決した（パワーを回復した）分離信号Yi(f,τ)(i={1,...,M})を生成し、メモリ１０の分離信号領域１３（図１）に格納する（ステップＳ９）。
ここで、式(4)を含む上記一連の操作により、分離信号Yi(f,τ)は、以下の２つの性質を持つ。第一に分離信号Yi(f,τ)が互いに無相関であれば、
【数５】
が成り立つ。すなわち、分離信号Yi(f,τ)のパワーの総和とセンサでの観測信号Xj(f,τ)のパワーの総和が等しくなる。さらに、分離信号Yi(f,τ)が互いに独立であれば、
【数６】
【００２４】
が成り立つ。なお、Sk(f,τ)(k={1,...,N})は源信号成分を示す。すなわち、ある分離信号Yi(f,τ)のパワーと、それに対応する源信号Sk(f,τ)をすべてのセンサで観測した際のパワーの総和とは等しくなる。分離信号Yi(f,τ)が互いに無相関、さらには互いに独立になることは、独立成分分析の目的であり、多くの場合この条件はほぼ満たされている。従って、上記一連の操作により、各分離信号Yi(f,τ)のパワーは、それに対応する源信号Sk(f,τ)がセンサで観測された際のパワーの総和に近くなる。
なお、式(4)の対角行列Λ(f)の代わりに、対角行列Λ(f)=diag[W(f)-1]を使用してもよく、より一般的にW(f)-1のｉ列ｊ行目の要素をｊ行目の対角成分とする対角行列Λを使用してもよい。この場合、各分離信号Yi(f,τ)のパワーは、対応する源信号Sk(f,τ)をあるセンサｊで観測したパワー、すなわち|Hjk(f)・Sk(f,τ)|2に近似する。
【００２５】
［判定処理］
判定処理では、スケーリング問題を解決した（パワーを回復した）分離信号Yi(f,τ)から、源信号の数を推定する。まず、制御部５０（図１）が変数ｉに１を代入し、レジスタ５１に格納する（ステップＳ１０）。
次に、パワー算出部３２の平均パワー算出部３２ａ（図２（ｂ））が、例えばメモリ１０の分離信号領域１３（図１）から各τに対する分離信号Yi(f,τ)を順次抽出し、そのパワー値｜Yi(f,τ)｜2を順次算出してレジスタ３２ｂに格納する。そして、平均パワー算出部３２ａは、レジスタ３２ｂに格納されたパワー値｜Yi(f,τ)｜2を読み出し、分離信号Yi(f,τ)の時間τに関する平均パワー値
σi2(f)←〈｜Yi(f,τ)｜2〉τ
を算出して、レジスタ３２ｂ（「記憶部」に相当）に格納する（ステップＳ１１）。
【００２６】
次に、エンベロープ相関算出部３３のエンベロープ算出部３３ａ（図３（ａ））が、例えばメモリ１０の分離信号領域１３（図１）から各τに対する分離信号Yi(f,τ)を順次抽出し、その絶対値｜Yi(f,τ)｜を順次算出してレジスタ３３ｂに格納する。次に、エンベロープ算出部３３ａは、レジスタ３２ｂに格納された絶対値｜Yi(f,τ)｜を読み出し、時間τに関する平均が０になるように分離信号Yi(f,τ)の絶対値｜Yi(f,τ)｜を正規化したエンベロープ
ｖｉ(f,τ)←｜Yi(f,τ)｜-〈｜Yi(f,τ)｜〉τ
を算出してレジスタ３３ｂ（「記憶部」に相当）に格納する（ステップＳ１２）。
【００２７】
次に、制御部５０（図１）が、レジスタ５１に格納された変数ｉがＭであるか否かを判断する（ステップ１３）。ここでｉ＝Ｍでなければ、制御部５０がｉに１を加算し、その値を新たなｉとしレジスタ５１に格納し（ステップＳ１４）、ステップＳ１１に戻る。一方、ｉ＝Ｍであれば、制御部５０は、この変数ｉに１を代入してレジスタ５１に格納し（ステップＳ１５）、以下の処理を実行する。
まず、パワー算出部３２のパワー正規化部３２ｃ（図２（ｂ））が、レジスタ３２ｂから平均パワー値σ12(f),...,σM2(f)を抽出し、平均パワー値σi2(f)を正規化した正規化パワー値
【数７】
を算出して、メモリ１０のパワー値領域１４に格納する（ステップＳ１６）。
【００２８】
次に、制御部５０（図１）が、変数ｋに１を代入し、レジスタ５１に格納する（ステップＳ１７）。次に、エンベロープ相関算出部３３の相関算出部３３ｃ（図３（ａ））が、レジスタ３３ｂからエンベロープvi(f,τ)及びvk(f,τ)を抽出する。そして、相関算出部３３ｃは、これらのエンベロープvi(f,τ)及びvk(f,τ)を用い、時間差Δτ（例えばL/2やL/4）による分離信号Yi(f,τ)の分離信号Yk(f,τ)とのエンベロープ相関値
【数８】
を算出し、その演算結果Cori,k(f)をレジスタ３３ｂに格納する（ステップＳ１８）。なお、この例のΔτは、例えばプログラムコードに組み込まれた定数である。
【００２９】
次に、制御部５０（図１）はレジスタ５１に格納された変数ｋがＭであるか否かを判断する（ステップＳ１９）。ここで、ｋ＝Ｍでなかった場合、制御部５０がｋに１を加算し、その値を新たなｋとしてレジスタ５１に格納し、ステップＳ１８に戻る（ステップＳ２０）。一方、ｋ＝Ｍであった場合、エンベロープ相関算出部３３の最大値算出部３３ｄ（図３（ａ））は、レジスタ３３ｂ（図３（ａ））から、エンベロープ相関値Cori,1(f),...,Cori,M(f)を抽出する。そして、最大値算出部３３ｄは、これらを用い、エンベロープ相関値Cori,k(f)のｉごとの最大値maxCori(f)を算出し、メモリ１０のエンベロープ領域１５（図１）に格納する（ステップＳ２１）。
【００３０】
次に、判定部３４の比較部３４ａ（図３（ｂ））が、メモリ１０のパワー値領域１４、エンベロープ領域１５及びパラメータ領域１６から、平均パワー値の正規化値NPi(f)、エンベロープ相関値の最大値maxCori(f)並びに第１パラメータthnoise、第２パラメータthrev及び第３パラメータthcorを読み出す。そして、比較部３４ａは、以下の論理式により、分離信号Yi(f,τ)が、源信号に対応するか、ノイズや残響成分に対応するかを判定する（ステップＳ２２）。
【数９】
【００３１】
すなわち、ここでは３種類のパラメータthnoise、threv、thcorを用いている。そして、平均パワー値の正規化値NPi(f)が第１パラメータthnoise未満であればノイズ成分と判定し、平均パワー値の正規化値NPi(f)が第２パラメータthrev未満であり、さらにエンベロープ相関値の最大値maxCori(f)が第３パラメータthcorを超えれば残響成分と判定する。結局、sigi(f)が0になれば、分離信号Yi(f,τ)がノイズや残響成分に対応する（源信号成分でない）と判定されたことになり、sigi(f)=1になれば、分離信号Yi(f,τ)が源信号に対応すると判定されたことになる。そして、このように生成された判定結果sigiはレジスタ３４ｂ（図３（ｂ））に送られて格納される。なお、上記論理式中の「＜」の少なくとも一部を「≦」としてもよく、「＞」を「≧」としてもよい。
【００３２】
次に、制御部５０はレジスタ５１（図１）に格納されている変数ｉがＭであるか否かを判断する（ステップＳ２３）。ここで、ｉ＝Ｍでなければ、制御部５０がｉに１を加算し、その値を新たなｉとしてレジスタ５１に格納してステップＳ１６に戻る（ステップＳ２４）。一方、ｉ＝Ｍであれば、判定部３４の信号源数算出部３４ｃ（図３（ｂ））が、レジスタ３４ｂから判定結果sig1(f),...,sigM(f)を抽出し、信号源数推定値
EN(f)=Σi sigi(f)
を算出し、それをメモリ１０の信号源数領域１７（図１）に格納する（ステップＳ２５）。
【００３３】
次に、制御部５０は、レジスタ５１に格納された変数ｆが{(L-1)/L}fs（fsはサンプリング周波数）であるか否かを判断する（ステップＳ２６）。ここで、変数ｆが{(L-1)/L}fsでなかった場合、制御部５０が変数ｆにfs/Lを加算し、その値を新たな変数ｆとし、レジスタ５１に格納してステップＳ７の処理に戻る（ステップＳ２７）。一方、変数ｆが{(L-1)/L}fsであった場合、以下の結果統合処理を行う。
［結果統合処理］
まず結果統合部４０が、メモリ１０の信号源数領域１７から各周波数fで推定された信号源数推定値EN(0),...,EN({(L-1)/L}fs)を読み出し、これを元に、全体としての信号源数の推定値enを算出して出力する（ステップＳ２８）。この例では、単純に多数決で全体の推定値enを決定する。信頼できる周波数（例えば、高い周波数）に大きな重みを与えて、重みづけの多数決で全体の推定値enを決定しても良い。
【００３４】
［適用結果］
本形態の信号源数の推定方法を音源数の推定に適用した結果を示す。
図１０に一般的な実験条件を例示する。この実験条件は以下である。
・信号源：７秒間の音声
・残響時間：ＴＲ＝２００ ｍｓ
・背景ノイズパワー：−２１．８ ｄＢ
・サンプリング周波数：ｆｓ＝８０００ Ｈｚ
・部屋の大きさ：４．４５ ｍ×３．５５ ｍ×２．５０ ｍ
・音源数：１〜３個
・音源配置・間隔：４ｃｍの間隔で直線上に配置
・センサの数：３個
・中心音源と各センサとの距離：１．１ ｍ
・中心音源と各センサを結んだ直線と、各センサが配置される直線とがなす角度：４５°，９０°，１２０°
この図１０に示す条件で１〜３個の音源を鳴らし、３個のマイクでの観測信号を用いて鳴っている音源の数を推定した。
【００３５】
図６に従来手法と本形態の手法とによる推定結果の比較を示す。ここで、図６（ａ）は固有値に基づく従来手法による信号源数の推定結果を示しており、図６（ｂ）は本形態の手法よる信号源数の推定結果を示している。また、横軸は真の音源数、縦軸は音源数0,1,2,3としてそれぞれ推定した周波数ビンの数を示す。
この図に示すように、従来手法では、１音源や３音源の場合にも多数決によると２音源と推定してしまっている。このように従来手法が推定を誤る原因は、図１１を用いて説明したように、個々の音源やノイズのパワーが固有値に適切に現れていないことや、残響の影響を考慮されていないことである。一方、本形態の手法によるとすべての場合に正しく推定されている。
【００３６】
次に本形態の手法による推定が正確である理由を示す。
まず、パワーの回復（各分離信号が、各信号のパワーを適切に反映しているか）に関して考察する。図７は、３音源の場合にセンサで観測された真のパワー値（図７（ａ））と、これらの混合音を本形態の手法により分離した分離信号（図７（ｂ））のパワー値との比較を示すものである。なお、図７（ａ）の観測結果は、各音源を１つずつ鳴らして測定し、その結果を正規化したものである。これらの図に示すように、本形態の手法による各分離信号のパワー値は、各音源の観測値の真のパワー値に近似し、音源数を推定できる程度に正しくパワーが回復されていることがわかる。
【００３７】
次に、残響の影響ヘの対処に関して考察する。図８は、１音源の場合の１番目と２番目（i=1,2）の分離信号のパワー値（図８（ａ））とそれらのエンベロープの相関値（図８（ｂ））を示すものである。図８（ａ）のパワー値だけを見ると、２番目の分離信号のパワー値が決して十分には小さくないので、信号源なのか残響を含むノイズなのか判断し難い。しかし、右側に示す１番目と２番目の分離信号のエンベロープの相関値を見ると、その値が十分に大きいため、２番目の分離信号は１番目の分離信号の残響成分を多く含むノイズであることがわかる。すなわち、エンベロープの相関値は−１〜＋１の値をとり、信号間の相関性が低いほど０に近づく。図８（ｂ）の例では、エンベロープの相関値が０．６〜１の間に集中しており、１番目の分離信号と２番目の分離信号の相関性が高いことが分かる。そしてパワー値が弱い２番目の分離信号が１番目の信号の残響成分であることが推定できる。
【００３８】
そして、これらの判断に必要なノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorを適切に設定することにより、ノイズや残響の影響が無視できない実環境において、アクティブな源信号の数を精度良く推定することができる。
〔第２の実施の形態〕
次に、本発明における第２の実施の形態について説明する。
本形態は第１の実施の形態の変形例であり、ＩＣＡを用いた信号分離の代わりに固有値に基づく信号分離を行う形態である。以下では、第１の実施の形態との相違点を中心に説明を行い、第１の実施の形態と共通する事項については説明を省略する。
【００３９】
図９（ａ）は、本形態における信号分離部１３１の構成を例示したブロック図である。
なお、本形態の推定装置と第１の実施の形態の推定装置１との相違点は、信号分離部３１が信号分離部１３１になる点のみである。また、本形態の処理と第１の実施の形態の処理との相違点は、信号分離処理（図４：ステップＳ７〜９）と平均パワー算出処理（図４：ステップＳ１１）のみである。
【００４０】
［信号分離処理］
図９（ｂ）は、本形態の信号分離処理を説明するためのフローチャートである。
まず、信号分離部１３１の相関行列生成部１３１ａ（図９（ａ））が、メモリ１０の周波数毎の時系列データ領域１２（図１）から時系列データXj(f,τ)を順次抽出し、時系列ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tに対する相関行列R(f)=〈X(f,τ)・X(f,τ)H〉τを生成する（ステップＳ３１）。
生成された相関行列R(f)は固有値分解部１３１ｂ（図９（ａ））に送られ、固有値分解部１３１ｂはこの相関行列R(f)を、R(f)=V(f)・Λ(f)・V(f)Hの積に分解する（ステップＳ３２）。なお、V(f)=[v1(f),v2(f),...,vM(f)]とし、Λ(f)をλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列とし、vj(f)を固有ベクトルとし、λj(f)をこれに対応する固有値とする。生成された固有値λ1(f)，λ2(f)，...，λM(f)は対応するτに関連つけてメモリ１０（図１）に格納され（ステップＳ３２）、V(f)は積演算部１３１ｃ（図９（ａ））に送られる。
【００４１】
積演算部１３１ｃは、メモリ１０の周波数毎の時系列データ領域１２（図１）から時系列データXj(f,τ)を抽出し、[Y1(f,τ),...,YM(f,τ)]T=V(f)H・[X1(f,τ),...,XM(f,τ)]Tの演算によって、分離信号Yi(f,τ)(i={1,...,M})を生成してメモリ１０の分離信号領域１３に格納する（ステップＳ３３）。なお、〈｜Yi(f,τ)｜2〉τ=λi(f)が成立する。
【００４２】
［平均パワー算出処理］
本形態では、第１の実施の形態のステップＳ１１において、パワー算出部３２の平均パワー算出部３２ａ（図２（ｂ））が、分離信号Yi(f,τ)からパワー値｜Yi(f,τ)｜2を算出し、平均パワー値σi2(f)←〈｜Yi(f,τ)｜2〉τを算出してレジスタ３２ｂに格納していた代わりに、平均パワー算出部３２ａがメモリ１０（図１）から固有値λi(f)を順次抽出し、分離信号Yi(f,τ)の時間τに関する平均パワー値
σi2(f)←λi(f)
を算出して、レジスタ３２ｂに格納する。
なお、その他の処理については第１の実施の形態と同様である。
【００４３】
以上のような構成の場合、パワーが適切に回復されない問題（各信号のパワーが固有値に適切に現れていない問題）は解決されないが、第１の実施の形態と同様、残響の影響の問題は解決できる。そのため、パワーが適切に回復されない問題の影響が少ない周波数領域では、本形態でも正確な信号源数の推定ができる。
なお、本発明は上述の各実施の形態に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００４４】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００４５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００４６】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【００４７】
本発明の音信号に対する応用例としては、例えば、適応ビームフォーマやブラインド音源分離の前処理において、ある区間でのアクティブな音源数を推定する処理を例示できる。
【図面の簡単な説明】
【００４８】
【図１】第１の実施の形態における推定装置の全体を示すブロック図。
【図２】（ａ）は図１に例示した信号分離部の機能構成を例示したブロック図。（ｂ）はパワー算出部の機能構成を例示したブロック図。
【図３】（ａ）はエンベロープ相関算出部の機能構成を例示したブロック図。（ｂ）は判定部の機能構成を例示したブロック図。
【図４】第１の実施の形態における信号源数の推定方法を説明するためのフローチャート。
【図５】第１の実施の形態における信号源数の推定方法を説明するためのフローチャート。
【図６】（ａ）は固有値に基づく従来手法による信号源数の推定結果を示したグラフ。（ｂ）は本形態の手法よる信号源数の推定結果を示したグラフ。
【図７】（ａ）は３音源の場合にセンサで観測された真のパワー値を示した図。（ｂ）は混合音を第１の実施の形態の手法により分離した分離信号のパワー値を示した図。
【図８】（ａ）は１音源の場合の１番目と２番目の分離信号のパワー値を示した図。（ｂ）は、それらのエンベロープの相関値を示した図。
【図９】（ａ）は、第２の実施の形態における信号分離部の構成を例示したブロック図。（ｂ）は第２の実施の形態における信号分離処理を説明するためのフローチャート。
【図１０】実験条件を示した図。
【図１１】固有値に基づく方法でのパワー推定値を示した図。（ａ）は、１つの音源だけを鳴らした場合の各周波数における固有値の正規化パワー値。（ｂ）は、３音源すべてを鳴らした場合の各周波数における固有値の正規化パワー値。
【符号の説明】
【００４９】
１ 推定装置 ３１ 信号分離部
１０ メモリ ３２ パワー算出部
２０ 周波数領域変換部 ３３ エンベロープ相関算出部
３０ 信号減数推定部 ３４ 判定部
【技術分野】
【０００１】
本発明は、混合された未知数の信号を複数のセンサにより観測した観測値を用いて信号の数を推定する技術に関し、特に、実環境において信号源の数を正しく推定する技術に関する。
【背景技術】
【０００２】
複数の音源が混合した観測信号を短時間フーリエ変換し、各周波数ビンで空間相関行列の固有値を調べることにより、信号源の数を推定する方法が提案されている（例えば、非特許文献１参照。）。
［問題の定式化］
まず、この方法で取り扱う問題の定式化を行う。すべての信号はあるサンプリング周波数でサンプリングされ、離散的に表現されるものとする。Ｎ個の信号が混合されてＭ個のセンサで観測されたとする。以下では、信号の発生源からセンサまでに距離があり、信号が減衰・遅延し、かつ複数の経路を経てセンサに到達する状況を扱う。このような状況での混合は、信号源kからセンサjヘのインパルス応答hjk(l)による畳み込み混合
【数２】
となる。ここでPはインパルス応答hjk(l)の持続時間を、nj(t)はセンサでのノイズを表す。具体的な例としては、音信号が室内で混合される場合、音源からマイクまでの距離により音が減衰・遅延し、また壁などの反射により残響が発生し、さらにマイクに背景ノイズが付加される。
【０００３】
［固有値に基づく方法］
次に、非特許文献１で提案された信号源数の推定方法を、順を追って説明する。なお、センサの数は信号源の数と同等かそれ以上、すなわちN≦Mを仮定する。
まず、センサｊでの観測信号xj(t)にＬ点の短時間離散フーリエ変換を適用して周波数毎の時間系列
【数３】
を求める。ここでfは周波数であり、f=0, (1/L)fs,...,{(L-1)/L}fsと離散化されている（fsはサンプリング周波数)。g(l)は窓関数であり、ハニング窓g(l)=(1/2)(1+cos(2πl/L))などのg(0)にパワーの中心を持つ窓関数を用いることで、Xj(f,τ)は時刻t=τを中心とする観測信号xj(t)の周波数特性を表現する。Xj(f,τ)はLサンプルにわたる情報を含んでいるため、すべての時間τに対してXj(f,τ)を求める必要はなく、適当な間隔（例えばL/2やL/4）の時間τ毎にXj(f,τ)を求める。
【０００４】
畳み込み混合された信号には、周波数領域での操作が有効である。式(1)で示される時間領域での畳み込み混合が、周波数領域では
【数４】
と各周波数での単純混合に近似表現できるからである。ここで、Hjk(f)は信号源kからセンサｊまでの周波数応答、Sk(f,τ)やNj(f,τ)は式(2)と同様の式に従って源信号sk(t)やノイズnj(t)に短時間離散フーリエ変換を施したものである。
【０００５】
次に、X(f,τ)=[X1(f,τ),..., XM(f,τ)]Tに対して相関行列R(f)=〈X(f,τ) X(f,τ)H〉τを計算し、これをR(f)=V(f)・Λ(f)・V(f)Hのように固有値分解する。なお、V(f)=[v1(f),v2(f),...,vM(f)]であり、Λ(f)はλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列である。また、・Hは行列の共役転置を求める操作、〈・〉τは時間τに関する平均、vj(f)は固有ベクトル（M次元の縦ベクトル）、λj(f)はこれに対応する固有値であり、λ1(f)≧λ2(f)≧...≧λM(f)の順にソートされている。また、各固有値λj(f)は、[Y1(f,τ),...,YM(f,τ)]T←V(f)H・[X1(f,τ),...,XM(f,τ)]Tとしたときのｊ番目の信号Yj(f,τ)のパワー値を示す。
【０００６】
そして、分解された固有値のうち支配的な値を持つ固有値の個数Nを信号源の数と推定し、残りのM-N個の固有値の大きさをノイズのパワー値σn(f)2と推定する（λN+1(f)=…=λM(f)=σn(f)2）。
【非特許文献１】山本潔，W. F.G. van Rooijen， E. Y. Ling， 浅野太， 山田武志，北脇信彦，「ＳＶＭを用いた音源数推定法の音源分離システムヘの応用」，日本音響学会２００２年秋季研究発表会，２−５−１０，ｐｐ．５３７−５３８，２００２年９月」
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかし、従来技術の固有値に基づく方法では、実際の信号源の数を正しく推定できない場合があるという問題点がある。
例えば、上述の固有値に基づく方法を現実的な状況で用いる場合、以下に挙げる２つの問題を考慮しなければならない。
１つ目の問題は残響の影響である。一般に、残響の長さは短時間離散フーリエ変換のフレーム長Ｌよりも長いため、ある信号のある時刻の成分が複数のフレームに影響する。その結果、支配的な固有値の数が実際の信号の数よりも多く推定されることがある。
【０００８】
図１１（ａ）は、図１０に示す条件で１つの音源だけを鳴らした場合の各周波数における固有値の正規化パワー値である。この図に示すように、残響の影響により２番目に大きな値をとる固有値のパワー値が−２０ｄＢ程度になっている。上述の固有値に基づく方法の場合、所定のしきい値よりも値が大きな固有値の個数を信号源の数と判断することになるが、このしきい値が−２０ｄＢより小さかった場合、上述の２番目に値が大きな固有値も「支配的な固有値」の一つにカウントされ、音源の数が２個と推定されてしまう。すなわち、残響の影響から、このしきい値をある程度大きな値としなければ正確な音源数の推定はできない。
２つ目の問題は、各信号のパワーが固有値に適切に現れていない場合があるということである。この問題は特に位相差が小さくなる低周波数で顕著になる。
【０００９】
図１１（ｂ）は、図１０に示す条件で３音源すべてを鳴らした場合の各周波数における固有値の正規化パワー値である。この例の場合、音源数は３であるから３つの支配的な固有値が存在するはずである。しかし、この図に示すように、各音源のパワー値は同等に設定したにもかかわらず、固有値のパワー値は、２番目、３番目となるにつれ次第に小さくなっていく。この傾向は、低周波数になるほど顕著となる。そのため、この状況において多くの周波数で３音源が存在すると推定されるためには、上述のしきい値を−３０ｄＢ程度以下に設定しなければならない。しかし、しきい値を小さく設定すると、今度は残響に対応する固有値も「支配的な固有値」にカウントされ、例えば、図１１（ａ）の１音源の場合に２音源以上と推定されてしまう。
【００１０】
以上説明してきたように、従来技術の固有値に基づく方法では、残響の影響の問題と、各信号のパワーが固有値に適切に現れていない問題とにより、実際の信号源の数を正しく推定できないことがある。
本発明はこのような点に鑑みてなされたものであり、実環境でも信号源の数を正しく推定できる技術を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明では上記課題を解決するために、まず、Ｍ個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換し、この時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する。次に、上記の各分離信号Yi(f,τ)のパワー値を求めて記憶部に格納し、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープの相関値を算出して記憶部に格納する。なお、分離信号Yi(f,τ)のエンベロープとは、分離信号の絶対値の包絡線｜Yi(f,τ)｜を意味する。そして、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、記憶部に格納されている複数のパラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する。
【００１２】
ここで、源信号と残響信号とは相関を持つため、これらの間のエンベロープ相関値は高い。また、残響信号は対応する源信号よりもパワーが小さい。つまり、エンベロープ相関値が高く、パワー値が比較的小さいのが残響信号である。本発明ではこの特徴に着目し、分離信号のエンベロープ相関値やパワー値を各しきい値を示す複数のパラメータと比較して、その分離信号が源信号であるか否かを判断する。これにより、固有値のパワー値のみを指標として源信号を判別していた場合に比べ、実環境を考慮した信号源数の推定が可能となる。
【発明の効果】
【００１３】
以上のように、本発明では、分離信号のエンベロープ相関値とパワー値とを算出し、それらと複数のパラメータとを比較して源信号を判別することとしたため、実環境において信号源の数を正しく推定することが可能となる。
【発明を実施するための最良の形態】
【００１４】
以下、この発明の実施の形態を図面を参照して説明する。
〔第１の実施の形態〕
まず、本発明における第１の実施の形態について説明する。
＜全体の構成＞
図１は本形態における推定装置１の全体を示すブロック図である。
推定装置１は、例えば、ＣＰＵ（central processing unit）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、ハードディスク等がバスで接続されたノイマン型コンピュータに所定のプログラム（推定プログラム）を実行させることにより構築されるものである。
【００１５】
図１に例示するように、本形態の推定装置１は、メモリ１０、周波数領域変換部２０、信号源数推定部３０、結果統合部４０及び制御部５０を有している。ここで、信号源数推定部３０は、信号分離部３１、パワー算出部３２、エンベロープ相関算出部３３及び判定部３４を有し、メモリ１０は、観測信号領域１１、周波数毎の時系列データ領域１２、分離信号領域１３、パワー値領域１４、エンベロープ領域１５、パラメータ領域１６及び信号源数領域１７を有している。また、制御部５０はレジスタ５１を有し、推定装置１全体を制御する。また、この図における破線の矢印は理論上の情報の流れを示し、実線の矢印は現実のデータの流れ（同時に電気的或いは情報的な接続関係も）を示している。ただし、制御部５０における入出力データの表記は省略してある。
【００１６】
＜処理の概要＞
本形態では、源信号が混合された混合信号をＭ個のセンサで観測した観測信号x1(t),...,xM(t)から源信号の数を推定する。
本形態では、まず前処理として複数のパラメータ（ノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcor）を特定するデータをメモリ１０に格納しておく。入力された時間領域の観測信号xj(t)(j={1,...,M})は、それぞれ、周波数領域変換部２０で周波数毎の時系列データXj(f,τ)に変換され、信号源数推定部３０の信号分離部３１に送られる。信号分離部３１は、周波数ｆ毎にこの時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成する。そして、パワー算出部３２が、各分離信号Yi(f,τ)のパワー値を算出し、エンベロープ相関算出部３３が、異なる分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出する。これらが算出されると、判定部３４は、各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値とメモリ内の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断して、周波数ｆに対する信号源の数EN(f)を推定する。
【００１７】
その後、音声などの広帯域信号に対しては、最後に結果統合部４０において周波数毎の推定値が統合され、全体としての信号源数の推定値enを得る。一方、通信分野などで用いられる狭帯域信号に対しては、周波数毎の推定値を統合する必要はなく、着目する周波数fでの推定値EN(f)を得れば良い。
＜本形態の詳細＞
図２（ａ）は図１に例示した信号分離部３１の機能構成を、図２（ｂ）はパワー算出部３２の機能構成を、図３（ａ）はエンベロープ相関算出部３３の機能構成を、図３（ｂ）は判定部の機能構成を、それぞれ例示したブロック図である。また、図４及び図５は、本形態における信号源数の推定方法を説明するためのフローチャートである。
【００１８】
以下、図１〜図５を用い、本形態における構成・処理の詳細について説明する。
［前処理］
まず、前処理としてノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorを特定するデータを、メモリ１０（「記憶部」に相当）のパラメータ領域１６に格納する。なおパラメータとしては、例えば、thnoise=0.01、threv =0.2、thcor=0.5を例示できる。ただし、実際の測定時において、鳴っている音源数が分かるサンプルがあれば、その観測データをもとに各パラメータを調整していってもよい。具体的には、例えば、第１パラメータthnoiseを、ノイズ信号の正規化パワー値よりも大きく残響信号の正規化パワー値よりも小さくなるように調整し、第２パラメータthrevを、源信号の正規化パワー値よりも小さく残響信号の正規化パワー値よりも大きくなるように調整し、第３パラメータthcorを、源信号と残響信号とのエンベロープ相関値より小さくなるように調整する。なお、正規化パワー値やエンベロープ相関値の意味については後述する。
また、信号源数の推定対象となる時間領域の観測信号xj(t)(j={1,...,M})をメモリ１０の観測信号領域１１に書き込む。なお、この観測信号xj(t)はＭ個のセンサ（マイクロホン等）での観測信号であり、下付添字のjは、その観測信号xj(t)がj番目のセンサで観測されたことを示す。
【００１９】
［周波数領域への変換］
まず、制御部５０（図１）が変数jに１を代入し、それをレジスタ５１に格納する（ステップＳ１）。次に、周波数領域変換部２０が、メモリ１０の観測信号領域１１にアクセスし、観測信号xj(t)を読み込む（ステップＳ２）。観測信号xj(t)を読み込んだ周波数領域変換部２０は、それを周波数毎の時間系列データXj(f,τ)に変換してメモリ１０の周波数毎の時系列データ領域１２に格納する（ステップＳ３）。なお、この例では、サンプリング周波数fs、Ｌ点の短時間離散フーリエ変換を利用してこの変換を行う。
次に制御部５０は、レジスタ５１に格納された変数jがＭか否かを判断する（ステップＳ４）。ここでj＝Ｍでないと判断された場合、制御部５０がjに１を加算した値を新たなjとし（ステップＳ５）、それをレジスタ５１に格納してステップＳ２の処理に戻る。一方、j＝Ｍであると判断された場合、以下の信号源推定処理に移る。
【００２０】
［信号源推定処理］
まず、制御部５０（図１）が、変数ｆに０を代入してレジスタ５１に格納する（ステップＳ６）。
独立成分分析（ＩＣＡ：Independent Component Analysis）処理：
次に、独立成分分析（ＩＣＡ）部３１ａ（図２（ａ））が、メモリ１０の周波数毎の時系列データ領域１２から時系列データXj(f,τ)を抽出し、独立成分分析（ＩＣＡ）を用い、X(f,τ)=[X1(f,τ),...,XM(f,τ)]Tから、Ｍ×Ｍ行列の分離行列Ｗ(f)とＩＣＡ分離信号Z(f,τ)=[Z1(f,τ),...,ZM(f,τ)]Tとを生成してレジスタ３１ｂ（「記憶部」に相当）に格納する（ステップＳ７）。ここでＩＣＡによる信号分離は、ＩＣＡ分離信号Z(f,τ)の各要素が互いに独立になるようにZ(f,τ)=Ｗ(f)・X(f,τ)となるＷ(f) を算出する手法である。また、ＩＣＡのアルゴリズムは、A. Hyvarinen and J. Karhunen and E. Oja, "Independent Component Analysis," John Wiley & Sons, 2001, ISBN 0-471-40540,などに様々なものが示されている。
【００２１】
なおＩＣＡの解にはスケーリングの任意性がある。Z(f,τ)のある要素にあるスカラ値を掛けても、要素間の独立性は変化しないからである。従って、この段階では、センサで観測された源信号のパワーがＩＣＡ分離信号に正しく反映されていない可能性が高い。また、源信号の数Ｎがセンサ数Ｍより少なければ、ＩＣＡ分離信号Z(f,τ)のＮ個の要素は源信号に対応し、残りのM-N個の要素はノイズや残響成分に対応するが、この段階のノイズや残響に対応する要素の大きさは一般に増幅されている。そこで、次にスケーリング部３１ｃ（図２（ａ））において、このスケーリングの任意性の問題を解決する。
【００２２】
スケーリング問題解決処理：
スケーリング部３１ｃでは、スケーリングの任意性の問題を解決するため、以下に示す操作を行う。まず、対角行列生成部３１ｃａが、レジスタ３１ｂから分離行列Ｗ(f)を読み出し、この分離行列Ｗ(f)からスケーリング問題を解決するための対角行列Λ(f)を生成する（ステップＳ８）。この対角行列Λ(f)としては、例えば、
Λ(f)=sqrt(diag[(W(f)・W(f)H)-1]) …(4)
が例示できる。ここで、・-1は逆行列、・Hは共役転置行列、diagは対角成分以外を０にする操作、sqrtは各要素の平方根を計算する操作である。
【００２３】
生成された対角行列Λ(f)は、積演算部３１ｃｂに送られ、積演算部３１ｃｂは、これとレジスタ３１ｂから読み出したＩＣＡ分離信号Z(f,τ)とを用い、[Y1(f,τ),...,YM(f,τ)]T←Λ(f)・[Z1(f,τ),...,ZM(f,τ)]Tの演算によって、スケーリング問題を解決した（パワーを回復した）分離信号Yi(f,τ)(i={1,...,M})を生成し、メモリ１０の分離信号領域１３（図１）に格納する（ステップＳ９）。
ここで、式(4)を含む上記一連の操作により、分離信号Yi(f,τ)は、以下の２つの性質を持つ。第一に分離信号Yi(f,τ)が互いに無相関であれば、
【数５】
が成り立つ。すなわち、分離信号Yi(f,τ)のパワーの総和とセンサでの観測信号Xj(f,τ)のパワーの総和が等しくなる。さらに、分離信号Yi(f,τ)が互いに独立であれば、
【数６】
【００２４】
が成り立つ。なお、Sk(f,τ)(k={1,...,N})は源信号成分を示す。すなわち、ある分離信号Yi(f,τ)のパワーと、それに対応する源信号Sk(f,τ)をすべてのセンサで観測した際のパワーの総和とは等しくなる。分離信号Yi(f,τ)が互いに無相関、さらには互いに独立になることは、独立成分分析の目的であり、多くの場合この条件はほぼ満たされている。従って、上記一連の操作により、各分離信号Yi(f,τ)のパワーは、それに対応する源信号Sk(f,τ)がセンサで観測された際のパワーの総和に近くなる。
なお、式(4)の対角行列Λ(f)の代わりに、対角行列Λ(f)=diag[W(f)-1]を使用してもよく、より一般的にW(f)-1のｉ列ｊ行目の要素をｊ行目の対角成分とする対角行列Λを使用してもよい。この場合、各分離信号Yi(f,τ)のパワーは、対応する源信号Sk(f,τ)をあるセンサｊで観測したパワー、すなわち|Hjk(f)・Sk(f,τ)|2に近似する。
【００２５】
［判定処理］
判定処理では、スケーリング問題を解決した（パワーを回復した）分離信号Yi(f,τ)から、源信号の数を推定する。まず、制御部５０（図１）が変数ｉに１を代入し、レジスタ５１に格納する（ステップＳ１０）。
次に、パワー算出部３２の平均パワー算出部３２ａ（図２（ｂ））が、例えばメモリ１０の分離信号領域１３（図１）から各τに対する分離信号Yi(f,τ)を順次抽出し、そのパワー値｜Yi(f,τ)｜2を順次算出してレジスタ３２ｂに格納する。そして、平均パワー算出部３２ａは、レジスタ３２ｂに格納されたパワー値｜Yi(f,τ)｜2を読み出し、分離信号Yi(f,τ)の時間τに関する平均パワー値
σi2(f)←〈｜Yi(f,τ)｜2〉τ
を算出して、レジスタ３２ｂ（「記憶部」に相当）に格納する（ステップＳ１１）。
【００２６】
次に、エンベロープ相関算出部３３のエンベロープ算出部３３ａ（図３（ａ））が、例えばメモリ１０の分離信号領域１３（図１）から各τに対する分離信号Yi(f,τ)を順次抽出し、その絶対値｜Yi(f,τ)｜を順次算出してレジスタ３３ｂに格納する。次に、エンベロープ算出部３３ａは、レジスタ３２ｂに格納された絶対値｜Yi(f,τ)｜を読み出し、時間τに関する平均が０になるように分離信号Yi(f,τ)の絶対値｜Yi(f,τ)｜を正規化したエンベロープ
ｖｉ(f,τ)←｜Yi(f,τ)｜-〈｜Yi(f,τ)｜〉τ
を算出してレジスタ３３ｂ（「記憶部」に相当）に格納する（ステップＳ１２）。
【００２７】
次に、制御部５０（図１）が、レジスタ５１に格納された変数ｉがＭであるか否かを判断する（ステップ１３）。ここでｉ＝Ｍでなければ、制御部５０がｉに１を加算し、その値を新たなｉとしレジスタ５１に格納し（ステップＳ１４）、ステップＳ１１に戻る。一方、ｉ＝Ｍであれば、制御部５０は、この変数ｉに１を代入してレジスタ５１に格納し（ステップＳ１５）、以下の処理を実行する。
まず、パワー算出部３２のパワー正規化部３２ｃ（図２（ｂ））が、レジスタ３２ｂから平均パワー値σ12(f),...,σM2(f)を抽出し、平均パワー値σi2(f)を正規化した正規化パワー値
【数７】
を算出して、メモリ１０のパワー値領域１４に格納する（ステップＳ１６）。
【００２８】
次に、制御部５０（図１）が、変数ｋに１を代入し、レジスタ５１に格納する（ステップＳ１７）。次に、エンベロープ相関算出部３３の相関算出部３３ｃ（図３（ａ））が、レジスタ３３ｂからエンベロープvi(f,τ)及びvk(f,τ)を抽出する。そして、相関算出部３３ｃは、これらのエンベロープvi(f,τ)及びvk(f,τ)を用い、時間差Δτ（例えばL/2やL/4）による分離信号Yi(f,τ)の分離信号Yk(f,τ)とのエンベロープ相関値
【数８】
を算出し、その演算結果Cori,k(f)をレジスタ３３ｂに格納する（ステップＳ１８）。なお、この例のΔτは、例えばプログラムコードに組み込まれた定数である。
【００２９】
次に、制御部５０（図１）はレジスタ５１に格納された変数ｋがＭであるか否かを判断する（ステップＳ１９）。ここで、ｋ＝Ｍでなかった場合、制御部５０がｋに１を加算し、その値を新たなｋとしてレジスタ５１に格納し、ステップＳ１８に戻る（ステップＳ２０）。一方、ｋ＝Ｍであった場合、エンベロープ相関算出部３３の最大値算出部３３ｄ（図３（ａ））は、レジスタ３３ｂ（図３（ａ））から、エンベロープ相関値Cori,1(f),...,Cori,M(f)を抽出する。そして、最大値算出部３３ｄは、これらを用い、エンベロープ相関値Cori,k(f)のｉごとの最大値maxCori(f)を算出し、メモリ１０のエンベロープ領域１５（図１）に格納する（ステップＳ２１）。
【００３０】
次に、判定部３４の比較部３４ａ（図３（ｂ））が、メモリ１０のパワー値領域１４、エンベロープ領域１５及びパラメータ領域１６から、平均パワー値の正規化値NPi(f)、エンベロープ相関値の最大値maxCori(f)並びに第１パラメータthnoise、第２パラメータthrev及び第３パラメータthcorを読み出す。そして、比較部３４ａは、以下の論理式により、分離信号Yi(f,τ)が、源信号に対応するか、ノイズや残響成分に対応するかを判定する（ステップＳ２２）。
【数９】
【００３１】
すなわち、ここでは３種類のパラメータthnoise、threv、thcorを用いている。そして、平均パワー値の正規化値NPi(f)が第１パラメータthnoise未満であればノイズ成分と判定し、平均パワー値の正規化値NPi(f)が第２パラメータthrev未満であり、さらにエンベロープ相関値の最大値maxCori(f)が第３パラメータthcorを超えれば残響成分と判定する。結局、sigi(f)が0になれば、分離信号Yi(f,τ)がノイズや残響成分に対応する（源信号成分でない）と判定されたことになり、sigi(f)=1になれば、分離信号Yi(f,τ)が源信号に対応すると判定されたことになる。そして、このように生成された判定結果sigiはレジスタ３４ｂ（図３（ｂ））に送られて格納される。なお、上記論理式中の「＜」の少なくとも一部を「≦」としてもよく、「＞」を「≧」としてもよい。
【００３２】
次に、制御部５０はレジスタ５１（図１）に格納されている変数ｉがＭであるか否かを判断する（ステップＳ２３）。ここで、ｉ＝Ｍでなければ、制御部５０がｉに１を加算し、その値を新たなｉとしてレジスタ５１に格納してステップＳ１６に戻る（ステップＳ２４）。一方、ｉ＝Ｍであれば、判定部３４の信号源数算出部３４ｃ（図３（ｂ））が、レジスタ３４ｂから判定結果sig1(f),...,sigM(f)を抽出し、信号源数推定値
EN(f)=Σi sigi(f)
を算出し、それをメモリ１０の信号源数領域１７（図１）に格納する（ステップＳ２５）。
【００３３】
次に、制御部５０は、レジスタ５１に格納された変数ｆが{(L-1)/L}fs（fsはサンプリング周波数）であるか否かを判断する（ステップＳ２６）。ここで、変数ｆが{(L-1)/L}fsでなかった場合、制御部５０が変数ｆにfs/Lを加算し、その値を新たな変数ｆとし、レジスタ５１に格納してステップＳ７の処理に戻る（ステップＳ２７）。一方、変数ｆが{(L-1)/L}fsであった場合、以下の結果統合処理を行う。
［結果統合処理］
まず結果統合部４０が、メモリ１０の信号源数領域１７から各周波数fで推定された信号源数推定値EN(0),...,EN({(L-1)/L}fs)を読み出し、これを元に、全体としての信号源数の推定値enを算出して出力する（ステップＳ２８）。この例では、単純に多数決で全体の推定値enを決定する。信頼できる周波数（例えば、高い周波数）に大きな重みを与えて、重みづけの多数決で全体の推定値enを決定しても良い。
【００３４】
［適用結果］
本形態の信号源数の推定方法を音源数の推定に適用した結果を示す。
図１０に一般的な実験条件を例示する。この実験条件は以下である。
・信号源：７秒間の音声
・残響時間：ＴＲ＝２００ ｍｓ
・背景ノイズパワー：−２１．８ ｄＢ
・サンプリング周波数：ｆｓ＝８０００ Ｈｚ
・部屋の大きさ：４．４５ ｍ×３．５５ ｍ×２．５０ ｍ
・音源数：１〜３個
・音源配置・間隔：４ｃｍの間隔で直線上に配置
・センサの数：３個
・中心音源と各センサとの距離：１．１ ｍ
・中心音源と各センサを結んだ直線と、各センサが配置される直線とがなす角度：４５°，９０°，１２０°
この図１０に示す条件で１〜３個の音源を鳴らし、３個のマイクでの観測信号を用いて鳴っている音源の数を推定した。
【００３５】
図６に従来手法と本形態の手法とによる推定結果の比較を示す。ここで、図６（ａ）は固有値に基づく従来手法による信号源数の推定結果を示しており、図６（ｂ）は本形態の手法よる信号源数の推定結果を示している。また、横軸は真の音源数、縦軸は音源数0,1,2,3としてそれぞれ推定した周波数ビンの数を示す。
この図に示すように、従来手法では、１音源や３音源の場合にも多数決によると２音源と推定してしまっている。このように従来手法が推定を誤る原因は、図１１を用いて説明したように、個々の音源やノイズのパワーが固有値に適切に現れていないことや、残響の影響を考慮されていないことである。一方、本形態の手法によるとすべての場合に正しく推定されている。
【００３６】
次に本形態の手法による推定が正確である理由を示す。
まず、パワーの回復（各分離信号が、各信号のパワーを適切に反映しているか）に関して考察する。図７は、３音源の場合にセンサで観測された真のパワー値（図７（ａ））と、これらの混合音を本形態の手法により分離した分離信号（図７（ｂ））のパワー値との比較を示すものである。なお、図７（ａ）の観測結果は、各音源を１つずつ鳴らして測定し、その結果を正規化したものである。これらの図に示すように、本形態の手法による各分離信号のパワー値は、各音源の観測値の真のパワー値に近似し、音源数を推定できる程度に正しくパワーが回復されていることがわかる。
【００３７】
次に、残響の影響ヘの対処に関して考察する。図８は、１音源の場合の１番目と２番目（i=1,2）の分離信号のパワー値（図８（ａ））とそれらのエンベロープの相関値（図８（ｂ））を示すものである。図８（ａ）のパワー値だけを見ると、２番目の分離信号のパワー値が決して十分には小さくないので、信号源なのか残響を含むノイズなのか判断し難い。しかし、右側に示す１番目と２番目の分離信号のエンベロープの相関値を見ると、その値が十分に大きいため、２番目の分離信号は１番目の分離信号の残響成分を多く含むノイズであることがわかる。すなわち、エンベロープの相関値は−１〜＋１の値をとり、信号間の相関性が低いほど０に近づく。図８（ｂ）の例では、エンベロープの相関値が０．６〜１の間に集中しており、１番目の分離信号と２番目の分離信号の相関性が高いことが分かる。そしてパワー値が弱い２番目の分離信号が１番目の信号の残響成分であることが推定できる。
【００３８】
そして、これらの判断に必要なノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorを適切に設定することにより、ノイズや残響の影響が無視できない実環境において、アクティブな源信号の数を精度良く推定することができる。
〔第２の実施の形態〕
次に、本発明における第２の実施の形態について説明する。
本形態は第１の実施の形態の変形例であり、ＩＣＡを用いた信号分離の代わりに固有値に基づく信号分離を行う形態である。以下では、第１の実施の形態との相違点を中心に説明を行い、第１の実施の形態と共通する事項については説明を省略する。
【００３９】
図９（ａ）は、本形態における信号分離部１３１の構成を例示したブロック図である。
なお、本形態の推定装置と第１の実施の形態の推定装置１との相違点は、信号分離部３１が信号分離部１３１になる点のみである。また、本形態の処理と第１の実施の形態の処理との相違点は、信号分離処理（図４：ステップＳ７〜９）と平均パワー算出処理（図４：ステップＳ１１）のみである。
【００４０】
［信号分離処理］
図９（ｂ）は、本形態の信号分離処理を説明するためのフローチャートである。
まず、信号分離部１３１の相関行列生成部１３１ａ（図９（ａ））が、メモリ１０の周波数毎の時系列データ領域１２（図１）から時系列データXj(f,τ)を順次抽出し、時系列ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tに対する相関行列R(f)=〈X(f,τ)・X(f,τ)H〉τを生成する（ステップＳ３１）。
生成された相関行列R(f)は固有値分解部１３１ｂ（図９（ａ））に送られ、固有値分解部１３１ｂはこの相関行列R(f)を、R(f)=V(f)・Λ(f)・V(f)Hの積に分解する（ステップＳ３２）。なお、V(f)=[v1(f),v2(f),...,vM(f)]とし、Λ(f)をλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列とし、vj(f)を固有ベクトルとし、λj(f)をこれに対応する固有値とする。生成された固有値λ1(f)，λ2(f)，...，λM(f)は対応するτに関連つけてメモリ１０（図１）に格納され（ステップＳ３２）、V(f)は積演算部１３１ｃ（図９（ａ））に送られる。
【００４１】
積演算部１３１ｃは、メモリ１０の周波数毎の時系列データ領域１２（図１）から時系列データXj(f,τ)を抽出し、[Y1(f,τ),...,YM(f,τ)]T=V(f)H・[X1(f,τ),...,XM(f,τ)]Tの演算によって、分離信号Yi(f,τ)(i={1,...,M})を生成してメモリ１０の分離信号領域１３に格納する（ステップＳ３３）。なお、〈｜Yi(f,τ)｜2〉τ=λi(f)が成立する。
【００４２】
［平均パワー算出処理］
本形態では、第１の実施の形態のステップＳ１１において、パワー算出部３２の平均パワー算出部３２ａ（図２（ｂ））が、分離信号Yi(f,τ)からパワー値｜Yi(f,τ)｜2を算出し、平均パワー値σi2(f)←〈｜Yi(f,τ)｜2〉τを算出してレジスタ３２ｂに格納していた代わりに、平均パワー算出部３２ａがメモリ１０（図１）から固有値λi(f)を順次抽出し、分離信号Yi(f,τ)の時間τに関する平均パワー値
σi2(f)←λi(f)
を算出して、レジスタ３２ｂに格納する。
なお、その他の処理については第１の実施の形態と同様である。
【００４３】
以上のような構成の場合、パワーが適切に回復されない問題（各信号のパワーが固有値に適切に現れていない問題）は解決されないが、第１の実施の形態と同様、残響の影響の問題は解決できる。そのため、パワーが適切に回復されない問題の影響が少ない周波数領域では、本形態でも正確な信号源数の推定ができる。
なお、本発明は上述の各実施の形態に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００４４】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００４５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００４６】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【００４７】
本発明の音信号に対する応用例としては、例えば、適応ビームフォーマやブラインド音源分離の前処理において、ある区間でのアクティブな音源数を推定する処理を例示できる。
【図面の簡単な説明】
【００４８】
【図１】第１の実施の形態における推定装置の全体を示すブロック図。
【図２】（ａ）は図１に例示した信号分離部の機能構成を例示したブロック図。（ｂ）はパワー算出部の機能構成を例示したブロック図。
【図３】（ａ）はエンベロープ相関算出部の機能構成を例示したブロック図。（ｂ）は判定部の機能構成を例示したブロック図。
【図４】第１の実施の形態における信号源数の推定方法を説明するためのフローチャート。
【図５】第１の実施の形態における信号源数の推定方法を説明するためのフローチャート。
【図６】（ａ）は固有値に基づく従来手法による信号源数の推定結果を示したグラフ。（ｂ）は本形態の手法よる信号源数の推定結果を示したグラフ。
【図７】（ａ）は３音源の場合にセンサで観測された真のパワー値を示した図。（ｂ）は混合音を第１の実施の形態の手法により分離した分離信号のパワー値を示した図。
【図８】（ａ）は１音源の場合の１番目と２番目の分離信号のパワー値を示した図。（ｂ）は、それらのエンベロープの相関値を示した図。
【図９】（ａ）は、第２の実施の形態における信号分離部の構成を例示したブロック図。（ｂ）は第２の実施の形態における信号分離処理を説明するためのフローチャート。
【図１０】実験条件を示した図。
【図１１】固有値に基づく方法でのパワー推定値を示した図。（ａ）は、１つの音源だけを鳴らした場合の各周波数における固有値の正規化パワー値。（ｂ）は、３音源すべてを鳴らした場合の各周波数における固有値の正規化パワー値。
【符号の説明】
【００４９】
１ 推定装置 ３１ 信号分離部
１０ メモリ ３２ パワー算出部
２０ 周波数領域変換部 ３３ エンベロープ相関算出部
３０ 信号減数推定部 ３４ 判定部
【特許請求の範囲】
【請求項１】
観測信号から信号源の数を推定する信号源数の推定方法であって、
複数のパラメータを特定するデータが記憶部に格納されており、
周波数領域変換部が、M個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換して記憶部に格納する手順と、
信号分離部が、上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する手順と、
パワー算出部が、上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順と、
エンベロープ相関算出部が、異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納する手順と、
判定部が、上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項２】
請求項１記載の信号源数の推定方法であって、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)を生成して記憶部に格納する手順は、
独立成分分析部が、独立成分分析を用い、上記時系列データX1(f,τ),...,XM(f,τ)からＭ×Ｍ行列の分離行列Ｗ(f)とＩＣＡ分離信号[Z1(f,τ),...,ZM(f,τ)]Tとを生成して記憶部に格納する手順と、
対角行列生成部が、上記分離行列Ｗ(f)から、スケーリング問題を解決するための対角行列Λ(f)を生成する手順と、
積演算部が、[Y1(f,τ),...,YM(f,τ)]T←Λ(f)・[Z1(f,τ),...,ZM(f,τ)]Tの演算によって、上記分離信号Yi(f,τ)を生成して記憶部に格納する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項３】
請求項１記載の信号源数の推定方法であって、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)を生成して記憶部に格納する手順は、
相関行列生成部が、時系列ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tに対する相関行列R(f)←〈X(f,τ)・X(f,τ)H〉τを生成する手順と、
固有値分解部が、上記相関行列R(f)を、V(f)=[v1(f),v2(f),...,vM(f)]とし、Λ(f)をλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列とし、vj(f)を固有ベクトルとし、λj(f)をこれに対応する固有値とした場合における、R(f)=V(f)・Λ(f)・V(f)Hの積に分解する手順と、
積演算部が、[Y1(f,τ),...,YM(f,τ)]T←V(f)H・[X1(f,τ),...,XM(f,τ)]Tの演算によって、上記分離信号Yi(f,τ)を生成して記憶部に格納する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項４】
請求項１記載の信号源数の推定方法であって、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順は、
平均パワー算出部が、各分離信号Yi(f,τ)の時間τに関する平均パワー値を算出する手順と、
パワー正規化部が、上記平均パワー値を正規化する手順と、を有し、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記平均パワー値の正規化値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順である、
ことを特徴とする信号源数の推定方法。
【請求項５】
請求項１記載の信号源数の推定方法であって、
上記エンベロープ相関値を算出して記憶部に格納する手順は、
エンベロープ算出部が、時間τに関する平均が０になるように上記の各分離信号Yi(f,τ)の絶対値｜Yi(f,τ)｜を正規化したエンベロープｖｉ(f,τ)を算出して記憶部に格納する手順と、
相関算出部が、
【数１】
の演算によってエンベロープ相関値Cori,k(f)を算出して記憶部に格納する手順と、
最大値算出部が、ｉごとに上記エンベロープ相関値Cori,k(f)の最大値maxCori(f)を算出する手順と、を有し、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記の各分離信号Yi(f,τ)のパワー値及び上記最大値maxCori(f)と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順である、
ことを特徴とする信号源数の推定方法。
【請求項６】
請求項１記載の信号源数の推定方法であって、
上記複数のパラメータは、
ノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorであり、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記分離信号Yi(f,τ)の上記パワー値が上記第１パラメータthnoise以下若しくは未満である場合、又は上記分離信号Yi(f,τ)の上記パワー値が上記第２パラメータthrev以下若しくは未満であって上記エンベロープ相関値が上記第３パラメータthcor以上若しくは超える場合、当該分離信号Yi(f,τ)は源信号成分でないと判断する手順を有する、
ことを特徴とする信号源数の推定方法。
【請求項７】
観測信号から信号源の数を推定する推定装置であって、
複数のパラメータを特定するデータが格納された記憶部と、
M個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換して記憶部に格納する周波数領域変換部と、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納するする信号分離部と、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納するパワー算出部と、
異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納するエンベロープ相関算出部と、
上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する判定部と、
を有することを特徴とする推定装置。
【請求項８】
Ｍ個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換する手順と、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する手順と、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順と、
異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納する手順と、
上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、記憶部に格納された各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順と、
をコンピュータに実行させるための推定プログラム。
【請求項９】
請求項８記載の推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１】
観測信号から信号源の数を推定する信号源数の推定方法であって、
複数のパラメータを特定するデータが記憶部に格納されており、
周波数領域変換部が、M個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換して記憶部に格納する手順と、
信号分離部が、上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する手順と、
パワー算出部が、上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順と、
エンベロープ相関算出部が、異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納する手順と、
判定部が、上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項２】
請求項１記載の信号源数の推定方法であって、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)を生成して記憶部に格納する手順は、
独立成分分析部が、独立成分分析を用い、上記時系列データX1(f,τ),...,XM(f,τ)からＭ×Ｍ行列の分離行列Ｗ(f)とＩＣＡ分離信号[Z1(f,τ),...,ZM(f,τ)]Tとを生成して記憶部に格納する手順と、
対角行列生成部が、上記分離行列Ｗ(f)から、スケーリング問題を解決するための対角行列Λ(f)を生成する手順と、
積演算部が、[Y1(f,τ),...,YM(f,τ)]T←Λ(f)・[Z1(f,τ),...,ZM(f,τ)]Tの演算によって、上記分離信号Yi(f,τ)を生成して記憶部に格納する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項３】
請求項１記載の信号源数の推定方法であって、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)を生成して記憶部に格納する手順は、
相関行列生成部が、時系列ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tに対する相関行列R(f)←〈X(f,τ)・X(f,τ)H〉τを生成する手順と、
固有値分解部が、上記相関行列R(f)を、V(f)=[v1(f),v2(f),...,vM(f)]とし、Λ(f)をλ1(f)，λ2(f)，...，λM(f)を対角要素とするＭ行Ｍ列の対角行列とし、vj(f)を固有ベクトルとし、λj(f)をこれに対応する固有値とした場合における、R(f)=V(f)・Λ(f)・V(f)Hの積に分解する手順と、
積演算部が、[Y1(f,τ),...,YM(f,τ)]T←V(f)H・[X1(f,τ),...,XM(f,τ)]Tの演算によって、上記分離信号Yi(f,τ)を生成して記憶部に格納する手順と、
を有することを特徴とする信号源数の推定方法。
【請求項４】
請求項１記載の信号源数の推定方法であって、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順は、
平均パワー算出部が、各分離信号Yi(f,τ)の時間τに関する平均パワー値を算出する手順と、
パワー正規化部が、上記平均パワー値を正規化する手順と、を有し、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記平均パワー値の正規化値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順である、
ことを特徴とする信号源数の推定方法。
【請求項５】
請求項１記載の信号源数の推定方法であって、
上記エンベロープ相関値を算出して記憶部に格納する手順は、
エンベロープ算出部が、時間τに関する平均が０になるように上記の各分離信号Yi(f,τ)の絶対値｜Yi(f,τ)｜を正規化したエンベロープｖｉ(f,τ)を算出して記憶部に格納する手順と、
相関算出部が、
【数１】
の演算によってエンベロープ相関値Cori,k(f)を算出して記憶部に格納する手順と、
最大値算出部が、ｉごとに上記エンベロープ相関値Cori,k(f)の最大値maxCori(f)を算出する手順と、を有し、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記の各分離信号Yi(f,τ)のパワー値及び上記最大値maxCori(f)と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順である、
ことを特徴とする信号源数の推定方法。
【請求項６】
請求項１記載の信号源数の推定方法であって、
上記複数のパラメータは、
ノイズレベルのしきい値を示す第１パラメータthnoise、残響レベルのしきい値を示す第２パラメータthrev、及びエンベロープ相関値のしきい値を示す第３パラメータthcorであり、
上記分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順は、
上記分離信号Yi(f,τ)の上記パワー値が上記第１パラメータthnoise以下若しくは未満である場合、又は上記分離信号Yi(f,τ)の上記パワー値が上記第２パラメータthrev以下若しくは未満であって上記エンベロープ相関値が上記第３パラメータthcor以上若しくは超える場合、当該分離信号Yi(f,τ)は源信号成分でないと判断する手順を有する、
ことを特徴とする信号源数の推定方法。
【請求項７】
観測信号から信号源の数を推定する推定装置であって、
複数のパラメータを特定するデータが格納された記憶部と、
M個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換して記憶部に格納する周波数領域変換部と、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納するする信号分離部と、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納するパワー算出部と、
異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納するエンベロープ相関算出部と、
上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、上記の各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する判定部と、
を有することを特徴とする推定装置。
【請求項８】
Ｍ個のセンサでの観測信号xj(t)(j={1,...,M})を周波数毎の時系列データXj(f,τ)に変換する手順と、
上記時系列データXj(f,τ)から分離信号Yi(f,τ)(i={1,...,M})を生成して記憶部に格納する手順と、
上記の各分離信号Yi(f,τ)のパワー値を算出して記憶部に格納する手順と、
異なる上記分離信号Yi(f,τ)間の時間差Δτに対するエンベロープ相関値を算出して記憶部に格納する手順と、
上記の各分離信号Yi(f,τ)のパワー値及びエンベロープ相関値と、記憶部に格納された各パラメータとを比較し、当該分離信号Yi(f,τ)が源信号成分であるか否かを判断する手順と、
をコンピュータに実行させるための推定プログラム。
【請求項９】
請求項８記載の推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２００６−５８０６５（Ｐ２００６−５８０６５Ａ）
【公開日】平成１８年３月２日（２００６．３．２）
【国際特許分類】
【出願番号】特願２００４−２３８１７４（Ｐ２００４−２３８１７４）
【出願日】平成１６年８月１８日（２００４．８．１８）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成１８年３月２日（２００６．３．２）
【国際特許分類】
【出願日】平成１６年８月１８日（２００４．８．１８）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]