残響除去方法、残響除去装置、プログラム
【課題】ロバストにブラインド残響しノイズの増幅を抑える残響除去方法を提供する。
【解決方法】Ｎチャネル（Ｎ≧２）の収音信号間の相関関数からＬ×Ｌ相関行列（Ｌは逆フィルタの１チャネルあたりのタップ数）を生成し、当該相関行列からＮＬ×ＮＬ多チャネル相関行列を算出する多チャネル相関行列算出ステップと、多チャネル相関行列を特異値分解して特異値を出力する特異値分解サブステップと、特異値の減衰傾向に基づいて信号部分空間の次元を推定する信号部分空間次元推定サブステップと、信号部分空間の次元に基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定サブステップと、疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定サブステップと、逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解しＮチャネルの各収音信号に適用する逆フィルタステップとを有する。
【解決方法】Ｎチャネル（Ｎ≧２）の収音信号間の相関関数からＬ×Ｌ相関行列（Ｌは逆フィルタの１チャネルあたりのタップ数）を生成し、当該相関行列からＮＬ×ＮＬ多チャネル相関行列を算出する多チャネル相関行列算出ステップと、多チャネル相関行列を特異値分解して特異値を出力する特異値分解サブステップと、特異値の減衰傾向に基づいて信号部分空間の次元を推定する信号部分空間次元推定サブステップと、信号部分空間の次元に基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定サブステップと、疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定サブステップと、逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解しＮチャネルの各収音信号に適用する逆フィルタステップとを有する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は残響とノイズを含む多チャネル収音信号から音源の発する本来の信号を抽出する残響除去方法、残響除去装置、プログラムに関する。
【背景技術】
【０００２】
マイクロホンで音声等を収音する場合、マイクロホンに到達する音は、音源から発せられて直接到達した直接音と、壁等に１回反射して到達した反射音と、複数回反射して到達した反射音の混ざったものになる。マイクロホンが音源に近接していない場合は、この反射音からなる残響成分が多くなる。話者音声を離れた位置で収音する場合や残響時間が長い部屋で収音する場合には、残響により音声の明瞭度が低下して聞き取り難くなり、例えば音声認識システムの音声認識率も著しく低下することが知られている。話者位置からマイクロホン位置までのインパルス応答が予め測定されている場合には、ＭＩＮＴ（Ｍｕｌｔｉ−ｉｎｐｕｔ／ｏｕｔｐｕｔ Ｉｎｖｅｒｓｅ Ｔｈｅｏｒｅｍ）法（非特許文献１）を用いてインパルス応答から逆フィルタ係数を求め、収音信号に逆フィルタを適用することで残響に含まれる収音信号から原音を回復することができる。
【０００３】
ＭＩＮＴ法を簡単に説明する。マイクロホンがＮ個あり、目的音源からｉ番目のマイクロホンまでの音響経路のインパルス応答をｈｉ（０）〜ｈｉ（Ｋ）とする。Ｌタップの逆フィルタを求めるとき、ＭＩＮＴ法は目的インパルス応答ベクトルをｂ、ｉ番目のマイクロホン収音信号に適用するＬ次元の逆フィルタ係数ベクトルをｇｎ（ｎ＝１〜Ｎ）としてサイズ（Ｌ＋Ｋ）×Ｌのインパルス応答畳み込み行列（式（１））を用いて、式（２）若しくは式（３）の連立方程式を解いて逆フィルタ係数ベクトルｇｎを求める。
【０００４】
【数１】
【０００５】
完全な逆フィルタを求めるには、Ｌ＋Ｋ≧ＮＬが満たされている必要がある。Ｌ＋Ｋ<ＮＬの関係を満たすＬをもちいた場合、近似の逆フィルタが求まる。Ｌを小さくしていくと、その近似精度は下がってく。目的インパルス応答ベクトルｂはＬ＋Ｋ次元のベクトルで、無遅延のパルス（式（４））や遅延付きのパルスが（式（５））用いられる。
【０００６】
【数２】
【０００７】
通常は、音源位置（話者位置）からマイクロホン位置までのインパルス応答を予め測定することが出来ないため、残響に含まれる収音信号のみから原音を回復するブラインド残響除去法が用いられる。ブラインド残響除去の一方法としてｂｌｉｎｄ−ＭＩＮＴ法（特許文献１、非特許文献２）がある。
【０００８】
図１、図２を参照してｂｌｉｎｄ−ＭＩＮＴ法を簡単に説明する。図１は従来技術（ｂｌｉｎｄ−ＭＩＮＴ法）に係る残響除去装置９００の構成を示すブロック図である。図２は従来技術（ｂｌｉｎｄ−ＭＩＮＴ法）に係る残響除去装置９００の動作を示すフローチャートである。残響除去装置９００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、逆フィルタ推定部９０と、逆フィルタ部７０とを備える。
【０００９】
目的音源１はスペクトルがフラットな音を発し、Ｎ個のマイクロホン１０−１〜１０−Ｎで収音される。収音信号ｙ１（ｋ）〜ｙＮ（ｋ）は、多チャネル相関行列算出部２０に入力される。ここで、収音信号ｙ１（ｋ）〜ｙＮ（ｋ）は、例えば標本化周波数は１６ｋＨｚで離散値化されたディジタル信号であり、ｋは時刻を表す。なお、収音信号をディジタル化するＡ／Ｄ変換器等については省略している。多チャネル相関行列算出部２０は、収音信号ｙ１（ｋ）〜ｙＮ（ｋ）間の相関関数ｒｉｊ（τ）を式（６）で計算する。
【００１０】
【数３】
【００１１】
ここでτは、２時点の時間間隔である。多チャネル相関行列算出部２０は、相関関数ｒｉｊ（τ）からＬ×Ｌ相関行列Ｒｉｊ（式（７））を生成し、ＮＬ×ＮＬ多チャネル相関行列Ｒ（式（８））を生成する（Ｓ２０）。
【００１２】
【数４】
【００１３】
但し、Ｌは逆フィルタの１チャネル当りのタップ数である。逆フィルタ推定部９０は、式（９）の連立一次方程式を解いて逆フィルタ係数ベクトルｇを推定する（Ｓ９０）。
【００１４】
【数５】
【００１５】
ここで、目的インパルス応答ベクトルｂ（式（１０））はＬＮ次元ベクトルであり、Ｌ次元縦ベクトルｂ１〜ｂＮが縦に並んでいる（式（１１））。
【００１６】
【数６】
【００１７】
Ｌ次元ベクトルｂｎ中のδｎは、第ｎ番目のマイクロホン１０−ｎと目的音源との位置関係で決まる。第ｎ番目のマイクロホン１０−ｎが全マイクロホンの中で目的音源１に最も近いときにδｎは１を取り、それ以外のときは０を取る。逆フィルタ係数ベクトルｇは、ｎ番目のマイクロホンの収音信号ｙｎ（ｋ）に適用するＬ次元の逆フィルタ係数ベクトル（式（１２））を用いると、式（１３）と書ける。
【００１８】
【数７】
【００１９】
逆フィルタ部７０は、式（１４）に示すように各収音信号ｙ１（ｋ）〜ｙＮ（ｋ）に逆フィルタ係数ｇを適用して最終出力信号ｘ＾（ｋ）を得る。この出力信号は、残響の除去された目的信号になる。
【００２０】
【数８】
【００２１】
従来のブラインド残響除去装置９００は、目的音源１から各マイクロホンまでの未知音響経路インパルス応答ｈｉ（０）〜ｈｉ（Ｌ−１）から成る行列Ｈと収音信号から求めた相関行列Ｒとの間に式（１５）の関係があることを利用し、逆フィルタ係数を求めている。αは比例定数（スカラー量）である。
【００２２】
【数９】
【先行技術文献】
【特許文献】
【００２３】
【特許文献１】特開２００６−６６９８９号公報
【非特許文献】
【００２４】
【非特許文献１】大賀寿郎、山崎芳男、金田豊著、「音響システムとディジタル処理」、電子情報通信学会、1995年4月
【非特許文献２】K. Furuya and A. Kataoka, "Robust Speech Dereverberatio Using Multichannel Blind Deconvolution With Spectral Subtraction", IEEE Transactions on Audio, Speech, and Language Processing, 2007 July, vol. 15, no. 5, pp. 1579-1591.
【発明の概要】
【発明が解決しようとする課題】
【００２５】
しかしながら、従来法では、音源に由来する直接音と反射音のみが収音されることが前提になっている。上記のほかに、収音信号にノイズが含まれている場合には、従来法で求めた逆フィルタは目的音以上にノイズを大幅に増幅してしまう。またノイズに由来する推定誤差のために、残響の除去も保証されない。そこで本発明では、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、ノイズの増幅を抑えることができる残響除去方法を提供することを目的とする。
【課題を解決するための手段】
【００２６】
本発明の残響除去方法は、多チャネル相関行列算出ステップと、信号部分空間抽出ステップと、逆フィルタ推定ステップと、逆フィルタステップとを有する。信号部分空間抽出ステップは、特異値分解サブステップと、信号部分空間次元推定サブステップとを有する。逆フィルタ推定ステップは、疑似逆行列推定サブステップと、逆フィルタ係数推定サブステップとを有する。
【００２７】
多チャネル相関行列算出ステップは、Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する。特異値分解サブステップは、算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する。信号部分空間次元推定サブステップは、出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する。疑似逆行列推定サブステップは、推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する。逆フィルタ係数推定サブステップは、推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する。逆フィルタステップは、推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、Ｎチャネルの収音信号のそれぞれに適用する。
【発明の効果】
【００２８】
本発明の残響除去方法によれば、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、ノイズの増幅を抑えることができる。
【図面の簡単な説明】
【００２９】
【図１】従来技術に係る残響除去装置の構成を示すブロック図。
【図２】従来技術に係る残響除去装置の動作を示すフローチャート。
【図３】実施例１に係る残響除去装置の構成を示すブロック図。
【図４】実施例１に係る残響除去装置の動作を示すフローチャート。
【図５】実施例２に係る残響除去装置の構成を示すブロック図。
【図６】実施例２に係る残響除去装置の動作を示すフローチャート。
【図７】実施例３に係る残響除去装置の構成を示すブロック図。
【図８】実施例３に係る残響除去装置の動作を示すフローチャート。
【図９】実施例４に係る残響除去装置の構成を示すブロック図。
【図１０】実施例４に係る残響除去装置の動作を示すフローチャート。
【図１１】特異値番号と特異値の関係を収音ＳＮ比毎に示す図。
【図１２】特異値の最大値から最小値に一定比率で減衰する関数ｐ（ｊ）を示す図。
【図１３】関数ｐ（ｊ）とｊ番目の特異値σｊの比で定義される関数ｑ（ｊ）を示す図。
【図１４】特異値番号と所定値以上の特異値逆数を一定値に置き換えた特異値逆数の関係を示す図。
【図１５】特異値の最大値から最小値に一定比率で減衰する関数ｐ〜（ｉ）を示す図。
【図１６】関数ｐ〜（ｉ）とｉ番目の特異値σｉの比で定義される関数ｑ〜（ｉ）を示す図。
【図１７】シュミレーションに用いたインパルス応答と残響曲線を示す図。
【図１８】従来技術の残響除去装置のシュミレーション結果を示す図。
【図１９】実施例１の残響除去装置のシュミレーション結果を示す図。
【図２０】実施例２の残響除去装置のシュミレーション結果を示す図。
【図２１】実施例３（共通多項式除算後）の残響除去装置のシュミレーション結果を示す図。
【図２２】実施例１（共通多項式除算前）の残響除去装置のシュミレーション結果を示す図。
【発明を実施するための形態】
【００３０】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【実施例１】
【００３１】
図３、図４を参照して本発明の最も基本的な構成である実施例１の残響除去装置について説明する。図３は本実施例に係る残響除去装置１００の構成を示すブロック図である。図４は本実施例に係る残響除去装置１００の動作を示すフローチャートである。本実施例の残響除去装置１００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部４０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部４０は、疑似逆行列推定手段４１と、逆フィルタ係数推定手段４２とを備える。
【００３２】
Ｎ個のマイクロホン１０−１〜１０−Ｎから収音される収音信号ｙ１（ｋ）〜ｙＮ（ｋ）には、目的音源１から放出された音およびノイズが含まれている。この収音信号ｙ１（ｋ）〜ｙＮ（ｋ）から以下のステップで残響除去信号を得る。多チャネル相関行列算出部２０は、式（６）（７）（８）をもちいてＮチャネル（Ｎは２以上の整数）の収音信号間ｙ１（ｋ）〜ｙＮ（ｋ）の相関関数ｒｉｊ（τ）から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列Ｒｉｊを生成し、当該相関行列ＲｉｊからＮＬ行ＮＬ列の多チャネル相関行列Ｒを算出する（Ｓ２０）。次に、信号部分空間抽出部３０において、ステップＳ２０で求めた多チャネル相関行列Ｒから、観測空間を信号部分空間とノイズ部分空間に分割する（Ｓ３０）。そして逆フィルタ推定部４０においてノイズ部分空間に含まれる信号成分のキャンセル、および信号部分空間に含まれる信号成分の残響除去を同時におこなう逆フィルタを求める（Ｓ４０）。具体的には、特異値分解手段３１は、算出された多チャネル相関行列Ｒを特異値分解して、ＮＬ個の特異値を出力する（ＳＳ３１）。より具体的には、特異値分解手段３１は、ＮＬ×ＮＬ多チャネル相関行列Ｒを式（１６）に示すように特異値分解する。行列Ｕは多チャネル相関行列Ｒの左特異値ベクトルを並べた行列である。σ１、…σＮＬは特異値を表す。ＶＴは、多チャネル相関行列Ｒの右特異値ベクトルを並べた行列の転置行列である。
【００３３】
【数１０】
【００３４】
次に、信号部分空間次元推定手段３２は、出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する（ＳＳ３２）。行列Ｖの各縦ベクトルは、観測空間の基底ベクトルになる。信号部分空間次元推定手段３２は、特異値σ１〜σＮＬのカーブから、信号部分空間を構成する基底ベクトルの本数Ｄ、すなわち信号部分空間の次元Ｄを推定する。推定された信号部分空間の次元Ｄをもちいて、観測空間を信号部分空間Ｖｓとノイズ部分空間Ｖｎに分割し、式（１７−１）のように信号部分空間成分のみからなる多チャネル相関行列Ｒ２を算出することができる。
【００３５】
【数１１】
【００３６】
ただし行列ＶＳのサイズはＮＬ×Ｄであり、その縦ベクトルが信号部分空間の基底ベクトルになる。また行列ＶＮの縦ベクトルがノイズ部分空間の基底ベクトルになる。行列ＶＳとＶＮはそれぞれ正規直交基底であり、行列ＶＳと行列ＶＮは直交する。
【００３７】
以下、信号部分空間次元推定手段３２が信号部分空間の次元Ｄを推定する方法を詳細に説明する。そのために、まず収音ＳＮ比の特異値への影響を示す。図１１に収音信号に付加する無相関白色雑音のレベルを上げていったときの、収音信号の多チャネル相関行列Ｒの特異値カーブを示す。図１１は特異値番号と特異値の関係を収音ＳＮ比毎に示す図である。多チャネル相関行列Ｒは、逆フィルタのフィルタ長を１チャネル当り１０００タップに設定して、４チャネル収音信号から求めた。４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測し、５５００タップに打ち切ったインパルス応答を白色雑音に畳み込んで生成した。この図からは、（１）ＳＮ比が低くノイズが多いほど、最小特異値は大きくなること、（２）ノイズが多いほど、特異値の減少は早く飽和すること、（３）特異値のカーブは、信号に由来する急峻な減衰カーブＡとノイズに由来する比較的フラットなカーブＢからなること、（４）カーブＡとカーブＢが交錯するところで、特異値カーブの傾きが大きく変化していること、以上（１）〜（４）の特徴を読み取ることができる。
【００３８】
以下、図１２、１３を参照して信号部分空間の次元Ｄを求める方法を記載する。図１２は特異値の最大値から最小値に一定比率で減衰する関数ｐ（ｊ）を示す図である。図１３は関数ｐ（ｊ）とｊ番目の特異値σｊの比で定義される関数ｑ（ｊ）を示す図である。まず、σ１からσＮＬに一定比率で減衰するカーブｐ（ｊ）を求める（図１２参照）。
【００３９】
【数１２】
【００４０】
ｐ（ｊ）とσｊの比ｑ（ｊ）をとり、その最大値の位置から変曲点位置Ｊｏを求める（図１３参照）。
【００４１】
【数１３】
【００４２】
変曲点での値ｑ（Ｊｏ）よりβ［ｄＢ］小さい値をとる点Ｊ−、Ｊ＋（Ｊ−＜Ｊ＋）求め、信号部分空間の次元をＤ＝Ｊ＋とする。βは例えば０〜８［ｄＢ］の範囲で指定する。以上のステップにより信号部分空間の抽出が可能である。
【００４３】
次に、疑似逆行列推定手段４１は、推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列Ｒ２の疑似逆行列Ｒ２＃を推定する（ＳＳ４１）。より具体的には、疑似逆行列推定手段４１は、式（２０）をもちいて、疑似逆行列すなわち信号部分空間の逆行列Ｒ２＃を求める。
【００４４】
【数１４】
【００４５】
逆フィルタ係数推定手段４２は、推定された疑似逆行列Ｒ２＃に基づいて逆フィルタ係数ベクトルを推定する（ＳＳ４２）。より具体的には、逆フィルタ係数推定手段４２は、サブステップＳＳ４１で求めた疑似逆行列Ｒ２＃をもちい、式（２１）により逆フィルタ係数ベクトルを求める。
【００４６】
【数１５】
【００４７】
逆フィルタ部７０は、推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、Ｎチャネルの収音信号のそれぞれに適用する（Ｓ７０）。より具体的には、逆フィルタ部７０では、ステップＳ４０（ＳＳ４１、ＳＳ４２）で求めた逆フィルタ係数ベクトルｇを式（１２）（１３）により各チャネルのフィルタ係数ベクトルｇ１〜ｇＮに分解し、式（１４）により、Ｎチャネル収音信号に適用する。上記で算出した逆フィルタの効果を見るために、逆フィルタ後の信号の特性を調べる。Ｌサンプル分のＮチャネル収音信号ベクトルｙは、信号部分空間とノイズ部分空間に射影することで、式（２２）のように信号部分空間の成分とノイズ部分空間の成分に分解できる。
【００４８】
【数１６】
【００４９】
これに逆フィルタを適用すると行列ＶＳとＶＮが直交しているために、式（２３）のように信号部分空間の成分のみが残り、ノイズ部分空間の成分はキャンセルされる。
【００５０】
【数１７】
【００５１】
本実施例の残響除去装置１００によれば、信号部分空間抽出部３０が特異値のカーブの減衰傾向から信号部分空間の次元Ｄを推定して観測空間を信号部分空間とノイズ部分空間に分割し、逆フィルタ推定部４０が、推定された次元Ｄに基づいて推定した逆フィルタにより、信号部分空間に含まれる信号成分の残響を除去すると同時に、ノイズ部分空間に含まれる信号成分をキャンセルすることができるため、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、かつノイズの増幅を抑えることができる。
【実施例２】
【００５２】
次に、図５、図６を参照して実施例１の残響除去装置における疑似逆行列推定手段に変更を加えた実施例２の残響除去装置について説明する。図５は本実施例に係る残響除去装置２００の構成を示すブロック図である。図６は本実施例に係る残響除去装置２００の動作を示すフローチャートである。本実施例の残響除去装置２００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部５０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部５０は、疑似逆行列推定手段５１と、逆フィルタ係数推定手段４２とを備える。本実施例の残響除去装置２００と、実施例１の残響除去装置１００の違いは、実施例１の残響除去装置１００の疑似逆行列推定手段４１が本実施例の残響除去装置２００において疑似逆行列推定手段５１に置き換えられている点、および当該疑似逆行列推定手段５１を有する逆フィルタ推定部５０について実施例１の残響除去装置１００において対応する逆フィルタ推定部４０と番号を変えて表示した点のみである。従って、実施例１の残響除去装置１００と共通する番号を有する各構成部は、実施例１の残響除去装置１００において同一番号を有する各構成部と全く同じ機能を有するため説明を割愛する。
【００５３】
以下、実施例１との相違点であるサブステップＳＳ５１について説明する。疑似逆行列推定手段５１は、信号部分空間成分のみからなる多チャネル相関行列Ｒ２の疑似逆行列Ｒ♯２の対角行列の要素である特異値の逆数１／σ１、１／σ２、…、１／σＤのうち、予め定めたしきい値以上となる要素１／σＤ−、…、１／σＤを１／σ'Ｚ＜１／σＺ（Ｄ−≦Ｚ≦Ｄ）を充たす要素１／σ'Ｄ−、…、１／σ'Ｄに置き換えて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する（ＳＳ５１）。
【００５４】
参考非特許文献１（参考非特許文献１：浅野太著、「（音響テクノロジーシリーズ１６）音のアレイ信号処理−音源の定位・追跡と分離−」、コロナ社、2011年1月27日、p.111-p.113）によれば、特異値は対応する特異値ベクトルの張る空間における信号パワーと雑音パワーの和である。信号パワーは偏在する傾向があり、雑音パワーは均一に分布する傾向がある。そのため、小さい特異値に対応する特異ベクトルの張る空間では、雑音パワーの占める比率がより大きくなる。このことを念頭に、疑似逆行列推定手段５１は、上記実施例で使用している疑似逆行列
【００５５】
【数１８】
【００５６】
を変形する。式（２４）によれば、疑似逆行列の対角項は特異値の逆数になっている。このため、疑似逆行列は雑音パワー比の大きい成分に、より大きい値を乗算し、ノイズを増幅しやすいことが分かる。そこで、雑音パワー比の大きい成分に特異値の逆数よりも小さい値を乗算することで、ノイズ成分の増幅を抑える。そのために、疑似逆行列推定手段５１は、疑似逆行列Ｒ２＃の代わりに行列Ｑを用いる。
行列Ｑの対角項は、
【００５７】
【数１９】
【００５８】
に設定されている。これにより、残響除去の性能をほぼ維持しつつ、逆フィルタによるノイズ増幅を抑えることができる。行列Ｑの一具体形として、式（２５）のように１／σＤ−をもちいる方法が考えられる。
【００５９】
【数２０】
【００６０】
なおＤ−は、Ｊ−〜Ｊ＋の範囲から指定する。図１４に、特異値、特異値の逆数、行列Ｑの対角項の関係を示す。図１４（ａ）は特異値番号と特異値の関係を示す図である。図１４（ｂ）は特異値番号と所定値以上の特異値逆数を一定値に置き換えた特異値逆数（行列Ｑの対角項）の関係を示す図である。
【００６１】
以下、実施例１の残響除去装置１００、実施例２の残響除去装置２００について、シミュレーションにより検証した効果について図１７、図１８、図１９、図２０を参照して説明する。図１７はシュミレーションに用いたインパルス応答（ａ）と残響曲線（ｂ）を示す図である。図１８は従来技術の残響除去装置のシュミレーション結果を示す図である。図１９は実施例１の残響除去装置１００のシュミレーション結果を示す図である。図２０は実施例２の残響除去装置２００のシュミレーション結果を示す図である。シュミレーションに用いた音源は白色雑音であり、４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測したインパルス応答を５５００タップに打ち切って使用した。各マイクロホンには無相関の白色ノイズをＳＮ比２５ｄＢになるよう加えている。逆フィルタのフィルタ長は１チャネル当り１０００タップに設定した。相関関数を算出するときには、２０ｓの収音信号を使用した。図１８、図１９、図２０はいずれも（ａ）音響経路＋逆フィルタのインパルス応答、（ｂ）収音経路＋逆フィルタのインパルス応答の残響曲線（実線）、収音経路インパルス応答の残響曲線（点線）、（ｃ）逆フィルタ後の信号、（ｄ）逆フィルタ後のノイズである。なお、図１８〜図２２のグラフ（ａ）は縦軸を信号強度、横軸をタップ、グラフ（ｂ）は縦軸をエネルギー、横軸を時間、グラフ（ｃ）は縦軸を信号強度、横軸を時間、グラフ（ｄ）は縦軸を信号強度、横軸を時間、としてそれぞれプロットしたものである。図１８の従来法の処理結果では、逆フィルタをうまく求められていないため、音響経路＋逆フィルタのインパルス応答は、ピーク位置が拡散している。その残響曲線もほとんど減衰していない。また逆フィルタ後は、ノイズの振幅が信号の４倍にも達している。図１９の実施例１の処理結果では、逆フィルタがうまく求められており、音響経路＋逆フィルタのインパルス応答は、ピークが１つ立っている。その残響曲線も、８０サンプル目のところで−１０ｄＢ減衰している。そして逆フィルタ後は、ノイズの振幅が信号振幅の半分程度に抑えられている。図２０の実施例２の処理結果では、より良い逆フィルタが求められており、音響経路＋逆フィルタインパルス応答の残響曲線の減衰は、実施例１の処理結果と比較して２ｄＢ改善している。そして逆フィルタ後は、ノイズの振幅が信号振幅の１／３程度に抑えられている。
【００６２】
本実施例の残響除去装置２００によれば、実施例１の残響除去装置１００の効果に加え、疑似逆行列推定手段５１が雑音パワー比の大きい成分に乗算される特異値の逆数を当該特異値の逆数よりも小さな値に置き換えて疑似逆行列を推定するため、残響除去の性能をほぼ維持しつつ、逆フィルタによるノイズ増幅を抑えることができる。
【実施例３】
【００６３】
次に、図７、図８、図１５、図１６を参照して実施例１の残響除去装置１００に共通多項式を除去する機能を追加した実施例３に係る残響除去装置について説明する。図７は本実施例に係る残響除去装置３００の構成を示すブロック図である。図８は本実施例に係る残響除去装置３００の動作を示すフローチャートである。図１５は特異値の最大値から最小値に一定比率で減衰する関数ｐ〜（ｉ）を示す図である。図１６は関数ｐ〜（ｉ）とｉ番目の特異値σｉの比で定義される関数ｑ〜（ｉ）を示す図である。本実施例の残響除去装置３００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部４０と、共通多項式除去部６０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部４０は、疑似逆行列推定手段４１と、逆フィルタ係数推定手段４２とを備える。共通多項式除去部６０は、シルベスタ行列生成手段６１と、左ゼロ空間抽出手段６２と、共通多項式次元推定手段６３と、共通多項式フィルタ係数推定手段６４と、共通多項式除算手段６５とを備える。本実施例の残響除去装置３００と、実施例１の残響除去装置１００の違いは、実施例１の残響除去装置１００が備えない共通多項式除去部６０を、本実施例の残響除去装置３００が備える点のみである。従って、実施例１の残響除去装置１００と共通する番号を有する各構成部は、実施例１の残響除去装置１００において同一番号を有する各構成部と全く同じ機能を有するため説明を割愛する。
【００６４】
実施例１の残響除去装置１００は、音源信号が白色であること、すなわち音源信号のスペクトルがフラットなことを前提としている。音源信号が音声のように有色な場合には、実施例１の逆フィルタは音響経路に由来する残響を除去すると同時に、音源信号の周波数特性の平坦化もはかってしまう。有色音源信号に対応するには、逆フィルタｇ１（ｋ’）〜ｇＮ（ｋ’）で音源周波数特性を平坦化している部分、すなわち共通多項式を求めて除去すればよい。以下、共通多項式の次元を推定し、この推定次元をもとに参考非特許文献２のＱｉｕの手法で共通多項式を求め、第一の逆フィルタから共通多項式を除算する処理を説明する（参考非特許文献２：W. Qiu, Y. Hua, and K. Abed-Meraim, " Subspace Method for the Computation of the GCD of Polynomials" Automatica,1997,vol. 33, no. 4, pp. 741-743）。なお次元が既知のとき共通多項式を求めるアルゴリズムは、参考非特許文献２以外に参考非特許文献３などの手法があり、これらを代わりに使用してもよい（参考非特許文献３：R. M. Corless, S. M. Watt, and L. Zhi, QR "Factoring to Compute the GCD of Univariate Approximate Polynomials", IEEE Transactions on Signal Processing, 2004, vol. 52, no. 12, pp. 3394-3402）。
【００６５】
まず、シルベスタ行列生成手段６１は、逆フィルタ係数推定手段４２において推定された逆フィルタ係数ベクトル（以下、第１の逆フィルタ係数ベクトルという）から（２Ｌ−１）行ＮＬ列のシルベスタ行列を生成する（ＳＳ６１）。より具体的には、シルベスタ行列生成手段６１は、逆フィルタ係数推定手段４２で求められた第一の逆フィルタ係数ベクトルから、式（２６）（２７）で（２Ｌ−１）×ＮＬシルベスタ行列を生成する。
【００６６】
【数２１】
【００６７】
左ゼロ空間抽出手段６２は、生成されたシルベスタ行列を特異値分解して（２Ｌ−１）個の特異値を出力する（ＳＳ６２）。共通多項式の次元がＤ〜のとき、特異値分解は式（２８）のように分割される（非特許文献６）。
【００６８】
【数２２】
【００６９】
行列Ｕ〜０の縦ベクトルが行列Ｇの左ゼロ空間の基底ベクトルになり、その本数がＤ〜になる。特異値カーブは図１５のようになる。共通多項式次元推定手段６３は、左ゼロ空間抽出手段６２が出力した特異値の減衰傾向に基づいて共通多項式の次元を推定する（ＳＳ６３）。より具体的には、共通多項式次元推定手段６３は、行列Ｇの特異値σ〜１からσ〜２Ｌ−１に一定比率で減衰するカーブｐ〜（ｉ）を求める（図１５参照）。
【００７０】
【数２３】
【００７１】
次に、ｐ〜（ｉ）とσ〜ｉの比ｑ〜（ｉ）をとり、その最小値の位置から変曲点位置Ｉｏを求める（図１６参照）。
【００７２】
【数２４】
【００７３】
変曲点での値ｑ〜（Ｉ０）よりγ［ｄＢ］小さい値をとる二点のうち大きいほうの点Ｉ＋を求める（図１６参照）。そして共通多項式の次元をＤ〜＝２Ｌ−１−Ｉ＋とする。なおγは３〜２０ｄＢの範囲で指定する。共通多項式フィルタ係数推定手段６４は、シルベスタ行列の左ゼロ空間の基底ベクトルから共通多項式のフィルタ係数を推定する（ＳＳ６４）。より具体的には、共通多項式フィルタ係数推定手段６４は、左ゼロ空間の基底ベクトルから、式（３１−１〜３１−３）で行列Ｗを算出する。行列Ｗの最小固有値に対応する固有ベクトルが共通多項式のフィルタ係数ｃとなる。
【００７４】
【数２５】
【００７５】
共通多項式除算手段６５は、第１の逆フィルタ係数ベクトルをチャネル毎の共通多項式で除算して第２の逆フィルタ係数ベクトルを推定する（ＳＳ６５）。より具体的には、共通多項式除算手段６５は、逆フィルタ推定手段４２で求めた第一の逆フィルタを、チャネルごとに共通多項式で除算する。そのために式（３２）を解いて、チャネルごとに第２の逆フィルタ係数
【００７６】
【数２６】
【００７７】
を求める。
【００７８】
【数２７】
【００７９】
逆フィルタ部７０では、式（３３）により収音信号ｙ１（ｋ）〜ｙＮ（ｋ）に逆フィルタ係数を適用して出力信号ｘ＾（ｋ）を得る。
【００８０】
【数２８】
【００８１】
実施例３の残響除去装置３００について、シミュレーションにより効果を検証した。音源は音声であり、４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測したインパルス応答を５５００タップに打ち切って使用した。各マイクロホンには無相関の白色ノイズをＳＮ比２０ｄＢになるよう加えている。逆フィルタのフィルタ長は１チャネル当り１０００タップに設定した。相関関数を算出するときには、２０ｓの収音信号を使用した。実施例２の手法を適用した結果を図２１に、実施例１の手法を適用した結果を図２２に示す。いずれも（ａ）音響経路＋逆フィルタのインパルス応答、（ｂ）収音経路＋逆フィルタのインパルス応答の残響曲線（実線）、収音経路インパルス応答の残響曲線（点線）、（ｃ）逆フィルタ後の信号、（ｄ）逆フィルタ後のノイズである。
【００８２】
実施例１の手法を適用した場合には、図２２（ａ）インパルス応答ではきれいなピーク形状が得られず、ピーク後部に成分が残留している。また図２２（ｂ）の残響曲線でも８０サンプル付近からなだらかに減衰しており、残響がうまく取れていないことが分かる。
【００８３】
しかし実施例３の手法を適用することで、図２１（ａ）の音響経路＋逆フィルタのインパルス応答ではきれいなピークが得られており、その残響曲線も８０サンプル付近で０ｄＢから−１０ｄＢに急激に減衰している、このことから、音源が音声信号であっても実施例３の手法により残響を効果的に除去できることが分かる。
【００８４】
本実施例の残響除去装置３００によれば、実施例１の残響除去装置１００の効果に加えて、音源信号が有色である場合であっても、共通多項式除去部６０が共通多項式を求めて除去することにより、有色音源信号にも対応することができる。
【実施例４】
【００８５】
以下、図９、図１０を参照して実施例２の残響除去装置２００に実施例３の残響除去装置３００が備える共通多項式除去部６０を加えた実施例４の残響除去装置について説明する。図９は本実施例に係る残響除去装置４００の構成を示すブロック図である。図１０は本実施例に係る残響除去装置４００の動作を示すフローチャートである。本実施例の残響除去装置４００は、実施例２の残響除去装置２００と共通してＮチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部５０と、逆フィルタ部７０とを備える。また、本実施例の残響除去装置４００は、実施例３の残響除去装置３００と共通して共通多項式除去部６０を備える。従って、実施例２の残響除去装置２００と共通する番号を有する各構成部の機能は、実施例２の記載を参照されたい。また、実施例３の残響除去装置３００と共通する番号を有する各構成部の機能は、実施例３の記載を参照されたい。
【００８６】
また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００８７】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
【００８８】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【００８９】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００９０】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００９１】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【技術分野】
【０００１】
本発明は残響とノイズを含む多チャネル収音信号から音源の発する本来の信号を抽出する残響除去方法、残響除去装置、プログラムに関する。
【背景技術】
【０００２】
マイクロホンで音声等を収音する場合、マイクロホンに到達する音は、音源から発せられて直接到達した直接音と、壁等に１回反射して到達した反射音と、複数回反射して到達した反射音の混ざったものになる。マイクロホンが音源に近接していない場合は、この反射音からなる残響成分が多くなる。話者音声を離れた位置で収音する場合や残響時間が長い部屋で収音する場合には、残響により音声の明瞭度が低下して聞き取り難くなり、例えば音声認識システムの音声認識率も著しく低下することが知られている。話者位置からマイクロホン位置までのインパルス応答が予め測定されている場合には、ＭＩＮＴ（Ｍｕｌｔｉ−ｉｎｐｕｔ／ｏｕｔｐｕｔ Ｉｎｖｅｒｓｅ Ｔｈｅｏｒｅｍ）法（非特許文献１）を用いてインパルス応答から逆フィルタ係数を求め、収音信号に逆フィルタを適用することで残響に含まれる収音信号から原音を回復することができる。
【０００３】
ＭＩＮＴ法を簡単に説明する。マイクロホンがＮ個あり、目的音源からｉ番目のマイクロホンまでの音響経路のインパルス応答をｈｉ（０）〜ｈｉ（Ｋ）とする。Ｌタップの逆フィルタを求めるとき、ＭＩＮＴ法は目的インパルス応答ベクトルをｂ、ｉ番目のマイクロホン収音信号に適用するＬ次元の逆フィルタ係数ベクトルをｇｎ（ｎ＝１〜Ｎ）としてサイズ（Ｌ＋Ｋ）×Ｌのインパルス応答畳み込み行列（式（１））を用いて、式（２）若しくは式（３）の連立方程式を解いて逆フィルタ係数ベクトルｇｎを求める。
【０００４】
【数１】
【０００５】
完全な逆フィルタを求めるには、Ｌ＋Ｋ≧ＮＬが満たされている必要がある。Ｌ＋Ｋ<ＮＬの関係を満たすＬをもちいた場合、近似の逆フィルタが求まる。Ｌを小さくしていくと、その近似精度は下がってく。目的インパルス応答ベクトルｂはＬ＋Ｋ次元のベクトルで、無遅延のパルス（式（４））や遅延付きのパルスが（式（５））用いられる。
【０００６】
【数２】
【０００７】
通常は、音源位置（話者位置）からマイクロホン位置までのインパルス応答を予め測定することが出来ないため、残響に含まれる収音信号のみから原音を回復するブラインド残響除去法が用いられる。ブラインド残響除去の一方法としてｂｌｉｎｄ−ＭＩＮＴ法（特許文献１、非特許文献２）がある。
【０００８】
図１、図２を参照してｂｌｉｎｄ−ＭＩＮＴ法を簡単に説明する。図１は従来技術（ｂｌｉｎｄ−ＭＩＮＴ法）に係る残響除去装置９００の構成を示すブロック図である。図２は従来技術（ｂｌｉｎｄ−ＭＩＮＴ法）に係る残響除去装置９００の動作を示すフローチャートである。残響除去装置９００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、逆フィルタ推定部９０と、逆フィルタ部７０とを備える。
【０００９】
目的音源１はスペクトルがフラットな音を発し、Ｎ個のマイクロホン１０−１〜１０−Ｎで収音される。収音信号ｙ１（ｋ）〜ｙＮ（ｋ）は、多チャネル相関行列算出部２０に入力される。ここで、収音信号ｙ１（ｋ）〜ｙＮ（ｋ）は、例えば標本化周波数は１６ｋＨｚで離散値化されたディジタル信号であり、ｋは時刻を表す。なお、収音信号をディジタル化するＡ／Ｄ変換器等については省略している。多チャネル相関行列算出部２０は、収音信号ｙ１（ｋ）〜ｙＮ（ｋ）間の相関関数ｒｉｊ（τ）を式（６）で計算する。
【００１０】
【数３】
【００１１】
ここでτは、２時点の時間間隔である。多チャネル相関行列算出部２０は、相関関数ｒｉｊ（τ）からＬ×Ｌ相関行列Ｒｉｊ（式（７））を生成し、ＮＬ×ＮＬ多チャネル相関行列Ｒ（式（８））を生成する（Ｓ２０）。
【００１２】
【数４】
【００１３】
但し、Ｌは逆フィルタの１チャネル当りのタップ数である。逆フィルタ推定部９０は、式（９）の連立一次方程式を解いて逆フィルタ係数ベクトルｇを推定する（Ｓ９０）。
【００１４】
【数５】
【００１５】
ここで、目的インパルス応答ベクトルｂ（式（１０））はＬＮ次元ベクトルであり、Ｌ次元縦ベクトルｂ１〜ｂＮが縦に並んでいる（式（１１））。
【００１６】
【数６】
【００１７】
Ｌ次元ベクトルｂｎ中のδｎは、第ｎ番目のマイクロホン１０−ｎと目的音源との位置関係で決まる。第ｎ番目のマイクロホン１０−ｎが全マイクロホンの中で目的音源１に最も近いときにδｎは１を取り、それ以外のときは０を取る。逆フィルタ係数ベクトルｇは、ｎ番目のマイクロホンの収音信号ｙｎ（ｋ）に適用するＬ次元の逆フィルタ係数ベクトル（式（１２））を用いると、式（１３）と書ける。
【００１８】
【数７】
【００１９】
逆フィルタ部７０は、式（１４）に示すように各収音信号ｙ１（ｋ）〜ｙＮ（ｋ）に逆フィルタ係数ｇを適用して最終出力信号ｘ＾（ｋ）を得る。この出力信号は、残響の除去された目的信号になる。
【００２０】
【数８】
【００２１】
従来のブラインド残響除去装置９００は、目的音源１から各マイクロホンまでの未知音響経路インパルス応答ｈｉ（０）〜ｈｉ（Ｌ−１）から成る行列Ｈと収音信号から求めた相関行列Ｒとの間に式（１５）の関係があることを利用し、逆フィルタ係数を求めている。αは比例定数（スカラー量）である。
【００２２】
【数９】
【先行技術文献】
【特許文献】
【００２３】
【特許文献１】特開２００６−６６９８９号公報
【非特許文献】
【００２４】
【非特許文献１】大賀寿郎、山崎芳男、金田豊著、「音響システムとディジタル処理」、電子情報通信学会、1995年4月
【非特許文献２】K. Furuya and A. Kataoka, "Robust Speech Dereverberatio Using Multichannel Blind Deconvolution With Spectral Subtraction", IEEE Transactions on Audio, Speech, and Language Processing, 2007 July, vol. 15, no. 5, pp. 1579-1591.
【発明の概要】
【発明が解決しようとする課題】
【００２５】
しかしながら、従来法では、音源に由来する直接音と反射音のみが収音されることが前提になっている。上記のほかに、収音信号にノイズが含まれている場合には、従来法で求めた逆フィルタは目的音以上にノイズを大幅に増幅してしまう。またノイズに由来する推定誤差のために、残響の除去も保証されない。そこで本発明では、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、ノイズの増幅を抑えることができる残響除去方法を提供することを目的とする。
【課題を解決するための手段】
【００２６】
本発明の残響除去方法は、多チャネル相関行列算出ステップと、信号部分空間抽出ステップと、逆フィルタ推定ステップと、逆フィルタステップとを有する。信号部分空間抽出ステップは、特異値分解サブステップと、信号部分空間次元推定サブステップとを有する。逆フィルタ推定ステップは、疑似逆行列推定サブステップと、逆フィルタ係数推定サブステップとを有する。
【００２７】
多チャネル相関行列算出ステップは、Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する。特異値分解サブステップは、算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する。信号部分空間次元推定サブステップは、出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する。疑似逆行列推定サブステップは、推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する。逆フィルタ係数推定サブステップは、推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する。逆フィルタステップは、推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、Ｎチャネルの収音信号のそれぞれに適用する。
【発明の効果】
【００２８】
本発明の残響除去方法によれば、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、ノイズの増幅を抑えることができる。
【図面の簡単な説明】
【００２９】
【図１】従来技術に係る残響除去装置の構成を示すブロック図。
【図２】従来技術に係る残響除去装置の動作を示すフローチャート。
【図３】実施例１に係る残響除去装置の構成を示すブロック図。
【図４】実施例１に係る残響除去装置の動作を示すフローチャート。
【図５】実施例２に係る残響除去装置の構成を示すブロック図。
【図６】実施例２に係る残響除去装置の動作を示すフローチャート。
【図７】実施例３に係る残響除去装置の構成を示すブロック図。
【図８】実施例３に係る残響除去装置の動作を示すフローチャート。
【図９】実施例４に係る残響除去装置の構成を示すブロック図。
【図１０】実施例４に係る残響除去装置の動作を示すフローチャート。
【図１１】特異値番号と特異値の関係を収音ＳＮ比毎に示す図。
【図１２】特異値の最大値から最小値に一定比率で減衰する関数ｐ（ｊ）を示す図。
【図１３】関数ｐ（ｊ）とｊ番目の特異値σｊの比で定義される関数ｑ（ｊ）を示す図。
【図１４】特異値番号と所定値以上の特異値逆数を一定値に置き換えた特異値逆数の関係を示す図。
【図１５】特異値の最大値から最小値に一定比率で減衰する関数ｐ〜（ｉ）を示す図。
【図１６】関数ｐ〜（ｉ）とｉ番目の特異値σｉの比で定義される関数ｑ〜（ｉ）を示す図。
【図１７】シュミレーションに用いたインパルス応答と残響曲線を示す図。
【図１８】従来技術の残響除去装置のシュミレーション結果を示す図。
【図１９】実施例１の残響除去装置のシュミレーション結果を示す図。
【図２０】実施例２の残響除去装置のシュミレーション結果を示す図。
【図２１】実施例３（共通多項式除算後）の残響除去装置のシュミレーション結果を示す図。
【図２２】実施例１（共通多項式除算前）の残響除去装置のシュミレーション結果を示す図。
【発明を実施するための形態】
【００３０】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【実施例１】
【００３１】
図３、図４を参照して本発明の最も基本的な構成である実施例１の残響除去装置について説明する。図３は本実施例に係る残響除去装置１００の構成を示すブロック図である。図４は本実施例に係る残響除去装置１００の動作を示すフローチャートである。本実施例の残響除去装置１００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部４０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部４０は、疑似逆行列推定手段４１と、逆フィルタ係数推定手段４２とを備える。
【００３２】
Ｎ個のマイクロホン１０−１〜１０−Ｎから収音される収音信号ｙ１（ｋ）〜ｙＮ（ｋ）には、目的音源１から放出された音およびノイズが含まれている。この収音信号ｙ１（ｋ）〜ｙＮ（ｋ）から以下のステップで残響除去信号を得る。多チャネル相関行列算出部２０は、式（６）（７）（８）をもちいてＮチャネル（Ｎは２以上の整数）の収音信号間ｙ１（ｋ）〜ｙＮ（ｋ）の相関関数ｒｉｊ（τ）から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列Ｒｉｊを生成し、当該相関行列ＲｉｊからＮＬ行ＮＬ列の多チャネル相関行列Ｒを算出する（Ｓ２０）。次に、信号部分空間抽出部３０において、ステップＳ２０で求めた多チャネル相関行列Ｒから、観測空間を信号部分空間とノイズ部分空間に分割する（Ｓ３０）。そして逆フィルタ推定部４０においてノイズ部分空間に含まれる信号成分のキャンセル、および信号部分空間に含まれる信号成分の残響除去を同時におこなう逆フィルタを求める（Ｓ４０）。具体的には、特異値分解手段３１は、算出された多チャネル相関行列Ｒを特異値分解して、ＮＬ個の特異値を出力する（ＳＳ３１）。より具体的には、特異値分解手段３１は、ＮＬ×ＮＬ多チャネル相関行列Ｒを式（１６）に示すように特異値分解する。行列Ｕは多チャネル相関行列Ｒの左特異値ベクトルを並べた行列である。σ１、…σＮＬは特異値を表す。ＶＴは、多チャネル相関行列Ｒの右特異値ベクトルを並べた行列の転置行列である。
【００３３】
【数１０】
【００３４】
次に、信号部分空間次元推定手段３２は、出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する（ＳＳ３２）。行列Ｖの各縦ベクトルは、観測空間の基底ベクトルになる。信号部分空間次元推定手段３２は、特異値σ１〜σＮＬのカーブから、信号部分空間を構成する基底ベクトルの本数Ｄ、すなわち信号部分空間の次元Ｄを推定する。推定された信号部分空間の次元Ｄをもちいて、観測空間を信号部分空間Ｖｓとノイズ部分空間Ｖｎに分割し、式（１７−１）のように信号部分空間成分のみからなる多チャネル相関行列Ｒ２を算出することができる。
【００３５】
【数１１】
【００３６】
ただし行列ＶＳのサイズはＮＬ×Ｄであり、その縦ベクトルが信号部分空間の基底ベクトルになる。また行列ＶＮの縦ベクトルがノイズ部分空間の基底ベクトルになる。行列ＶＳとＶＮはそれぞれ正規直交基底であり、行列ＶＳと行列ＶＮは直交する。
【００３７】
以下、信号部分空間次元推定手段３２が信号部分空間の次元Ｄを推定する方法を詳細に説明する。そのために、まず収音ＳＮ比の特異値への影響を示す。図１１に収音信号に付加する無相関白色雑音のレベルを上げていったときの、収音信号の多チャネル相関行列Ｒの特異値カーブを示す。図１１は特異値番号と特異値の関係を収音ＳＮ比毎に示す図である。多チャネル相関行列Ｒは、逆フィルタのフィルタ長を１チャネル当り１０００タップに設定して、４チャネル収音信号から求めた。４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測し、５５００タップに打ち切ったインパルス応答を白色雑音に畳み込んで生成した。この図からは、（１）ＳＮ比が低くノイズが多いほど、最小特異値は大きくなること、（２）ノイズが多いほど、特異値の減少は早く飽和すること、（３）特異値のカーブは、信号に由来する急峻な減衰カーブＡとノイズに由来する比較的フラットなカーブＢからなること、（４）カーブＡとカーブＢが交錯するところで、特異値カーブの傾きが大きく変化していること、以上（１）〜（４）の特徴を読み取ることができる。
【００３８】
以下、図１２、１３を参照して信号部分空間の次元Ｄを求める方法を記載する。図１２は特異値の最大値から最小値に一定比率で減衰する関数ｐ（ｊ）を示す図である。図１３は関数ｐ（ｊ）とｊ番目の特異値σｊの比で定義される関数ｑ（ｊ）を示す図である。まず、σ１からσＮＬに一定比率で減衰するカーブｐ（ｊ）を求める（図１２参照）。
【００３９】
【数１２】
【００４０】
ｐ（ｊ）とσｊの比ｑ（ｊ）をとり、その最大値の位置から変曲点位置Ｊｏを求める（図１３参照）。
【００４１】
【数１３】
【００４２】
変曲点での値ｑ（Ｊｏ）よりβ［ｄＢ］小さい値をとる点Ｊ−、Ｊ＋（Ｊ−＜Ｊ＋）求め、信号部分空間の次元をＤ＝Ｊ＋とする。βは例えば０〜８［ｄＢ］の範囲で指定する。以上のステップにより信号部分空間の抽出が可能である。
【００４３】
次に、疑似逆行列推定手段４１は、推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列Ｒ２の疑似逆行列Ｒ２＃を推定する（ＳＳ４１）。より具体的には、疑似逆行列推定手段４１は、式（２０）をもちいて、疑似逆行列すなわち信号部分空間の逆行列Ｒ２＃を求める。
【００４４】
【数１４】
【００４５】
逆フィルタ係数推定手段４２は、推定された疑似逆行列Ｒ２＃に基づいて逆フィルタ係数ベクトルを推定する（ＳＳ４２）。より具体的には、逆フィルタ係数推定手段４２は、サブステップＳＳ４１で求めた疑似逆行列Ｒ２＃をもちい、式（２１）により逆フィルタ係数ベクトルを求める。
【００４６】
【数１５】
【００４７】
逆フィルタ部７０は、推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、Ｎチャネルの収音信号のそれぞれに適用する（Ｓ７０）。より具体的には、逆フィルタ部７０では、ステップＳ４０（ＳＳ４１、ＳＳ４２）で求めた逆フィルタ係数ベクトルｇを式（１２）（１３）により各チャネルのフィルタ係数ベクトルｇ１〜ｇＮに分解し、式（１４）により、Ｎチャネル収音信号に適用する。上記で算出した逆フィルタの効果を見るために、逆フィルタ後の信号の特性を調べる。Ｌサンプル分のＮチャネル収音信号ベクトルｙは、信号部分空間とノイズ部分空間に射影することで、式（２２）のように信号部分空間の成分とノイズ部分空間の成分に分解できる。
【００４８】
【数１６】
【００４９】
これに逆フィルタを適用すると行列ＶＳとＶＮが直交しているために、式（２３）のように信号部分空間の成分のみが残り、ノイズ部分空間の成分はキャンセルされる。
【００５０】
【数１７】
【００５１】
本実施例の残響除去装置１００によれば、信号部分空間抽出部３０が特異値のカーブの減衰傾向から信号部分空間の次元Ｄを推定して観測空間を信号部分空間とノイズ部分空間に分割し、逆フィルタ推定部４０が、推定された次元Ｄに基づいて推定した逆フィルタにより、信号部分空間に含まれる信号成分の残響を除去すると同時に、ノイズ部分空間に含まれる信号成分をキャンセルすることができるため、収音信号にノイズが含まれている場合でも、ロバストにブラインド残響し、かつノイズの増幅を抑えることができる。
【実施例２】
【００５２】
次に、図５、図６を参照して実施例１の残響除去装置における疑似逆行列推定手段に変更を加えた実施例２の残響除去装置について説明する。図５は本実施例に係る残響除去装置２００の構成を示すブロック図である。図６は本実施例に係る残響除去装置２００の動作を示すフローチャートである。本実施例の残響除去装置２００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部５０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部５０は、疑似逆行列推定手段５１と、逆フィルタ係数推定手段４２とを備える。本実施例の残響除去装置２００と、実施例１の残響除去装置１００の違いは、実施例１の残響除去装置１００の疑似逆行列推定手段４１が本実施例の残響除去装置２００において疑似逆行列推定手段５１に置き換えられている点、および当該疑似逆行列推定手段５１を有する逆フィルタ推定部５０について実施例１の残響除去装置１００において対応する逆フィルタ推定部４０と番号を変えて表示した点のみである。従って、実施例１の残響除去装置１００と共通する番号を有する各構成部は、実施例１の残響除去装置１００において同一番号を有する各構成部と全く同じ機能を有するため説明を割愛する。
【００５３】
以下、実施例１との相違点であるサブステップＳＳ５１について説明する。疑似逆行列推定手段５１は、信号部分空間成分のみからなる多チャネル相関行列Ｒ２の疑似逆行列Ｒ♯２の対角行列の要素である特異値の逆数１／σ１、１／σ２、…、１／σＤのうち、予め定めたしきい値以上となる要素１／σＤ−、…、１／σＤを１／σ'Ｚ＜１／σＺ（Ｄ−≦Ｚ≦Ｄ）を充たす要素１／σ'Ｄ−、…、１／σ'Ｄに置き換えて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する（ＳＳ５１）。
【００５４】
参考非特許文献１（参考非特許文献１：浅野太著、「（音響テクノロジーシリーズ１６）音のアレイ信号処理−音源の定位・追跡と分離−」、コロナ社、2011年1月27日、p.111-p.113）によれば、特異値は対応する特異値ベクトルの張る空間における信号パワーと雑音パワーの和である。信号パワーは偏在する傾向があり、雑音パワーは均一に分布する傾向がある。そのため、小さい特異値に対応する特異ベクトルの張る空間では、雑音パワーの占める比率がより大きくなる。このことを念頭に、疑似逆行列推定手段５１は、上記実施例で使用している疑似逆行列
【００５５】
【数１８】
【００５６】
を変形する。式（２４）によれば、疑似逆行列の対角項は特異値の逆数になっている。このため、疑似逆行列は雑音パワー比の大きい成分に、より大きい値を乗算し、ノイズを増幅しやすいことが分かる。そこで、雑音パワー比の大きい成分に特異値の逆数よりも小さい値を乗算することで、ノイズ成分の増幅を抑える。そのために、疑似逆行列推定手段５１は、疑似逆行列Ｒ２＃の代わりに行列Ｑを用いる。
行列Ｑの対角項は、
【００５７】
【数１９】
【００５８】
に設定されている。これにより、残響除去の性能をほぼ維持しつつ、逆フィルタによるノイズ増幅を抑えることができる。行列Ｑの一具体形として、式（２５）のように１／σＤ−をもちいる方法が考えられる。
【００５９】
【数２０】
【００６０】
なおＤ−は、Ｊ−〜Ｊ＋の範囲から指定する。図１４に、特異値、特異値の逆数、行列Ｑの対角項の関係を示す。図１４（ａ）は特異値番号と特異値の関係を示す図である。図１４（ｂ）は特異値番号と所定値以上の特異値逆数を一定値に置き換えた特異値逆数（行列Ｑの対角項）の関係を示す図である。
【００６１】
以下、実施例１の残響除去装置１００、実施例２の残響除去装置２００について、シミュレーションにより検証した効果について図１７、図１８、図１９、図２０を参照して説明する。図１７はシュミレーションに用いたインパルス応答（ａ）と残響曲線（ｂ）を示す図である。図１８は従来技術の残響除去装置のシュミレーション結果を示す図である。図１９は実施例１の残響除去装置１００のシュミレーション結果を示す図である。図２０は実施例２の残響除去装置２００のシュミレーション結果を示す図である。シュミレーションに用いた音源は白色雑音であり、４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測したインパルス応答を５５００タップに打ち切って使用した。各マイクロホンには無相関の白色ノイズをＳＮ比２５ｄＢになるよう加えている。逆フィルタのフィルタ長は１チャネル当り１０００タップに設定した。相関関数を算出するときには、２０ｓの収音信号を使用した。図１８、図１９、図２０はいずれも（ａ）音響経路＋逆フィルタのインパルス応答、（ｂ）収音経路＋逆フィルタのインパルス応答の残響曲線（実線）、収音経路インパルス応答の残響曲線（点線）、（ｃ）逆フィルタ後の信号、（ｄ）逆フィルタ後のノイズである。なお、図１８〜図２２のグラフ（ａ）は縦軸を信号強度、横軸をタップ、グラフ（ｂ）は縦軸をエネルギー、横軸を時間、グラフ（ｃ）は縦軸を信号強度、横軸を時間、グラフ（ｄ）は縦軸を信号強度、横軸を時間、としてそれぞれプロットしたものである。図１８の従来法の処理結果では、逆フィルタをうまく求められていないため、音響経路＋逆フィルタのインパルス応答は、ピーク位置が拡散している。その残響曲線もほとんど減衰していない。また逆フィルタ後は、ノイズの振幅が信号の４倍にも達している。図１９の実施例１の処理結果では、逆フィルタがうまく求められており、音響経路＋逆フィルタのインパルス応答は、ピークが１つ立っている。その残響曲線も、８０サンプル目のところで−１０ｄＢ減衰している。そして逆フィルタ後は、ノイズの振幅が信号振幅の半分程度に抑えられている。図２０の実施例２の処理結果では、より良い逆フィルタが求められており、音響経路＋逆フィルタインパルス応答の残響曲線の減衰は、実施例１の処理結果と比較して２ｄＢ改善している。そして逆フィルタ後は、ノイズの振幅が信号振幅の１／３程度に抑えられている。
【００６２】
本実施例の残響除去装置２００によれば、実施例１の残響除去装置１００の効果に加え、疑似逆行列推定手段５１が雑音パワー比の大きい成分に乗算される特異値の逆数を当該特異値の逆数よりも小さな値に置き換えて疑似逆行列を推定するため、残響除去の性能をほぼ維持しつつ、逆フィルタによるノイズ増幅を抑えることができる。
【実施例３】
【００６３】
次に、図７、図８、図１５、図１６を参照して実施例１の残響除去装置１００に共通多項式を除去する機能を追加した実施例３に係る残響除去装置について説明する。図７は本実施例に係る残響除去装置３００の構成を示すブロック図である。図８は本実施例に係る残響除去装置３００の動作を示すフローチャートである。図１５は特異値の最大値から最小値に一定比率で減衰する関数ｐ〜（ｉ）を示す図である。図１６は関数ｐ〜（ｉ）とｉ番目の特異値σｉの比で定義される関数ｑ〜（ｉ）を示す図である。本実施例の残響除去装置３００は、Ｎチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部４０と、共通多項式除去部６０と、逆フィルタ部７０とを備える。信号部分空間抽出部３０は、特異値分解手段３１と、信号部分空間次元推定手段３２とを備える。逆フィルタ推定部４０は、疑似逆行列推定手段４１と、逆フィルタ係数推定手段４２とを備える。共通多項式除去部６０は、シルベスタ行列生成手段６１と、左ゼロ空間抽出手段６２と、共通多項式次元推定手段６３と、共通多項式フィルタ係数推定手段６４と、共通多項式除算手段６５とを備える。本実施例の残響除去装置３００と、実施例１の残響除去装置１００の違いは、実施例１の残響除去装置１００が備えない共通多項式除去部６０を、本実施例の残響除去装置３００が備える点のみである。従って、実施例１の残響除去装置１００と共通する番号を有する各構成部は、実施例１の残響除去装置１００において同一番号を有する各構成部と全く同じ機能を有するため説明を割愛する。
【００６４】
実施例１の残響除去装置１００は、音源信号が白色であること、すなわち音源信号のスペクトルがフラットなことを前提としている。音源信号が音声のように有色な場合には、実施例１の逆フィルタは音響経路に由来する残響を除去すると同時に、音源信号の周波数特性の平坦化もはかってしまう。有色音源信号に対応するには、逆フィルタｇ１（ｋ’）〜ｇＮ（ｋ’）で音源周波数特性を平坦化している部分、すなわち共通多項式を求めて除去すればよい。以下、共通多項式の次元を推定し、この推定次元をもとに参考非特許文献２のＱｉｕの手法で共通多項式を求め、第一の逆フィルタから共通多項式を除算する処理を説明する（参考非特許文献２：W. Qiu, Y. Hua, and K. Abed-Meraim, " Subspace Method for the Computation of the GCD of Polynomials" Automatica,1997,vol. 33, no. 4, pp. 741-743）。なお次元が既知のとき共通多項式を求めるアルゴリズムは、参考非特許文献２以外に参考非特許文献３などの手法があり、これらを代わりに使用してもよい（参考非特許文献３：R. M. Corless, S. M. Watt, and L. Zhi, QR "Factoring to Compute the GCD of Univariate Approximate Polynomials", IEEE Transactions on Signal Processing, 2004, vol. 52, no. 12, pp. 3394-3402）。
【００６５】
まず、シルベスタ行列生成手段６１は、逆フィルタ係数推定手段４２において推定された逆フィルタ係数ベクトル（以下、第１の逆フィルタ係数ベクトルという）から（２Ｌ−１）行ＮＬ列のシルベスタ行列を生成する（ＳＳ６１）。より具体的には、シルベスタ行列生成手段６１は、逆フィルタ係数推定手段４２で求められた第一の逆フィルタ係数ベクトルから、式（２６）（２７）で（２Ｌ−１）×ＮＬシルベスタ行列を生成する。
【００６６】
【数２１】
【００６７】
左ゼロ空間抽出手段６２は、生成されたシルベスタ行列を特異値分解して（２Ｌ−１）個の特異値を出力する（ＳＳ６２）。共通多項式の次元がＤ〜のとき、特異値分解は式（２８）のように分割される（非特許文献６）。
【００６８】
【数２２】
【００６９】
行列Ｕ〜０の縦ベクトルが行列Ｇの左ゼロ空間の基底ベクトルになり、その本数がＤ〜になる。特異値カーブは図１５のようになる。共通多項式次元推定手段６３は、左ゼロ空間抽出手段６２が出力した特異値の減衰傾向に基づいて共通多項式の次元を推定する（ＳＳ６３）。より具体的には、共通多項式次元推定手段６３は、行列Ｇの特異値σ〜１からσ〜２Ｌ−１に一定比率で減衰するカーブｐ〜（ｉ）を求める（図１５参照）。
【００７０】
【数２３】
【００７１】
次に、ｐ〜（ｉ）とσ〜ｉの比ｑ〜（ｉ）をとり、その最小値の位置から変曲点位置Ｉｏを求める（図１６参照）。
【００７２】
【数２４】
【００７３】
変曲点での値ｑ〜（Ｉ０）よりγ［ｄＢ］小さい値をとる二点のうち大きいほうの点Ｉ＋を求める（図１６参照）。そして共通多項式の次元をＤ〜＝２Ｌ−１−Ｉ＋とする。なおγは３〜２０ｄＢの範囲で指定する。共通多項式フィルタ係数推定手段６４は、シルベスタ行列の左ゼロ空間の基底ベクトルから共通多項式のフィルタ係数を推定する（ＳＳ６４）。より具体的には、共通多項式フィルタ係数推定手段６４は、左ゼロ空間の基底ベクトルから、式（３１−１〜３１−３）で行列Ｗを算出する。行列Ｗの最小固有値に対応する固有ベクトルが共通多項式のフィルタ係数ｃとなる。
【００７４】
【数２５】
【００７５】
共通多項式除算手段６５は、第１の逆フィルタ係数ベクトルをチャネル毎の共通多項式で除算して第２の逆フィルタ係数ベクトルを推定する（ＳＳ６５）。より具体的には、共通多項式除算手段６５は、逆フィルタ推定手段４２で求めた第一の逆フィルタを、チャネルごとに共通多項式で除算する。そのために式（３２）を解いて、チャネルごとに第２の逆フィルタ係数
【００７６】
【数２６】
【００７７】
を求める。
【００７８】
【数２７】
【００７９】
逆フィルタ部７０では、式（３３）により収音信号ｙ１（ｋ）〜ｙＮ（ｋ）に逆フィルタ係数を適用して出力信号ｘ＾（ｋ）を得る。
【００８０】
【数２８】
【００８１】
実施例３の残響除去装置３００について、シミュレーションにより効果を検証した。音源は音声であり、４チャネル収音信号は、残響時間４５０ｍｓの部屋において、サンプリング周波数８ｋＨｚで実測したインパルス応答を５５００タップに打ち切って使用した。各マイクロホンには無相関の白色ノイズをＳＮ比２０ｄＢになるよう加えている。逆フィルタのフィルタ長は１チャネル当り１０００タップに設定した。相関関数を算出するときには、２０ｓの収音信号を使用した。実施例２の手法を適用した結果を図２１に、実施例１の手法を適用した結果を図２２に示す。いずれも（ａ）音響経路＋逆フィルタのインパルス応答、（ｂ）収音経路＋逆フィルタのインパルス応答の残響曲線（実線）、収音経路インパルス応答の残響曲線（点線）、（ｃ）逆フィルタ後の信号、（ｄ）逆フィルタ後のノイズである。
【００８２】
実施例１の手法を適用した場合には、図２２（ａ）インパルス応答ではきれいなピーク形状が得られず、ピーク後部に成分が残留している。また図２２（ｂ）の残響曲線でも８０サンプル付近からなだらかに減衰しており、残響がうまく取れていないことが分かる。
【００８３】
しかし実施例３の手法を適用することで、図２１（ａ）の音響経路＋逆フィルタのインパルス応答ではきれいなピークが得られており、その残響曲線も８０サンプル付近で０ｄＢから−１０ｄＢに急激に減衰している、このことから、音源が音声信号であっても実施例３の手法により残響を効果的に除去できることが分かる。
【００８４】
本実施例の残響除去装置３００によれば、実施例１の残響除去装置１００の効果に加えて、音源信号が有色である場合であっても、共通多項式除去部６０が共通多項式を求めて除去することにより、有色音源信号にも対応することができる。
【実施例４】
【００８５】
以下、図９、図１０を参照して実施例２の残響除去装置２００に実施例３の残響除去装置３００が備える共通多項式除去部６０を加えた実施例４の残響除去装置について説明する。図９は本実施例に係る残響除去装置４００の構成を示すブロック図である。図１０は本実施例に係る残響除去装置４００の動作を示すフローチャートである。本実施例の残響除去装置４００は、実施例２の残響除去装置２００と共通してＮチャネルのマイクロホン１０−１〜１０−Ｎ（Ｎは２以上の整数）と、多チャネル相関行列算出部２０と、信号部分空間抽出部３０と、逆フィルタ推定部５０と、逆フィルタ部７０とを備える。また、本実施例の残響除去装置４００は、実施例３の残響除去装置３００と共通して共通多項式除去部６０を備える。従って、実施例２の残響除去装置２００と共通する番号を有する各構成部の機能は、実施例２の記載を参照されたい。また、実施例３の残響除去装置３００と共通する番号を有する各構成部の機能は、実施例３の記載を参照されたい。
【００８６】
また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００８７】
また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
【００８８】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【００８９】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００９０】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【００９１】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【特許請求の範囲】
【請求項１】
Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する多チャネル相関行列算出ステップと、
前記算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する特異値分解サブステップと、
前記出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する信号部分空間次元推定サブステップと、を有する信号部分空間抽出ステップと、
前記推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定サブステップと、
前記推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定サブステップと、を有する逆フィルタ推定ステップと、
前記推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、前記Ｎチャネルの収音信号のそれぞれに適用する逆フィルタステップと、
を有することを特徴とする残響除去方法。
【請求項２】
請求項１に記載の残響除去方法であって、
前記信号部分空間次元推定サブステップが、
前記ＮＬ個の特異値をσ１、σ２、…、σＮＬとし、ｊ（１≦ｊ≦ＮＬ）を変数として、
【数２９】
により定義される関数ｑ（ｊ）の変曲点位置ｊ＝Ｊ０から信号部分空間の次元を推定すること
を特徴とする残響除去方法。
【請求項３】
請求項１または２に記載の残響除去方法であって、
前記疑似逆行列推定サブステップが、
前記信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列の対角行列の要素である特異値の逆数１／σ１、１／σ２、…、１／σＤのうち、予め定めたしきい値以上となる要素１／σＤ−、…、１／σＤを１／σ'Ｚ＜１／σＺ（Ｄ−≦Ｚ≦Ｄ）を充たす要素１／σ'Ｄ−、…、１／σ'Ｄに置き換えて前記信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定すること
を特徴とする残響除去方法。
【請求項４】
請求項１から３の何れかに記載の残響除去方法であって、
共通多項式除去ステップをさらに有し、
前記共通多項式除去ステップが、
前記逆フィルタ係数推定サブステップにおいて推定された逆フィルタ係数ベクトル（以下、第１の逆フィルタ係数ベクトルという）から（２Ｌ−１）行ＮＬ列のシルベスタ行列を生成するシルベスタ行列生成サブステップと、
前記生成されたシルベスタ行列を特異値分解して（２Ｌ−１）個の特異値を出力する左ゼロ空間抽出サブステップと、
前記左ゼロ空間抽出サブステップが出力した特異値の減衰傾向に基づいて共通多項式の次元を推定する共通多項式次元推定サブステップと、
前記シルベスタ行列の左ゼロ空間の基底ベクトルから共通多項式のフィルタ係数を推定する共通多項式フィルタ係数推定サブステップと、
前記第１の逆フィルタ係数ベクトルを前記チャネル毎の共通多項式で除算して第２の逆フィルタ係数ベクトルを推定する共通多項式除算サブステップと、
を有することを特徴とする残響除去方法。
【請求項５】
請求項４に記載の残響除去方法であって、
前記共通多項式次元推定サブステップが、
前記シルベスタ行列を特異値分解して得られた（２Ｌ−１）個の特異値をσ〜１、σ〜２、…、σ〜２Ｌ−１とし、ｉ（１≦ｉ≦２Ｌ−１）を変数として、
【数３０】
により定義される関数ｑ〜（ｉ）の変曲点位置ｉ＝Ｉ０から共通多項式の次元を推定すること
を特徴とする残響除去方法。
【請求項６】
Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する多チャネル相関行列算出部と、
前記算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する特異値分解手段と、
前記出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する信号部分空間次元推定手段と、を備える信号部分空間抽出部と、
前記推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定手段と、
前記推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定手段と、を備える逆フィルタ推定部と、
前記推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、前記Ｎチャネルの収音信号のそれぞれに適用する逆フィルタ部と、
を有することを特徴とする残響除去装置。
【請求項７】
請求項１から５の何れかに記載の残響除去方法を実行すべき指令をコンピュータに対してするプログラム。
【請求項１】
Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する多チャネル相関行列算出ステップと、
前記算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する特異値分解サブステップと、
前記出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する信号部分空間次元推定サブステップと、を有する信号部分空間抽出ステップと、
前記推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定サブステップと、
前記推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定サブステップと、を有する逆フィルタ推定ステップと、
前記推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、前記Ｎチャネルの収音信号のそれぞれに適用する逆フィルタステップと、
を有することを特徴とする残響除去方法。
【請求項２】
請求項１に記載の残響除去方法であって、
前記信号部分空間次元推定サブステップが、
前記ＮＬ個の特異値をσ１、σ２、…、σＮＬとし、ｊ（１≦ｊ≦ＮＬ）を変数として、
【数２９】
により定義される関数ｑ（ｊ）の変曲点位置ｊ＝Ｊ０から信号部分空間の次元を推定すること
を特徴とする残響除去方法。
【請求項３】
請求項１または２に記載の残響除去方法であって、
前記疑似逆行列推定サブステップが、
前記信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列の対角行列の要素である特異値の逆数１／σ１、１／σ２、…、１／σＤのうち、予め定めたしきい値以上となる要素１／σＤ−、…、１／σＤを１／σ'Ｚ＜１／σＺ（Ｄ−≦Ｚ≦Ｄ）を充たす要素１／σ'Ｄ−、…、１／σ'Ｄに置き換えて前記信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定すること
を特徴とする残響除去方法。
【請求項４】
請求項１から３の何れかに記載の残響除去方法であって、
共通多項式除去ステップをさらに有し、
前記共通多項式除去ステップが、
前記逆フィルタ係数推定サブステップにおいて推定された逆フィルタ係数ベクトル（以下、第１の逆フィルタ係数ベクトルという）から（２Ｌ−１）行ＮＬ列のシルベスタ行列を生成するシルベスタ行列生成サブステップと、
前記生成されたシルベスタ行列を特異値分解して（２Ｌ−１）個の特異値を出力する左ゼロ空間抽出サブステップと、
前記左ゼロ空間抽出サブステップが出力した特異値の減衰傾向に基づいて共通多項式の次元を推定する共通多項式次元推定サブステップと、
前記シルベスタ行列の左ゼロ空間の基底ベクトルから共通多項式のフィルタ係数を推定する共通多項式フィルタ係数推定サブステップと、
前記第１の逆フィルタ係数ベクトルを前記チャネル毎の共通多項式で除算して第２の逆フィルタ係数ベクトルを推定する共通多項式除算サブステップと、
を有することを特徴とする残響除去方法。
【請求項５】
請求項４に記載の残響除去方法であって、
前記共通多項式次元推定サブステップが、
前記シルベスタ行列を特異値分解して得られた（２Ｌ−１）個の特異値をσ〜１、σ〜２、…、σ〜２Ｌ−１とし、ｉ（１≦ｉ≦２Ｌ−１）を変数として、
【数３０】
により定義される関数ｑ〜（ｉ）の変曲点位置ｉ＝Ｉ０から共通多項式の次元を推定すること
を特徴とする残響除去方法。
【請求項６】
Ｎチャネル（Ｎは２以上の整数）の収音信号間の相関関数から、Ｌ行Ｌ列（Ｌは逆フィルタの１チャネルあたりのタップ数）の相関行列を生成し、当該相関行列からＮＬ行ＮＬ列の多チャネル相関行列を算出する多チャネル相関行列算出部と、
前記算出された多チャネル相関行列を特異値分解して、ＮＬ個の特異値を出力する特異値分解手段と、
前記出力されたＮＬ個の特異値の減衰傾向に基づいて信号部分空間の次元Ｄを推定する信号部分空間次元推定手段と、を備える信号部分空間抽出部と、
前記推定された信号部分空間の次元Ｄに基づいて信号部分空間成分のみからなる多チャネル相関行列の疑似逆行列を推定する疑似逆行列推定手段と、
前記推定された疑似逆行列に基づいて逆フィルタ係数ベクトルを推定する逆フィルタ係数推定手段と、を備える逆フィルタ推定部と、
前記推定された逆フィルタ係数ベクトルを各チャネルのフィルタ係数ベクトルに分解し、前記Ｎチャネルの収音信号のそれぞれに適用する逆フィルタ部と、
を有することを特徴とする残響除去装置。
【請求項７】
請求項１から５の何れかに記載の残響除去方法を実行すべき指令をコンピュータに対してするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【公開番号】特開２０１３−１１３８６６（Ｐ２０１３−１１３８６６Ａ）
【公開日】平成２５年６月１０日（２０１３．６．１０）
【国際特許分類】
【出願番号】特願２０１１−２５７１３５（Ｐ２０１１−２５７１３５）
【出願日】平成２３年１１月２５日（２０１１．１１．２５）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成２５年６月１０日（２０１３．６．１０）
【国際特許分類】
【出願日】平成２３年１１月２５日（２０１１．１１．２５）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]