(57)【要約】【課題】汎用性が高く、混合音に基づいて、各音源から発生された音を的確に推定できる音響信号分析装置を提供する。【解決手段】  音響信号分析装置１０は、所定の第１の音を放音する放音手段と、前記放音された第１の音と、前記放音手段とは異なる音源から放音された第２の音とを含む混合音を収音する収音手段と、前記収音された混合音及び前記第１の音に基づいて、前記第２の音と、前記放音手段及び前記収音手段が設置された音場の特性とを同時にベイズ推定する推定手段を備える。【選択図】図４
2015055766.tif 000002
【特許請求の範囲】
【請求項１】
所定の第１の音を放音する放音手段と、
前記放音された第１の音と、前記放音手段とは異なる音源から放音された第２の音とを含む混合音を収音する収音手段と、
前記収音された混合音及び前記第１の音に基づいて、前記第２の音と、前記放音手段及び前記収音手段が設置された音場の特性とを同時にベイズ推定する推定手段を備えたことを特徴とする音響信号分析装置。
【請求項２】
請求項１に記載の音響信号分析装置において、
前記音場の特性は、前記第１の音の各周波数成分の強度に乗算される係数の集合として表され、
前記推定手段は、
前記第２の音のスペクトルの時系列及び前記係数の集合が所定の確率分布に従ってそれぞれ生成されることを表わす生成モデルの事後分布の下限である下限関数であって、前記第２の音及び前記音場の特性に関するパラメータを含む下限関数を設定するとともに前記パラメータを反復的に更新して前記下限関数を決定することにより、前記事後分布を近似的に推定することを特徴とする音響信号分析装置。
【請求項３】
請求項２に記載の音響信号分析装置において、
前記音場の特性に関するパラメータであって、所定の周波数成分の強度に乗算される前記係数に関するパラメータは、前記第１の音の発音開始から現在までの前記第１の音のスペクトルの前記所定の周波数成分の強度の総和及び前記混合音の発音開始から現在までの前記混合音のスペクトルの前記所定の周波数成分の強度の総和にのみ更新回数に応じた重み付け係数が乗算されるように設定された更新式に基づいて更新されることを特徴とする音響信号分析装置。
【請求項４】
請求項１乃至３のうちのいずれか１つに記載の音響信号分析装置において、
前記混合音は、前記第１の音と前記音場の特性とが畳み込まれた音、及び前記第２の音からなることを特徴とする音響信号分析装置。
【請求項５】
請求項１乃至３のうちのいずれか１つに記載の音響信号分析装置において、
前記混合音は、前記第１の音と前記音場の特性とが畳み込まれた音、及び前記第２の音と前記音場の特性とが畳み込まれた音からなることを特徴とする音響信号分析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の音源からそれぞれ発生された複数の音が混合された混合音に基づいて、各音源から発生された音をそれぞれ推定する音響信号分析装置に関する。
【背景技術】
【０００２】
従来から、例えば、下記非特許文献１及び２に示されているように、各音源から発生された音を混合音から抽出する音響信号分析装置は知られている。非特許文献１に記載の音響信号分析装置では、第１の音（話者の声）のモデルと、第２の音（背景音楽）のモデルとを学習しておき、その学習結果を用いて、混合音から第１の音と第２の音を抽出する。
【０００３】
また、非特許文献２に記載の音響信号分析装置では、独立成分分析法を用いて、混合音を外乱とその他の成分に分離する。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Ｈａｋａｎ  Ｅｒｄｏｇａｎ，Ｅｍａｄ  Ｍ．  Ｇｒａｉｓ、「Ｓｅｍｉ−Ｂｌｉｎｄ  Ｓｐｅｅｃｈ−Ｍｕｓｉｃ  Ｓｅｐａｒａｔｉｏｎ  Ｕｓｉｎｇ  Ｓｐａｒｓｉｔｙ  ａｎｄ  Ｃｏｎｔｉｎｕｉｔｙ  Ｐｒｉｏｒｓ」、Ｐａｔｔｅｒｎ  Ｒｅｃｏｇｎｉｔｉｏｎ（ＩＣＰＲ），２０１０，２０ｔｈ  Ｉｎｔｅｒｎａｔｉｏｎａｌ  Ｃｏｎｆｅｒｅｎｃｅ  ｏｎ、Ａｕｇ．  ２０１０、ｐ．４５７３−ｐ．４５７６
【非特許文献２】澤田紘志，Ｊａｎｉ  Ｅｖｅｎ，猿渡洋，鹿野清宏，高谷智哉、「内部雑音抑圧型ロボット音声対話システムにおけるマイクロホンアレー配置の検討」、ＡＩチャレンジ研究会、２００９年  １１月、ｐ．２６−ｐ．３１
【発明の概要】
【０００５】
上記非特許文献１の音響信号分析装置によれば、学習したモデルに関する音からなる混合音にしか対応できない。つまり、非特許文献１の音響信号分析装置は、汎用性に欠ける。
【０００６】
また、非特許文献２の音響信号分析装置によれば、音源の数に応じた複数の収音装置（マイク）が必要である。したがって、装置の構成が複雑になる。
【０００７】
また、上記非特許文献１及び２の音響信号分析装置を用いれば、音響信号分析装置が設置された部屋の音場特性を最尤推定することができる。推定される音場特性の信頼度は、各音源から発生された音に依存する。すなわち、各音源から発生された音に含まれない周波数帯の特性の信頼度は低い。すべての音源から発生された音が既知であるわけではない（少なくとも１つの音源から発生された音は未知である）ので、推定された音場特性のうち、どの周波数帯の信頼度が高く、どの周波数帯の信頼度が低いのかを認定することができない。そのため、各音源から発生された音を的確に推定することが困難である。
【０００８】
本発明は上記問題に対処するためになされたもので、その目的は、汎用性が高く、混合音に基づいて、各音源から発生された音を的確に推定できる音響信号分析装置を提供することにある。なお、下記本発明の各構成要件の記載においては、本発明の理解を容易にするために、実施形態の対応箇所の符号を括弧内に記載しているが、本発明の各構成要件は、実施形態の符号によって示された対応箇所の構成に限定解釈されるべきものではない。
【０００９】
上記目的を達成するために、本発明の特徴は、所定の第１の音を放音する放音手段（１６１）と、前記放音された第１の音と、前記放音手段とは異なる音源から放音された第２の音とを含む混合音を収音する収音手段（１６２）と、前記収音された混合音及び前記第１の音に基づいて、前記第２の音と、前記放音手段及び前記収音手段が設置された音場の特性とを同時にベイズ推定する推定手段（Ｓ１３〜Ｓ１５、Ｓ２３〜Ｓ３１）を備えた音響信号分析装置（１０、２０）としたことにある。
【００１０】
前記音場の特性は、前記第１の音の各周波数成分の強度に乗算される係数（Ｈｆ，ｉ）の集合として表され、前記推定手段は、前記第２の音のスペクトルの時系列及び前記係数の集合が所定の確率分布に従ってそれぞれ生成されることを表わす生成モデルの事後分布の下限である下限関数であって、前記第２の音及び前記音場の特性に関するパラメータを含む下限関数を設定するとともに前記パラメータを反復的に更新して前記下限関数を決定することにより、前記事後分布を近似的に推定するとよい。
【００１１】
この場合、前記混合音は、前記第１の音と前記音場の特性とが畳み込まれた音、及び前記第２の音からなるとよい。
【００１２】
また、この場合、前記混合音は、前記第１の音と前記音場の特性とが畳み込まれた音、及び前記第２の音と前記音場の特性とが畳み込まれた音からなってもよい。
【００１３】
上記のように構成された音響信号分析装置によれば、混合音を構成する音のモデル（又は混合音を構成する音をそれぞれ発生する各音源のモデル）を予め学習しておく必要が無いので、どのような混合音であっても歌声及び音場特性を推定できる。つまり、音響信号分析装置１０は、上記非特許文献１の音響信号分析装置に比べて汎用性が高い。
【００１４】
また、推定される音場特性は、楽音（直接音）の周波数特性に大きく依存するが、本発明によれば、音場特性の事後分布が推定されるので、推定された音場特性の不確かさを認定できる。つまり、推定された音場特性の事後分布の分散が所定の閾値を超える周波数帯域の信頼度は低く、分散が前記所定の閾値以下の周波数帯域の信頼度は高いと認定できる。そして、信頼度が低い周波数帯域の音場特性を、イコライザーなどを用いて補正すれば、混合音から第２の音をより正確に抽出することができる。
【００１５】
また、本発明の他の特徴は、前記音場の特性に関するパラメータであって、所定の周波数成分の強度に乗算される前記係数に関するパラメータは、前記第１の音の発音開始から現在までの前記第１の音のスペクトルの前記所定の周波数成分の強度の総和及び前記混合音の発音開始から現在までの前記混合音のスペクトルの前記所定の周波数成分の強度の総和にのみ更新回数（ｎ）に応じた重み付け係数（ηｎ）が乗算されるように設定された更新式に基づいて更新されることにある。
【００１６】
これによれば、特定の変数についてのみ重みが付されるので、所謂「確定的アニーリング」という手法を採用する場合に比べて、下限関数を決定する際の反復計算の回数が増大することを抑制できる。
【図面の簡単な説明】
【００１７】
【図１】本発明の第１実施形態に係る混合音の生成過程を示すブロック図である。
【図２】図１の音場特性の構成を示すブロック図である。
【図３】本発明の第１及び第２実施形態に係る音響信号分析装置の構成を示すブロック図である。
【図４】音響信号をバッチ処理する場合の手順を示すフローチャートである。
【図５】音響信号を実時間処理する場合の手順を示すフローチャートである。
【図６】本発明の第２実施形態に係る混合音の生成過程を示すブロック図である。
【図７】図６の音場特性の構成を示すブロック図である。
【発明を実施するための形態】
【００１８】
（第１実施形態）
本発明の第１実施形態に係る音響信号分析装置１０について説明する。まず、音響信号分析装置１０の概略について説明する。音響信号分析装置１０は、所定の音（本実施形態では所定の楽曲とする）を放音装置（スピーカ１６１：図３参照）から放音するとともに、収音装置（マイク１６２：図３参照）を用いて音響信号分析装置１０の周囲の音（本実施形態では歌手の歌声とする）を収音する。なお、本実施形態では、放音装置と収音装置とが互いに遠く離れた位置に設置されている。よって、放音装置及び収音装置が設置された部屋の音響的特性（以下、音場特性と呼ぶ）と楽音とが畳み込まれた音も収音装置によって収音される。つまり、放音手段から放音された音（直接音）のみならず、部屋の壁や床などで反射した反射音（残響）も収音される。また、収音装置は歌手に近い位置に設置されている。よって、収音装置で収音される歌声は音場の影響を受けない。本実施形態においては、残響を含む楽音と、歌声（直接音のみ）とが混合された音を混合音と呼ぶ。
【００１９】
混合音のパワースペクトルＹ、楽音（直接音）のパワースペクトルＸ、音場特性Ｈ、及び歌声のパワースペクトルＳの関係は、図１及び図２に示すようなブロック図として表わすことができる。このモデルは、下記の式（１）のように定式化することができる。音響特性Ｈ及び歌声のパワースペクトルＳは直接的には観測されないので、このモデルにおける潜在変数である。なお、楽音が本発明の第１の音に相当し、歌声が本発明の第２の音に相当する。音響信号分析装置１０は、混合音を収録（サンプリング）し、前記収録した混合音を観測データとして、歌声と音場の特性とを同時にベイズ推定する。
【数１】
2015055766.tif 000003
【００２０】
なお、「ｆ」は周波数ビンのインデックス（ｆ＝１，２，・・・，Ｆ）を表わし、「ｔ」は時間フレーム（以下、単にフレームと呼ぶ）のインデックス（ｔ＝１，２，・・・，Ｔ）を表わす。したがって、「Ｘｆ，ｔ」は、楽曲のｔ番目のフレームにおけるｆ番目の周波数ビンの強度（振幅）を表わす。「Ｙｆ，ｔ」は、混合音のｔ番目のフレームにおけるｆ番目の周波数ビンの強度（振幅）を表わす。「Ｓｆ，ｔ」は、歌声のｔ番目のフレームにおけるｆ番目の周波数ビンの強度（振幅）を表わす。また、音場特性Ｈは、係数Ｈｆ，ｉの集合として表わされる。「ｉ」は係数のインデックス（ｉ＝１，２，・・・，Ｉ）を表わす。すなわち、「Ｈｆ，ｉ」は、ｆ番目の周波数ビンの強度（振幅）であってｉ回（ｉ個のフレームの時間）遅延された強度（振幅）に乗算される係数を表わす。
【００２１】
次に音響信号分析装置１０の構成について説明する。音響信号分析装置１０は、図３に示すように、入力操作子１１、コンピュータ部１２、表示器１３、記憶装置１４、外部インターフェース回路１５、及びサウンドシステム１６を備えており、これらがバスＢＳを介して接続されている。
【００２２】
入力操作子１１は、オン・オフ操作に対応したスイッチ（例えば数値を入力するためのテンキー）、回転操作に対応したボリューム又はロータリーエンコーダ、スライド操作に対応したボリューム又はリニアエンコーダ、マウス、タッチパネルなどから構成される。これらの操作子は、演奏者の手によって操作されて、処理開始又は停止、音響信号の分析に関する各種パラメータの設定などに用いられる。入力操作子１１を操作すると、その操作内容を表す操作情報が、バスＢＳを介して、後述するコンピュータ部１２に供給される。
【００２３】
コンピュータ部１２は、バスＢＳにそれぞれ接続されたＣＰＵ１２ａ、ＲＯＭ１２ｂ及びＲＡＭ１２ｃからなる。ＣＰＵ１２ａは、混合音に基づいて歌声及び音場の特性を推定する手順を表わしたプログラムをＲＯＭ１２ｂから読み出して実行する。ＲＯＭ１２ｂには、前記プログラムに加えて、初期設定パラメータ、表示器１３に表示される画像を表わす表示データを生成するための図形データ及び文字データなどの各種データが記憶されている。ＲＡＭ１２ｃには、前記プログラムの実行時に必要なデータが一時的に記憶される。例えば、後述するマイク１６２で収音された混合音を所定のサンプリング周期（例えば１／４４１００ｓｅｃ）でサンプリングして得られた複数のサンプル値からなる混合音データがＲＡＭ１２ｃに記憶される。
【００２４】
表示器１３は、液晶ディスプレイ（ＬＣＤ）によって構成される。コンピュータ部１２は、図形データ、文字データなどを用いて表示すべき内容を表わす表示データを生成して表示器１３に供給する。表示器１３は、コンピュータ部１２から供給された表示データに基づいて画像を表示する。
【００２５】
また、記憶装置１４は、ＨＤＤ、ＦＤＤ、ＣＤ、ＤＶＤなどの大容量の不揮発性記録媒体と、同各記録媒体に対応するドライブユニットから構成されている。記憶装置１４には、前記所定の楽曲を表わす楽曲データが記憶されている。楽曲データは、前記所定の楽曲の演奏を所定のサンプリング周期（例えば１／４４１００ｓｅｃ）でサンプリングして得られた複数のサンプル値からそれぞれなり、各サンプル値が記憶装置１４における連続するアドレスに順に記録されている。楽曲データには、楽曲のタイトルを表わすタイトル情報、容量を表わすデータサイズ情報なども含まれている。楽曲データは予め記憶装置１４に記憶されていてもよいし、後述する外部インターフェース回路１５を介して外部から取り込んでもよい。
【００２６】
外部インターフェース回路１５は、音響信号分析装置１０を電子音楽装置、パーソナルコンピュータなどの外部機器に接続可能とする接続端子を備えている。音響信号分析装置１０は、外部インターフェース回路１５を介して、ＬＡＮ（Ｌｏｃａｌ  Ａｒｅａ  Ｎｅｔｗｏｒｋ）、インターネットなどの通信ネットワークにも接続可能である。
【００２７】
サウンドシステム１６は、楽曲データをアナログ音信号に変換するＤ／Ａ変換器、変換したアナログ音信号を増幅するアンプ、及び増幅されたアナログ音信号を音響信号に変換して放音するスピーカ１６１を備えている。また、サウンドシステム１６は、混合音を収音するためのマイク１６２、及び収音されたアナログ音信号としての混合音をデジタル音信号に変換するＡ／Ｄ変換器も備えている。なお、上記非特許文献２の音響信号分析装置においては複数のマイクを備えているが、本実施形態では１つのマイク１６２のみを備えている。
【００２８】
次に、上記のように構成した音響信号分析装置１０の動作（歌声及び音場特性の推定手順）について説明する。図４に示すように、ステップＳ１０にて歌声及び音場特性の推定処理が開始される。次に、ステップＳ１１にて、各種変数（後述する補助変数や事後分布のパラメータなど）が初期化される。次に、ステップＳ１２にて、楽曲データがサウンドシステム１６に供給されてスピーカ１６１から楽曲の放音が開始されるとともに、マイク１６２で収音された混合音のサンプリングが開始される。サンプリングされた混合音データは、ＲＡＭ１２ｃに記憶される。楽曲の放音が終了すると、以下説明するように、ＲＡＭ１２ｃに記憶された混合音データを観測データとして、歌声及び音場特性が同時に（一体的に）ベイズ推定される。
【００２９】
歌声のスペクトル、楽音のスペクトル及び混合音のスペクトル（短時間フーリエ変換）が複素正規分布から生成されるとすると、下記の式（２）乃至（４）で表されるような生成モデルを構築できる。なお、式（２）乃至（４）における「ＧＩＧ（ａ，ｂ，ｃ）」は、母数ａ，ｂ，ｃによって定義される一般化逆ガウス分布を表わす。また、以下の説明においては、必要に応じ、パラメータがいずれの要素（Ｓ、Ｈなど）に関するものであるかを示すために、右上の括弧内に要素を表わす変数名を記載する。例えば、「ａ（Ｈ）」と記載されたパラメータは、音場特性に関するパラメータである。
【数２】
2015055766.tif 000004
【数３】
2015055766.tif 000005
【数４】
2015055766.tif 000006
【００３０】
上記のモデルの事後分布を、変分ベイズ法を用いて計算する。ここで、対数同時分布は下記の式（５）のように表わされる。なお、式（５）においては、定数項を無視している。
【数５】
2015055766.tif 000007
【００３１】
しかし、式（５）に変分ベイズ法を適用することができないので、補助関数を用いて下限を定める。具体的には、下記の式（６）のような下限関数を設定し、新たに導入された補助変数Ｍ及び補助変数Φが更新される。
【数６】
2015055766.tif 000008
【００３２】
具体的には、ステップＳ１３にて、式（７）乃至（１０）によって定義される条件下で補助変数Ｍが最適化され、式（１１）によって定義される条件下で補助変数Φが最適化される。
【数７】
2015055766.tif 000009
【数８】
2015055766.tif 000010
【数９】
2015055766.tif 000011
【数１０】
2015055766.tif 000012
【数１１】
2015055766.tif 000013
【００３３】
次に、ステップＳ１４にて、下記の式（１２）及び式（１３）を用いて事後分布のパラメータが更新される。なお、式（１２）及び式（１３）中の各パラメータは下記の式（１４）乃至（１９）のように定義されている。
【数１２】
2015055766.tif 000014
【数１３】
2015055766.tif 000015
【数１４】
2015055766.tif 000016
【数１５】
2015055766.tif 000017
【数１６】
2015055766.tif 000018
【数１７】
2015055766.tif 000019
【数１８】
2015055766.tif 000020
【数１９】
2015055766.tif 000021
【００３４】
次にステップＳ１５にて、下限関数が収束したか否かが判定される。すなわち、補助変数Ｍ及び補助変数Φ、並びに事後分布の各パラメータが収束したか否かが判定される。下限関数が収束していない場合には「Ｎｏ」と判定され、ステップＳ１３及びステップＳ１４にて補助変数Ｍ及び補助変数Φ、並びに事後分布の各パラメータがそれぞれ更新される。一方、下限関数が収束した場合には、「Ｙｅｓ」と判定され、ステップＳ１６にて、歌声及び音場特性の推定処理が終了する。上記のようにして、補助変数Ｍ及び補助変数Φ、並びに事後分布の各パラメータが反復的に更新されて下限関数が決定されることにより、事後分布が近似的に計算される。これにより、歌声及び音場特性が同時に（一体的に）ベイズ推定される。なお、下記の式（２０）で示されるマスクを、ｔ番目のフレームにおける混合音のスペクトル（短時間フーリエ変換）に適用し、その逆フーリエ変換を計算することにより、混合音から歌声を抽出することができる。
【数２０】
2015055766.tif 000022
【００３５】
上記のように構成された音響信号分析装置１０によれば、混合音を構成する音のモデル（又は混合音を構成する音をそれぞれ発生する各音源のモデル）を予め学習しておく必要が無いので、どのような混合音であっても歌声及び音場特性を推定できる。つまり、音響信号分析装置１０は、上記非特許文献１の音響信号分析装置に比べて汎用性が高い。
【００３６】
また、上記非特許文献２の音響信号分析装置とは異なり、１つのマイクを備えていればよいので、装置の構成が単純である。
【００３７】
また、推定される音場特性は、楽音（直接音）の周波数特性に大きく依存するが、本実施形態では、音場特性の事後分布が推定されるので、推定された音場特性の不確かさを認定できる。つまり、推定された音場特性の事後分布の分散が所定の閾値（例えば、予め定められた値、又は楽音のパワースペクトルＸにおける特定の周波数ビンのパワーに比例した値）を超える周波数帯域の信頼度は低く、分散が前記所定の閾値以下の周波数帯域の信頼度は高いと認定できる。そして、信頼度が低い周波数帯域の音場特性を、イコライザーなどを用いて補正すれば、混合音から歌声をより正確に抽出することができる。
【００３８】
上記実施形態における式（２）乃至（４）によって表わされる生成モデルに代えて、次の式（２１）乃至（２３）によって表わされる生成モデルを採用しても良い。
【数２１】
2015055766.tif 000023
【数２２】
2015055766.tif 000024
【数２３】
2015055766.tif 000025
【００３９】
この場合、Ｐｏｉｓｓｏｎ分布の再生性を用いるとともに、補助変数Ｍ（Ｓ），Ｍ（Ｈ）を下記の式（２４）を満たすように制約すれば、上記の式（２１）乃至（２３）によって表わされる生成モデルは、下記の式（２５）乃至（２８）によって表わされる生成モデルと等価である。
【数２４】
2015055766.tif 000026
【数２５】
2015055766.tif 000027
【数２６】
2015055766.tif 000028
【数２７】
2015055766.tif 000029
【数２８】
2015055766.tif 000030
【００４０】
そして、対数同時分布は、下記の式（２９）のように表わされる。
【数２９】
2015055766.tif 000031
【００４１】
事後分布は、下記の式（３０）乃至（３３）を用いて更新される。
【数３０】
2015055766.tif 000032
【数３１】
2015055766.tif 000033
【数３２】
2015055766.tif 000034
【数３３】
2015055766.tif 000035
【００４２】
なお、上記式（３１）におけるＺｆ，ｔは、下記の式（３４）に示すような正規化係数である。また、その他のパラメータは、下記の式（３５）乃至（４１）に示すように定義される。
【数３４】
2015055766.tif 000036
【数３５】
2015055766.tif 000037
【数３６】
2015055766.tif 000038
【数３７】
2015055766.tif 000039
【数３８】
2015055766.tif 000040
【数３９】
2015055766.tif 000041
【数４０】
2015055766.tif 000042
【数４１】
2015055766.tif 000043
【００４３】
この場合、「Ｓｆ，ｔ」の平均値に混合音の短時間フーリエ変換の位相を付与することにより、混合音から歌声を抽出することができる。
【００４４】
上記実施形態では、「Ｉ」を固定する必要がある。しかし、「Ｉ」は音場の残響時間に依存するので、様々な音場に対して頑健であるためには、「Ｉ」に対する依存性が弱められていることが望ましい。そこで、下記の式（４２）のように定義される変数ｇｊを上記の式（２１）に導入し、下記の式（４３）のようにモデルを定式化しても良い。つまり、変数ｇｊは、現在のフレームの歌声のゲイン、及び過去のフレームの楽音のゲインを表わす。
【数４２】
2015055766.tif 000044
【数４３】
2015055766.tif 000045
【００４５】
これによれば、「Ｉ」の値が音場に応じて動的に変化するような挙動が得られる。これにより、混合音から歌声をより正確に抽出できる。
【００４６】
また、上記第１実施形態及びその変形例では、フレームの総和（「ｔ」に関する総和）を表わす項が含まれる数式を用いているため、歌声及び音場特性を実時間で推定することができない。そこで、式（２１）乃至（４３）を用いて説明した歌声及び音場特性の推定手順を変形し、歌声及び音場特性を実時間で推定する手順について説明する。
【００４７】
まず、事後分布は各周波数ビンに関して独立であるとみなし、１つの周波数ビンに注目する。すると、事後分布は下記の式（４４）及び式（４５）のように表わされる。
【数４４】
2015055766.tif 000046
【数４５】
2015055766.tif 000047
【００４８】
上記の事後分布を、変分ベイズ法を用いて計算する。具体的には、下記の式（４６）に示すような目的関数を最適化することにより、事後分布を近似的に計算する。なお、式（４６）中の「Ｈｑ（θ）」は、下記の式（４７）に示すように定義される。
【数４６】
2015055766.tif 000048
【数４７】
2015055766.tif 000049
【００４９】
ここで、式（４８）に示すように、事後分布における変数のうち、フレームごとに独立である変数を先に最適化する。式（４８）に示す目的関数Ｊ´を最適化することは、式（４６）に示す目的関数Ｊを最適化することと等価である。
【数４８】
2015055766.tif 000050
【００５０】
さらに、「０」から「Ｔ」の間で一様に分布する確率変数τ〜Ｕｎｉｆｏｒｍ（０，Ｔ）を導入し、目的関数Ｊτを下記の式（４９）のように定義する。目的関数Ｊτは、全ての観測データが、フレームτのものであったとみなしたときの、目的関数Ｊ´の目的関数である。
【数４９】
2015055766.tif 000051
【００５１】
なお、観測データは一様に分布しているので、下記の式（５０）が成立する。フレームのインデックスである「ｔ」を「１」ずつ増加させ、各フレームにおいて目的関数Ｊτ（θ）を評価し、その最適値を累積すれば、ＪＩ（θ）の平均値を実時間で更新することができる。
【数５０】
2015055766.tif 000052
【００５２】
上記の歌声及び音場特性の実時間推定処理の手順を、図５を用いて説明する。ステップＳ２０にて歌声及び音場特性の実時間推定処理が開始される。ステップＳ２１にて、各種変数が初期化される。次に、ステップＳ２２にて、楽曲の再生が開始される。次に、ステップＳ２３にて、ｔ番目のフレームに相当する混合音のスペクトルが計算される。具体的には、所定のサンプリング周期でバッファにサンプル値が逐次的に記憶されており、前記バッファに記憶されたサンプル値のうちのｔ番目のフレームに相当する複数のサンプル値を用いて、ｔ番目のフレームに相当する混合音のスペクトルが計算される。最初、「ｔ」は「０」に初期化されているので、最初のフレームに相当するスペクトルが計算される。次に、ステップＳ２４にて、補助変数Φが更新される。次に、ステップＳ２５にて、潜在変数Ｓｔが更新される。次に、ステップＳ２６にて、補助変数Φが収束したか否かが判定される。補助変数Φが収束していない場合、「Ｎｏ」と判定され、ステップＳ２４及びステップ２５を実行し、再び補助変数Φ及び潜在変数Ｓｔが更新される。一方、補助変数Φが収束した場合、「Ｙｅｓ」と判定され、ステップＳ２７にて、潜在変数Ｓｔの平均値に基づいて歌声の音響信号が復元される。そして、ステップＳ２８にて、事後分布のパラメータが更新される。次に、ステップＳ２９にて、事後分布が更新される。ただし、「ρｔ」は、０≦ρｔ≦１を満たす係数である。また、下記の式（５１）を満たす場合には、変数Φの収束が保証される。
【数５１】
2015055766.tif 000053
【００５３】
次に、ステップＳ３０にて、フレームのインデックスが更新（インクリメント）される。次に、ステップＳ３１にて、最終フレームの処理を終了したか否かが判定される。つまり、フレームのインデックス（ｔ）が最終フレームのインデックス（Ｔ）を超えていれば、最終フレームを既に処理したと判定し、ステップＳ３２にて実時間処理が終了する。一方、フレームのインデックス（ｔ）が最終フレームのインデックス（Ｔ）以下であれば、ステップＳ２３乃至Ｓ３１からなる処理が再び実行される。
【００５４】
上記の実時間処理において、「Ｈ」を階層化し、事前分布の超パラメータを最適化してもよい。これによれば、各フレーム間において「Ｈ」が独立であるため、現在のフレームのデータのみから「Ｈ」を更新できる。ただし、このままではサンプル数が少なくなりすぎ、歌声の推定精度が低下する虞がある。そこで、各フレームにて推定される「Ｈ」の事後分布から適切な「Ｈ」の事前分布を更新する。つまり、式（２２）に代えて、下記の式（５２）を用いる。
【数５２】
2015055766.tif 000054
【００５５】
これによれば、「Ｈｆ，ｔ，ｉ」の更新式には、フレームの総和（「ｔ」に関する総和）を表わす項が含まれない。したがって、「Ｈｆ，ｔ，ｉ」は各フレームにおいて独立して更新される。
【００５６】
次に、下記の式（５３）に示す目的関数を最大化することにより、超パラメータを最適化する。なお、任意の「ｆ」及び「ｉ」に関して、式（５４）乃至（５７）が成立する。
【数５３】
2015055766.tif 000055
【数５４】
2015055766.tif 000056
【数５５】
2015055766.tif 000057
【数５６】
2015055766.tif 000058
【数５７】
2015055766.tif 000059
【００５７】
ここで、Ｇ（ａ，ｂ）を「ｂ」に関して偏微分した結果が「０」に等しいとすると、「ｂ＝ａ／ｃ」という関係式が得られる。また、Ｇ（ａ，ｂ）を「ａ」に関して偏微分した結果が「０」に等しいとすると、「ｄ＋ｌｏｇｂ−ψ（ａ）＝０」という関係式が得られる。「ａ」が求まれば「ｂ」も求まるが、「ａ」を求めるのは困難である。そこで、上記の２つの関係式を組み合わせた下記の式（５８）を解く。
【数５８】
2015055766.tif 000060
【００５８】
ここで、ｅｘｐ（ａ´）＝ａとすると、下記の式（５９）に示す目的関数が導出される。
【数５９】
2015055766.tif 000061
【００５９】
この式（５９）にＮｅｗｔｏｎ  Ｒａｐｈｓｏｎ法を適用すると、下記の式（６０）に示す更新式が導出される。
【数６０】
2015055766.tif 000062
【００６０】
つまり、「ａ」は、下記の式（６１）のように表わされる。
【数６１】
2015055766.tif 000063
【００６１】
ディガンマ関数の近似式を用いれば、下記の式（６２）及び式（６３）が成立する。
【数６２】
2015055766.tif 000064
【数６３】
2015055766.tif 000065
【００６２】
したがって、下記の式（６４）に示すような更新式が導出される。この更新式を用いてパラメータを更新しても良い。
【数６４】
2015055766.tif 000066
【００６３】
なお、式（４１）及び式（４２）に鑑みれば、単純に履歴を蓄積していけばよいとも思われる。すなわち、例えば、下記の式（６５）のような更新式を用いれば良いとも思われる。
【数６５】
2015055766.tif 000067
【００６４】
しかし、変分ベイズ法は反復計算を利用するものであり、「Ｓ」の初期値はスパースである（「０」に近い値を優遇する）。一方、観測データが多くなるほど、「ａ（Ｈ）」及び「ｂ（Ｈ）」の確信度が高くなる。すなわち、初期段階では観測データ数が少ないので確信度が低い。また、式（６５）を用いた場合、新しい観測は古い観測に対して弱い影響しか与えない。そこで、下記の式（６６）及び（６７）に示すように、履歴の蓄積に対して重みを付ける変数ηｎを導入する。すなわち、変数ηｎは、楽音（第１の音）の発音開始から現在までの楽音のスペクトルの所定の周波数成分の強度の総和及び混合音の発音開始から現在までの混合音のスペクトルの所定の周波数成分の強度の総和にのみ更新回数に応じて乗算される重み付け係数に相当する。なお、「ｎ」は反復計算の回数を表わす。初期段階（すなわち、「ｎ」が小さいとき）では、変数ηｎを「０」に近い値に設定し、反復回数が増すにつれて変数ηｎを徐々に「１」に近づける。
【数６６】
2015055766.tif 000068
【数６７】
2015055766.tif 000069
【００６５】
これによれば、特定の変数についてのみ重みが付されるので、所謂「確定的アニーリング」という手法を採用する場合に比べて、下限関数を決定する際の反復計算の回数が増大することを抑制できる。
【００６６】
なお、式（２１）乃至（４３）を用いて説明した歌声及び音場特性の推定手順を変形し、歌声及び音場特性を実時間で推定する手順について説明したが、式（２）乃至（２０）を用いて説明した歌声及び音場特性の推定手順に関しても同様に変形し、歌声及び音場特性を実時間で推定することが可能である。
【００６７】
（第２実施形態）
次に、本発明の第２実施形態に係る音響信号分析装置２０について説明する。まず、音響信号分析装置２０の概略について説明する。音響信号分析装置２０は、音響信号分析装置１０と同様に、所定の楽曲を放音装置（スピーカ１６１）から放音するとともに、収音装置（マイク１６２）を用いて歌手の歌声を収音する。ただし、本実施形態では、第１実施形態とは異なり、放音装置と収音装置とが近い位置に設置されている。また、収音装置は歌手及び放音装置から遠く離れた位置に設置されている。よって、放音装置及び収音装置が設置された部屋の音響特性（音場特性）と楽音及び歌声とが畳み込まれた音が収音装置によって収音される。つまり、放音装置から放音された音（直接音）及び歌手の歌声（直接音）のみならず、部屋の壁や床などで反射した反射音（反響）も収音される。本実施形態においては、残響を含む楽曲の音と、残響を含む歌声とが混合された音を混合音と呼ぶ。
【００６８】
混合音のパワースペクトルＹ、楽音（直接音）のパワースペクトルＸ、音場特性Ｈ、及び歌声のパワースペクトルＳの関係は、図６及び図７に示すようなブロック図として表わすことができ、このモデルは、下記の式（６８）のように定式化することができる。音響特性Ｈ及び歌声のパワースペクトルＳは直接的には観測されないので、このモデルにおける潜在変数である。
【数６８】
2015055766.tif 000070
【００６９】
音響信号分析装置２０は、混合音を収録（サンプリング）し、前記収録した混合音を観測データとして、歌声と音場の特性とを同時に（一体的に）ベイズ推定する。音響信号分析装置２０のその他の構成は音響信号分析装置１０と同様であるので、その説明を省略する。
【００７０】
次に、上記のように構成した音響信号分析装置２０の動作（歌声及び音場特性の推定手順）について説明する。第１実施形態と同様に、図４に示すように、ステップＳ１０にて歌声及び音場特性の推定処理が開始される。次に、ステップＳ１１にて、各種変数（後述する補助変数や事後分布のパラメータなど）が初期化される。次に、ステップＳ１２にて、楽曲データがサウンドシステム１６に供給されてスピーカ１６１から楽曲の放音が開始されるとともに、マイク１６２で収音された混合音のサンプリングが開始される。サンプリングされた混合音データは、ＲＡＭ１２ｃに記憶される。楽曲の放音が終了すると、以下説明するように、ＲＡＭ１２ｃに記憶された混合音データを観測データとして、歌声及び音場特性が同時に（一体的に）ベイズ推定される。
【００７１】
歌声のスペクトル、楽音のスペクトル及び混合音のスペクトル（短時間フーリエ変換）が複素正規分布から生成されるとすると、下記の式（６９）乃至（７１）で表されるような生成モデルを構築できる。
【数６９】
2015055766.tif 000071
【数７０】
2015055766.tif 000072
【数７１】
2015055766.tif 000073
【００７２】
上記のモデルの事後分布を、変分ベイズ法を用いて計算する。ここで、対数同時分布は次の式（７２）のように表わされる。なお、式（７２）においては、定数項を無視している。
【数７２】
2015055766.tif 000074
【００７３】
しかし、式（７２）に変分ベイズ法を適用することができないので、補助関数を用いて下限を定める。具体的には、下記の式（７３）のような下限関数を設定し、新たに導入された補助変数Ｍ及び補助変数Φが更新される。
【数７３】
2015055766.tif 000075
【００７４】
具体的には、ステップＳ１３にて、式（７４）乃至（７７）によって定義される条件下で補助変数Ｍが最適化され、式（７８）によって定義される条件下で補助変数Φが最適化される。
【数７４】
2015055766.tif 000076
【数７５】
2015055766.tif 000077
【数７６】
2015055766.tif 000078
【数７７】
2015055766.tif 000079
【数７８】
2015055766.tif 000080
【００７５】
次に、ステップＳ１４にて、下記の式（７９）及び式（８０）を用いて事後分布のパラメータが更新される。なお、式（７９）及び式（８０）中の各パラメータは下記の式（８１）乃至（８６）のように定義されている。
【数７９】
2015055766.tif 000081
【数８０】
2015055766.tif 000082
【数８１】
2015055766.tif 000083
【数８２】
2015055766.tif 000084
【数８３】
2015055766.tif 000085
【数８４】
2015055766.tif 000086
【数８５】
2015055766.tif 000087
【数８６】
2015055766.tif 000088
【００７６】
次にステップＳ１５にて、下限関数が収束したか否かが判定される。すなわち、補助変数Ｍ及び補助変数Φ、並びに事後分布の各パラメータが収束したか否かが判定される。下限関数が収束していない場合には「Ｎｏ」と判定され、ステップＳ１３及びステップＳ１４にて補助変数Ｍ及び補助変数Φ、並びに事後分布の各パラメータがそれぞれ更新される。一方、下限関数が収束した場合には、「Ｙｅｓ」と判定され、ステップＳ１６にて、歌声及び音場特性の推定処理が終了する。上記のようにして、事後分布が近似的に計算される。これにより、歌声及び音場特性が同時に（一体的に）ベイズ推定される。なお、下記の式（８７）で示されるマスクを、ｔ番目のフレームにおける混合音のスペクトル（短時間フーリエ変換）に適用し、その逆フーリエ変換を計算することにより、混合音から歌声を抽出することができる。
【数８７】
2015055766.tif 000089
【００７７】
上記のように構成された音響信号分析装置２０によっても、第１実施形態と同様の効果が得られる。
【００７８】
上記実施形態における式（６９）乃至（７１）によって表わされる生成モデルに代えて、次の式（８８）乃至（９０）によって表わされる生成モデルを採用しても良い。
【数８８】
2015055766.tif 000090
【数８９】
2015055766.tif 000091
【数９０】
2015055766.tif 000092
【００７９】
この場合、Ｐｏｉｓｓｏｎ分布の再生性を用いるとともに、補助変数Ｍ（Ｓ），Ｍ（Ｈ）を下記の式（９１）を満たすように制約すれば、上記の式（８８）乃至（９０）によって表わされる生成モデルは、下記の式（９２）乃至（９５）によって表わされる生成モデルと等価である。
【数９１】
2015055766.tif 000093
【数９２】
2015055766.tif 000094
【数９３】
2015055766.tif 000095
【数９４】
2015055766.tif 000096
【数９５】
2015055766.tif 000097
【００８０】
この場合、対数同時分布は、下記の式（９６）のように表わされる。
【数９６】
2015055766.tif 000098
【００８１】
事後分布は、下記の式（９７）乃至（１００）を用いて更新される。
【数９７】
2015055766.tif 000099
【数９８】
2015055766.tif 000100
【数９９】
2015055766.tif 000101
【数１００】
2015055766.tif 000102
【００８２】
なお、上記式（９８）におけるＺｆ，ｔは、下記の式（１０１）に示すような正規化係数である。また、その他のパラメータは、下記の式（１０２）乃至（１０８）に示すように定義される。
【数１０１】
2015055766.tif 000103
【数１０２】
2015055766.tif 000104
【数１０３】
2015055766.tif 000105
【数１０４】
2015055766.tif 000106
【数１０５】
2015055766.tif 000107
【数１０６】
2015055766.tif 000108
【数１０７】
2015055766.tif 000109
【数１０８】
2015055766.tif 000110
【００８３】
この場合、「Ｓｆ，ｔ」の平均値に混合音の短時間フーリエ変換の位相を付与することにより、混合音から歌声を抽出することができる。
【００８４】
また、第２実施形態においても、第１実施形態の変形例と同様に数式を変形することにより、音場特性及び歌声を実時間推定することができる。
【符号の説明】
【００８５】
１０，２０・・・音響信号分析装置、１６１・・・スピーカ、１６２・・・マイク、Ｙ・・・混合音のパワースペクトル、Ｘ・・・楽音のパワースペクトル、Ｈ・・・音場特性、Ｓ・・・歌声のパワースペクトル
【図１】
2015055766.tif 000111
【図２】
2015055766.tif 000112
【図３】
2015055766.tif 000113
【図４】
2015055766.tif 000114
【図５】
2015055766.tif 000115
【図６】
2015055766.tif 000116
【図７】
2015055766.tif 000117