JP5700850B2 - 遅延推定方法とその方法を用いたエコー消去方法と、それらの装置とプログラムとその記録媒体 - Google Patents

Description

本発明は、例えばエコー消去装置等に利用することができる遅延推定方法とその方法を用いたエコー消去方法と、それらの装置とプログラムとその記録媒体に関する。

エコー消去装置は、受話信号をスピーカから音響信号として出力し、マイクロホンで収音した送話信号に回り込んでくるエコー信号成分を適応フィルタにより推定し、送話信号からエコー信号成分を消去するものである。ディジタルテレビやスマートフォンなどの機器では、音声を途切れないようにする目的で受話信号と送話信号のそれぞれにバッファが挿入される場合がある。また、ディジタルテレビの映像のコーデックで発生する映像と音声のずれを調整するリップシンク機能で用いられるバッファ等が知られている。

図２０に、受話信号と送話信号の双方にバッファを挿入したエコー消去装置９００の概略的な機能構成例を示す。エコー消去装置９００は、適応フィルタ９０１と、加算回路９０２と、受話信号バッファ９０３と、送話信号バッファ９０４と、を備える。受話信号バッファ９０３は、遠端話者からの受話信号を、数１０ｍｓ〜数１００ｍｓの間保存し、保存した順に受話信号をスピーカに出力する。送話信号バッファ９０４は、マイクロホンで収音した近端話者からの送話信号を、一時的に保存し、保存した順に送話信号を加算回路９０２に出力する。なお図２０において、ディジタル信号をアナログ信号に変換するＡ/Ｄ変換器や、その逆の変換をするＤ/Ａ変換器等の機能部の表記は省略している。

この受話信号バッファ９０３と送話信号バッファ９０４は、例えば複数のアプリケーションを並行して実行するためにＣＰＵの負荷が重くなり、エコー消去処理が間に合わない場合でも音声が途切れないようにする目的で設けられたものである。よって、このバッファ部分で遅延が発生する。一般にエコー消去装置の持つ適応フィルタは１００ｍｓ程度の長さであり、この適応フィルタの長さを超える遅延時間のエコー成分を消去することができない。

そこで、従来から図２１に示すような遅延推定を用いたエコー消去装置９５０が考えられている。エコー消去装置９５０は、上記したエコー消去装置９００に対して遅延推定部９５１と遅延挿入部９５２を備える点で異なる。遅延推定部９５１は、受話信号バッファ９０３に入力される前の受話信号である遅延前信号と、送話信号バッファ９０４から出力される送話信号である遅延後信号との間の遅延時間を推定し、推定した遅延時間を適応フィルタの前の遅延挿入部９５２に与えることで、大きな遅延時間があっても適応フィルタ９０１でエコー信号成分の消去を可能にするものである。

図２２に、特許文献１に開示された従来の遅延推定装置９６０の機能構成図を示して、その動作を簡単に説明する。ＦＦＴ部１１１_１〜１１１_Ｍは、マイクロホン１１_１〜１１_Ｍからの収音信号を周波数領域に変換する。白色化部１１２_１〜１１２_Ｍは、周波数領域に変換された収音信号を、周波数スペクトルで白色化（フラット）する。次に、マイクロホン対選択部１１３は、白色化部１１２_１〜１１２_Ｍの出力信号のうち２つを選択する。乗算部１１４は、マイクロホン対選択部１１３で選ばれた信号のうち一方だけ共役をとり、２つの信号を周波数成分ごとに乗算し、クロススペクトルを求める。乗算部１１４の出力信号をＩＦＦＴ部１１５により、時間領域に変換し、白色化相互相関を求める。次に、最大ピーク検出部１１６で、ＩＦＦＴ部１１５の出力の相互相関の最大ピークを検出し、その最大ピークの地点を収音信号間遅延時間差として出力する。

遅延推定装置９６０の考えを、上記した遅延推定部９５１に適用すると、遅延前信号が例えばＦＦＴ部１１１_１、遅延後信号が例えばＦＦＴ部１１１_２、にそれぞれ入力される関係になる。以下、数式を用いて更に遅延推定装置９６０の動作原理を説明する。

まず、遅延前信号の時間領域信号をｘ（ｔ）、遅延後信号の時間領域信号をｙ（ｔ）、推定したい遅延経路のインパルス応答をｈ（ｔ）とすれば式（１）の関係が成立する。

ここで＊は畳み込み演算を表す。

遅延経路のインパルス応答ｈ（ｔ）が求められれば、そのピークの位置から遅延を推定することができる。ｈ（ｔ）を求めるには、観測可能な遅延前信号の時間領域信号ｘ（ｔ）と、遅延後信号の時間領域信号ｙ（ｔ）を用いて、式（１）をｈ（ｔ）について解けばよい。しかし、式（１）は畳み込み演算となっているため、解くには信号点数の２乗のオーダの演算が必要となる。そこで、遅延推定装置９６０は、式（１）の信号を周波数領域の信号に変換して解くことで計算量を削減する。

遅延前信号の周波数領域信号をＸ（ω）、遅延後信号の周波数領域信号をＹ（ω）、推定したい遅延経路の周波数領域のインパルス応答をＨ（ω）とすれば式（１）は式（２）に変形される。

式（２）をＨ（ω）について解けば、式（３）となる。

式（３）によりＨ（ω）を求め、Ｈ（ω）を逆ＦＦＴすれば、ｈ（ｔ）を求めることができる。ＦＦＴの演算量のオーダはＮ・ｌｏｇ_２Ｎであるので、時間領域で解くよりも少ない演算量となる。

具体的な計算方法について説明する。上記したＮは、遅延前信号と遅延後信号をＦＦＴ処理するデータ点数であり、白色化相互相関を求める信号長である。乗算部１１４は、遅延前信号の周波数領域信号Ｘ（ω）の共役を取ったものと遅延後信号の周波数領域信号Ｙ（ω）とを周波数ビンごとに乗算して、式（３）の分子（クロススペクトル）を計算する。また、白色化部１１２において計算された遅延前信号の周波数領域信号Ｘ（ω）のノルムの２乗（式（３）の分母）の逆数を、クロススペクトルに乗算する。乗算部１１４の出力信号を、ＩＦＦＴ部１１５で、Ｎ点逆ＦＦＴすれば白色化相互相関が求められ、これが遅延経路のインパルス応答ｈ（ｔ）となる。最大ピーク検出部１１６は、その遅延経路のインパルス応答ｈ（ｔ）の最大ピーク位置を求め、推定遅延量として出力する。

このように、従来の遅延推定装置９６０は、Ｎ点のＦＦＴと逆ＦＦＴを用いて計算することで、時間領域で計算するよりも演算量を削減することができる。

特開２００７−８１４５５号公報

しかしながら、正確な遅延推定をするためには、数１００ｍｓ以上の信号長で相互相関を計算する必要があり、ＦＦＴを用いたとしても多くの演算量（具体的な計算量は後述する）が必要となる。

本発明は、この課題に鑑みてなされたものであり、正確な遅延推定を少ない演算量で実現する遅延推定方法とその方法を用いたエコー消去方法と、それらの装置とプログラムとその記録媒体を提供することを目的とする。

本発明の遅延推定方法は、第１の周波数分析過程と、第２の周波数分析過程と、第１のバッファ過程と、第２のバッファ過程と、クロススペクトル計算過程と、白色化係数計算過程と、白色化過程と、時間領域変換過程と、時間シフト加算過程と、最大ピーク位置探索過程と、を備える。第１の周波数分析過程は、白色化相互相関を求める時間幅の遅延前信号から採取されるサンプル点数をＮとした時に、上記遅延前信号を上記Ｎ点よりも少ない２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理をＭ回単位で繰り返す。第２の周波数分析過程は、遅延前信号を遅延させた遅延後信号を２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理を、上記Ｍフレームに遅延を推定する範囲であるＫフレームを加えたＭ＋Ｋ−１回の単位で繰り返す。第１のバッファ過程は、第１の周波数分析過程で周波数領域の信号に変換した遅延前信号を上記Ｍ個分保存する。第２のバッファ過程は、第２の周波数分析過程で周波数領域の信号に変換した遅延後信号を、上記Ｍ＋Ｋ−１個分保存する。クロススペクトル計算過程は、第１のバッファ過程で保存された遅延前信号の共役を取った信号と、第２のバッファ段階で保存された遅延後信号とを同一周波数毎に乗算し、当該乗算結果をフレームについて加算し上記Ｋ要素までのクロススペクトルを計算する。白色化係数計算過程は、第１のバッファ過程で保存された遅延前信号のノルムの２乗の逆数を白色化係数として計算する。白色化過程は、上記クロススペクトルに上記白色化係数を乗算する。時間領域変換過程は、白色化過程で白色化されたクロススペクトルを、時間領域に変換する。時間シフト加算過程は、時間領域変換過程の出力を時間シフトしながら第１要素から上記Ｋ要素まで加算して保存する。最大ピーク位置探索過程は、時間シフト加算過程で保存された信号の最大ピーク位置を推定遅延量として求める。

また、この発明のエコー消去方法は、上に示した遅延推定方法に加えて、遅延挿入過程と周波数領域適応フィルタ過程と、加算過程と、を備える。遅延挿入過程は、上記した第１の周波数分析過程で周波数領域の信号に変換された受話信号を、上記遅延推定方法で推定された推定遅延量分の時間遅延させる。周波数領域適応フィルタ過程は、遅延挿入過程で遅延が挿入された受話信号を入力として擬似エコー信号を生成する。加算過程は、周波数領域の遅延後信号から擬似エコー信号を除去する。

本発明の遅延推定方法によれば、白色化相互相関を求める時間長から得られるサンプル数Ｎよりも少ない数の２Ｎ/Ｍ点の周波数分析と時間領域変換の演算、つまり従来技術よりも少ない演算量で推定遅延量を求めることができる。

また、本発明のエコー消去方法によれば、この発明の遅延推定方法で推定した遅延時間を、適応フィルタの前に与えることで、長い遅延があってもエコー成分を消去することが可能である。また、周波数分析過程を共用することで演算量を大幅に削減することができる。

この発明の遅延推定装置１００の機能構成例を示す図。 遅延推定装置１００の動作フローを示す図。 サンプル点数Ｎと、Ｍとの関係を説明する図。 遅延推定装置１００のより具体的な動作フローを示す図。 白色化係数計算部４０の機能構成例を示す図。 クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）の第１要素の計算を示す図。 クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）の第２要素の計算を示す図。 時間シフト加算部７０の動作を示す図。 白色化係数計算部２４０の機能構成例を示す図。 白色化係数計算部２４０の動作フローを示す図。 遅延推定装置３００の機能構成例を示す図。 推定区間設定部３０１の機能構成例を示す図。 推定区間設定部３０１の動作フローを示す図。 マイクロホンで収音する収音信号のレベルを模式的に表した図。 従来の遅延推定装置（ＦＦＴ点数４０９６点）の推定結果の例を示す図。 従来の遅延推定装置（ＦＦＴ点数１０２４点）の推定結果の例を示す図。 この発明の遅延推定装置の推定結果の例を示す図。 遅延推定の正解率の比較を示す図。 この発明のエコー消去装置４００の機能構成例を示す図。 従来のエコー消去装置９００の機能構成を示す図。 従来のエコー消去装置９５０の機能構成を示す図。 従来の遅延推定装置９６０の機能構成を示す図。

以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。

図１に、この発明の遅延推定装置１００の機能構成例を示す。その動作フローを図２に示す。遅延推定装置１００は、第１の周波数分析部１０と、第２の周波数分析部１２と、第１のバッファ部２０と、第２のバッファ部２２と、クロススペクトル計算部３０と、白色化係数計算部４０と、白色化部５０と、時間領域変換部６０と、時間シフト加算部７０と、最大ピーク位置探索部８０と、制御部９０と、を備える。遅延推定装置１００は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。

第１の周波数分析部１０は、白色化相互相関を求める時間幅のサンプル点数をＮとした時に、遅延前信号を上記Ｎ点よりも少ない２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理をＭ回単位で繰り返す（ステップＳ１０）。第１のバッファ部２０は、第１の周波数分析部１０で周波数領域に変換された遅延前信号をＭ個分保存する（ステップＳ２０）。

第２の周波数分析部１２は、遅延前信号を遅延させた遅延後信号を２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理を、上記Ｍフレームに遅延を推定する範囲であるＫフレームを加えたＭ＋Ｋ−１回の単位で繰り返す（ステップＳ１２）。第２のバッファ部２２は、第２の周波数分析部１２で周波数領域に変換された遅延後信号を、Ｍ＋Ｋ−１個分保存する（ステップＳ２２）。

クロススペクトル計算部３０は、第１のバッファ部２０で保存された信号の共役を取った遅延前信号と第２のバッファ部２２で保存された遅延後信号とを、同一周波数毎に乗算し、当該乗算結果をＫ要素まで加算したクロススペクトルを計算する（ステップＳ３０）。白色化係数計算部４０は、第１のバッファ部２０で保存された遅延前信号のノルムの２乗の逆数を白色化係数として計算する（ステップＳ４０）。

白色化部５０は、クロススペクトルに白色化係数を乗算する（ステップＳ５０）。時間領域変換部６０は、白色化部５０で白色化されたクロススペクトルを、時間領域に変換する（ステップＳ６０）。時間シフト加算部７０は、時間領域変換部６０の出力信号を時間シフトしながら第１要素からＫ要素まで加算して保存する（ステップＳ７０）。最大ピーク位置探索部８０は、時間シフト加算部７０で保存された信号の最大ピーク位置を推定遅延量として求める（ステップＳ８０）。制御部９０は、遅延推定動作を終了するまで、或いは遅延前信号が無くなるまでステップＳ１０〜Ｓ８０の処理を繰り返すように、各部の動作を制御する。

ここで、図３を参照して第１の周波数分析部１０と第２の周波数分析部１２で周波数分析処理するサンプル点数Ｎ、及びＭについて説明する。周波数分析の手法は、ＦＦＴ（Fast Fourier Transform）、ＤＦＴ（Discrete Fourier Transformation）、ウエーブレット変換等の何れの周知な方法を用いることができる。以降の説明ではＦＦＴを用いた例で説明を行う。

白色化相互相関を取る信号時間幅を、例えば５１２ｍｓ、連続信号を離散値化するサンプリング周波数を８ｋＨｚと、仮定する。５１２ｍｓの時間幅で得られるサンプル数は０.５１２×８ｋ＝４０９６となる。従来技術の遅延推定法では、その仮定においてはサンプル数Ｎ＝４０９６のＦＦＴ処理を行っていた。

その従来技術に対してこの発明は、ＮをＭ毎に区切った（Ｎ/Ｍ）フレーム単位で、遅延量を推定する。つまり、サンプル数２Ｎ/Ｍで、シフト量Ｎ/ＭのＦＦＴ処理を、Ｍ回（この例ではＭ＝６４）繰り返す処理を行う点で、従来技術と異なる。

このように遅延推定装置１００では、信号長Ｎよりも短い（２Ｎ/Ｍ）点のＦＦＴをＭ回演算させることで、遅延前信号と遅延後信号を周波数領域の信号に変換する。ここでＭは整数であり、（２Ｎ/Ｍ）を、遅延推定装置１００を組み込むエコー消去装置で使用するＦＦＴと同じサンプル点数になるように、ＮとＭを設定する。ｍ回目（ｍはＭより小さい任意の整数）に周波数領域に変換された遅延前信号をＸ′_ｍ（ω）と表し、ｍ回目に周波数領域に変換された遅延後信号をＹ′_ｍ（ω）と表し、遅延経路のインパルス応答を（Ｎ/Ｍ）毎に区切って周波数領域に変換したものをＨ′_ｍ（ω）と表せば、次式の関係が成立する。

ここでＫは、遅延を推定するフレーム範囲を決定する整数であり、（ＫＮ/Ｍ）が、遅延の推定範囲のサンプル数となる。式（４）の各行列を式（５）〜式（７）で表せば、

式（４）は式（８）で表せる。

式（８）をＨ″（ω）について解けば式（９）となる。

^Ｈは行列の共役転置を表す。

Ｍ個の周波数領域に変換された遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）と、Ｍ＋Ｋ−１個の周波数領域に変換された遅延後信号Ｙ′_１（ω），…，Ｙ′_Ｍ（ω）から、式（９）を用いてＨ″（ω）を求めることができる。しかし、{Ｘ″^Ｈ（ω）Ｘ″′（ω）}^−１の計算については、Ｍ行Ｍ列の逆行列を計算する必要があり、演算量が多い。そこで、近似を用いて演算を省略する。Ｘ″^Ｈ（ω）Ｘ″′（ω）の各要素を計算すれば式（１０）となる。

ここで、式（１０）の対角成分は全て同じ値Ｘ′^＊ _１（ω）Ｘ′_１（ω）＋…＋Ｘ′^＊ _Ｍ（ω）Ｘ′_Ｍ（ω）である。この値はベクトル{Ｘ′_１（ω）＋…＋Ｘ′_Ｍ（ω）}のノルムの２乗に等しい。これに対して対角以外の成分については、同一時刻ではない信号の掛け算となっている。

遅延前信号ｘ（ｔ）の自己相関のＮ/Ｍサンプル以上の箇所が、信号パワーに比べ十分に小さいと仮定できれば、式（１０）の対角成分以外は対角成分に比べ十分に小さい値となり、式（１１）に示すように対角成分のみで近似することができる。一般に音声では、自己相関は数１０ｍｓで小さくなるので近似が可能である。

式（１１）を用いれば、式（９）の逆行列はスカラーの積で表すことができ式（１２）で近似される。

式（１２）を用いることで遅延経路の応答Ｈ″（ω）を、演算量を大幅に削減して求めることができる。遅延経路の応答Ｈ″（ω）を時間領域に戻すには、Ｈ″（ω）のＫ個の要素に対して、（２Ｎ/Ｍ）点の逆ＦＦＴをそれぞれ実施して（Ｎ/Ｍ）点ずつシフトして加算すればよい。

式（１２）の遅延経路の応答Ｈ″（ω）の具体的な計算方法と、遅延経路の応答Ｈ″（ω）を時間領域に変換する方法について図４に動作フローを示し、図１と対応させながら具体的に説明する。遅延前信号を周波数領域の信号に変換する第１の周波数分析部１０は、遅延前信号を２Ｎ/Ｍ点のＦＦＴで周波数領域に変換する（ステップＳ１０、図４）。信号のシフト量はＮ/Ｍサンプルであり、Ｍ回のＦＦＴを実施する。第１のバッファ部２０は周波数領域に変換された遅延前信号のＭ回分を保存する（ステップＳ２０）。これらの処理により、周波数領域に変換された遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）が、第１のバッファ部２０に保存される。

遅延後信号を周波数領域に変換する第２の周波数分析部１２は、遅延後信号を２Ｎ/Ｍ点のＦＦＴで周波数領域の信号に変換する（ステップＳ１２）。信号のシフト量はＮ/Ｍサンプルであり、Ｍ＋Ｋ−１回のＦＦＴを実施する。第２のバッファ部２２は周波数領域に変換された遅延後信号のＭ＋Ｋ−１回分を保存する（ステップＳ２２）。これらの処理により、周波数領域に変換された遅延後信号Ｙ′_１（ω），…，Ｙ′_{Ｍ＋Ｋ−１}（ω）が、第２のバッファ部２２に保存される。

白色化係数計算部４０は、第１のバッファ部２０に保存されている遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）から、白色化係数を計算する（ステップＳ４０）。図５に白色化係数計算部４０のより具体的な機能構成例を示す。白色化係数計算部４０は、ノルム計算手段４１と、逆数計算手段４２と、を備える。ノルム計算手段４１は、遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）のノルムの２乗Ｘ′^＊ _１（ω）Ｘ′_１（ω）＋…＋Ｘ′^＊ _Ｍ（ω）Ｘ′_Ｍ（ω）を計算する。逆数計算手段４２は、ノルムの２乗の逆数を計算して白色化係数として出力する。

クロススペクトル計算部３０は、第１のバッファ部２０に保存されている遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）と、第２のバッファ部２２に保存されている遅延後信号Ｙ′_１（ω），…，Ｙ′_{Ｍ＋Ｋ−１}（ω）から、クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）を計算する。

図６と図７を用いてクロススペクトルの計算方法を説明する。図６に、クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）の第１要素の計算を示す。第１要素とは、第１のバッファ部２０と第２のバッファ部２２の読み出し位置を先頭に設定した場合である（ステップＳ３１）。第１のバッファ部２０と第２のバッファ部２２の同じフレーム番号の信号について、周波数ビン毎にそれぞれ乗算（ステップＳ３２）し、フレーム番号について加算する（ステップＳ３３）。第１のバッファ部２０のデータは共役^＊を取ってから乗算される。なお、フレーム番号は、何回目のＦＦＴの結果であるかを表す番号である。図７に、クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）の第２要素の計算を示す。第２要素とは、第１のバッファ部２０の読み出し位置は変えずに、第２のバッファ部２２の読み出し位置を１フレームずらして、周波数ビン毎にそれぞれ乗算し、フレーム番号について加算をしたものである。この第２バッファ部２２の読み出し位置を１フレームずらす処理は、フレームをＭ＋Ｋ−１個ずらすまで繰り返される。

白色化部５０は、クロススペクトル計算部３０が出力するクロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）に白色化係数１/{Ｘ′^＊ _１（ω）Ｘ′_１（ω）＋…＋Ｘ′^＊ _Ｍ（ω）Ｘ′_Ｍ（ω）}を乗じて遅延経路の応答Ｈ″（ω）を計算する（ステップＳ５０）。

時間領域変換部６０は、遅延経路の応答Ｈ″（ω）の各要素を、２Ｎ/Ｍ点の逆ＦＦＴで時間領域の信号に変換する（ステップＳ６０）。時間シフト加算部７０は、時間領域変換部６０の出力信号を時間シフトしながら加算して保存する（ステップＳ７０）。図８に、時間シフト加算部７０の動作を示す。時間シフト加算部７０は、時間領域変換部６０が出力するクロススペクトルの第１要素〜第Ｋ要素を、Ｎ/Ｍサンプルずつずらして加算し、時間領域の白色化相互相関ｈ（ｔ）（推定した遅延経路のインパルス応答）を求める。

遅延を推定する範囲であるＫフレームについてのステップＳ３２〜Ｓ７０の処理が終了するまで、第２のバッファ部２２の読み出し位置を１フレームずつ進めながら（ステップＳ３４）その処理を繰り返す（ステップＳ９０のＮＯ）。クロススペクトルＸ″^Ｈ（ω）Ｙ″（ω）の計算は、フレームの第Ｍ＋Ｋ−１要素まで繰り返される（ステップＳ９０のＹＥＳ）。

最大ピーク位置探索部８０は、時間領域の白色化相互相関ｈ（ｔ）の最大ピーク位置を探索し、その位置を推定遅延量として出力する（ステップＳ８０）。

以上説明したようにこの発明の遅延推定方法では、第１の周波数分析部１０におけるサンプル点数２Ｎ/ＭのＭ回のＦＦＴ処理と、第２の周波数分析部１２におけるサンプル点数２Ｎ/ＭのＭ＋Ｋ−１回のＦＦＴ処理と、時間領域変換部６０におけるサンプル点数２Ｎ/ＭのＫ回の逆ＦＦＴと、の積和演算で推定遅延量を求めることができる。例えば、Ｎ＝４０９６、Ｍ＝１２８、Ｋ＝１６、とした時の各部で計算する積和演算回数は、第１の周波数分析部１０で４９１５２回、第２の周波数分析部１２で５４９１２回、時間領域変換部６０で６１４４回、それらを合計すると１１０２０８回の積和演算回数となる。

これに対して従来技術の遅延推定方法では、Ｎ点のＦＦＴが２回と、Ｎ点の逆ＦＦＴが１回必要なので、その積和演算回数は１４７４５６回（３Ｎｌｏｇ_２Ｎ）である。このようにこの発明の遅延推定方法は、積和演算の回数を約２５％削減することができる。

図９に、この発明の実施例２の遅延推定装置２００の白色化係数計算部２４０の機能構成を示す。遅延推定装置２００は、白色化係数計算部２４０の構成のみが遅延推定装置１００と異なり、雑音に対する耐性を向上させたものである。

遅延前信号の周波数成分のうちパワーの小さい成分は、雑音の影響を受け易いため、遅延量推定の推定精度の低下の原因となる。白色化係数計算部２４０は、このようなパワーの小さな周波数成分に対する白色化係数を小さく設定することで、雑音に対する耐性を向上させる。

白色化係数計算部２４０は、ノルム計算手段４１、フルバンドパワー計算手段２４１、係数乗算手段２４２、加算手段２４３、逆数計算手段４２、を備える。図１０に、各手段の動作ステップを示す。ノルム計算手段４１は、実施例１と同様に周波数領域に変換された遅延前信号Ｘ′_１（ω），…，Ｘ′_Ｍ（ω）のノルムの２乗Ｘ′^＊ _１（ω）Ｘ′_１（ω）＋…＋Ｘ′^＊ _Ｍ（ω）Ｘ′_Ｍ（ω）を計算する（ステップＳ４１）。

フルバンドパワー計算手段２４１は、ノルムの２乗を周波数について加算平均したフルバンド平均パワーを計算する（ステップＳ２４１）。係数乗算手段２４２は、フルバンドパワー計算手段２４１の出力に、予め設定した定数αを乗じる（ステップＳ２４２）。定数αを小さくすると白色化が厳密に行われるようになり、定数αを大きくすると雑音に対する耐性を向上させることができる。定数αは、白色化の効果と、雑音に対する耐性のバランスのとれた適切な値に設定する必要がある。

加算手段２４３は、ノルム計算手段４１の出力と係数乗算部の出力を加算する（ステップＳ２４３）。逆数計算手段４２は、加算手段２４３の出力の逆数を計算して、白色化係数として出力する（ステップＳ４２）。

白色化係数計算部２４０で計算される白色化係数β（ω）は式（１３）で表せる。

ここでＷは使用する周波数の最大値である。

式（１３）において、分母にフルバンドの平均パワーが加算されたことで、パワーの小さい周波数でも分母が極端に小さな値にならなくなる。それにより、パワーの小さい周波数成分に対して白色化係数が非常に大きくなることが防止でき、雑音に対する耐性を向上させることができる。

図１１に、この発明の実施例３の遅延推定装置３００の機能構成例を示す。遅延推定装置３００は、遅延推定装置１００，２００の構成に、推定区間設定部３０１を追加したものである。

推定区間設定部３０１は、遅延前信号を入力として遅延量推定に適した信号区間を設定する。遅延量推定に適した信号区間で遅延量の推定をすることで、推定誤りを減らして推定精度を高くすることができる。

まず、どのような信号区間が遅延量の推定に適した区間であるかについて説明する。遅延量の推定において、推定精度を低下させる要因としては、マイクロホンに収音される周囲騒音と部屋の残響音が挙げられる。これらの影響を少なくするためには、スピーカからの出力音量が大きく、なるべく残響音がマイクロホンに収音されない区間を選択すれば良い。

図１２に、マイクロホンに収音される収音信号のパワーを模式的に表す。図１２の横軸は経過時間、縦軸はパワーである。周囲騒音のパワーは、ほぼ一定レベルでマイクロホンに収音される（黒の塗つぶし部分）。周囲騒音の影響の少ないところを見つけるには、振幅の大きな部分を検出すれば良い。よって、雑音レベルを推定し、それに係数を掛け合わせた閾値を越える区間を抽出することで、周囲騒音の影響を受け難い区間を得ることができる。

残響音（細かい破線）は、スピーカから出てマイクロホンに直接到達する直接音（細い実線）から少し送れて到達する。したがって、残響音成分の少ない区間は、音声の立上り部分に存在する。その区間は、遅延前信号の振幅から残響音のレベルを推定し、その残響音のレベルよりも大きいレベルの区間である。

残響音成分が少なくて周囲騒音の影響を受け難い区間、つまり、遅延推定精度の高い区間は、収音信号のパワーが雑音レベルよりも大きく、且つ残響音信号のパワーよりも大きい区間である。推定区間設定部３０１は、その区間を設定する。

図１３に、推定区間設定部３０１の機能構成例を示す。その動作フローを図１４に示す。推定区間設定部３０１は、雑音レベル推定手段３０１０と、残響レベル推定手段３０１１と、閾値設定手段３０１２と、区間検出手段３０１３と、で構成される。

雑音レベル推定手段３０１０は、遅延前信号を入力として、当該遅延前信号の振幅のディップホールドを行って雑音レベルを推定する（ステップＳ３０１０）。雑音レベル推定手段３０１０は、遅延前信号のレベルＰ（ｔ）のディップホールドを行うことで、雑音レベルＮ（ｔ）を推定する。これは、例えば式（１４）を用いて行われる。

過去の推定ノイズレベルＮ（ｔ−１）が、遅延前信号のレベルＰ（ｔ）よりも大きい場合は、推定ノイズレベルに遅延前信号のレベルを代入する。過去の推定ノイズレベルＮ（ｔ−１）が、遅延前信号のレベルＰ（ｔ）よりも小さい場合は、過去の推定ノイズレベルＮ（ｔ−１）に１以上の定数ｕを乗じてわずかにノイズレベルを上昇させる。ここでｕは１以上の定数であり事前に設定する。ｕは推定ノイズレベルの上昇係数であり１に近いほど緩やかな上昇となり、ディップホールド効果が得られる。

残響レベル推定手段３０１１は、遅延前信号を入力として、当該遅延前信号の振幅のピークホールドを行って残響音レベルＲ（ｔ）を推定する（ステップＳ３０１１）。残響音レベルＲ（ｔ）の推定は、遅延前信号のレベルＰ（ｔ）をピークホールドすることで行う（式（１５））。

過去の推定残響音レベルＲ（ｔ−１）が、遅延前信号のレベルＰ（ｔ）よりも大きい場合は、推定残響音レベルＲ（ｔ−１）に１未満の定数ｖを乗じて残響音レベルを減少させる。過去の推定残響音レベルＲ（ｔ−１）が、遅延前信号のレベルＰ（ｔ）よりも小さい場合は、残響音レベルに遅延前信号のレベルを代入する。ここでｖは、推定残響音レベルの減衰係数であり残響時間に対応した値を設定する。

閾値設定手段３０１２は、推定雑音レベルＮ（ｔ）に１以上の定数を乗じた値と、推定残響音レベルＲ（ｔ）に１未満の定数を乗じた値と、を比較し大きい方を閾値として設定する（ステップＳ３０１２）。区間検出手段３０１３は、その閾値と遅延前信号の振幅とを比較して閾値よりも遅延前信号の振幅が大きい区間を、遅延推定する区間として検出する（ステップＳ３０１３）。

遅延推定装置１００，２００は、推定区間設定部３０１が出力する推定区間で遅延推定を行う。このように遅延推定装置３００は、遅延推定に適した遅延前信号から遅延推定を行うので、周囲雑音や残響音の影響を受けずに高い精度で遅延推定を行うことができる。
〔評価実験〕
この発明の遅延推定方法の効果を確認する目的で評価実験を行った。図１５に、従来の遅延推定装置で、ＦＦＴの点数を４０９６点で遅延推定をした結果を示す。横軸は時刻（ｍｓ）、縦軸は遅延推定結果（ｍｓ）である。この時の実際の遅延は６５ｍｓであり、多くの遅延推定結果（＋）が６５ｍｓ付近に集中しているが、推定誤りも多い。

図１６に、ＦＦＴの点数を１０２４点として従来の遅延推定方法で遅延推定した結果を示す。遅延推定結果（＋）が一様に分布してしまい推定が全く出来ていないことが分かる。この原因は、残響時間に対して相関計算に使用する信号長が短すぎることになる。このように、従来の遅延推定方法では、精度の高い推定を行うためには長い点数のＦＦＴが必要であり、多くの演算が必要となる。

図１７に、この発明の実施例３で説明した遅延推定方法で、ＦＦＴの点数１２８点で３５フレーム分の相関を計算して遅延推定を行った結果を示す。全ての遅延推定結果が実際の遅延である６５ｍｓ付近にあり、正確に遅延推定が出来ていることが分かる。

図１８に、遅延推定の正解率を従来技術と比較した結果を示す。縦軸は正解率、横軸は遅延推定方法を表す。遅延推定方法は、左からＦＦＴの点数４０９６点の従来法、ＦＦＴの点数１０２４点の従来法、ＦＦＴの点数１２８点で３５フレーム分は共通で実施例１と実施例２と実施例３と順に並んでいる。

この発明の実施例１で説明した遅延推定方法でも、従来技術のＦＦＴの点数４０９６点の正解率６５％を上回る６９％の正解率が得られた。また、実施例２の遅延推定方法では７５％、実施例３では１００％ぼ正解率が得られた。このようにこの発明の遅延推定方法によれば、周波数分析のデータ数が少なくても精度の良い遅延推定を行うことができる。

〔応用例〕
この発明の遅延推定装置１００，２００，３００は、エコー消去装置に用いることができる。図１９に、この発明の遅延推定装置を用いたエコー消去装置４００の機能構成例を示す。エコー消去装置４００は、遅延推定装置１００と、遅延挿入部９５２と、周波数領域適応フィルタ４０１と、加算回路９０２と、を具備する。遅延推定装置１００は、上記した遅延推定装置２００又は３００であっても良い。

遅延推定装置１００は、受話信号である時間領域の遅延前信号と、送話信号である時間領域の遅延後信号と、を入力としてその両信号間の推定遅延量を求める。遅延挿入部９５２は、遅延推定装置１００の内部の第１の周波数分析部１０で、サンプル点数をＮとした時に、Ｎよりも少ない２Ｎ/Ｍ点の周波数分析によって周波数領域の信号に変換された遅延前信号を、上記遅延推定装置で推定された推定遅延量分の時間遅延させる。

周波数領域適応フィルタ部４０１は、遅延挿入部９５２で遅延が挿入された遅延前信号を入力として擬似エコー信号を生成する。加算部９０２は、遅延推定装置１００の内部の第２の周波数分析部１２で、２Ｎ/Ｍ点の周波数分析によって周波数領域の信号に変換された遅延後信号から擬似エコー信号を除去する。擬似エコー信号が除去された周波数領域の遅延後信号は、２Ｎ/Ｍ点ＩＦＦＴ部４０２において、２Ｎ/Ｍ点の逆ＦＦＴがＫ回実行されて時間領域の送話信号に変換される。

エコー消去装置４００は、エコーキャンセラの為の周波数分析部と、遅延推定の為の周波数分析部と、を共通に用いるようにした構成に特徴がある。この構成にすることで、従来の遅延推定装置を用いたエコー消去装置では、Ｎ点のＦＦＴが２回と、Ｎ点の逆ＦＦＴが１回必要である。この部分の積和演算は１４７４５６回である。

これに対してこの発明のエコー消去装置４００では、信号変換に必要な演算量は、２Ｎ/Ｍ点の逆ＦＦＴをＫ回行うのみである。ここでＮ＝４０９６、Ｍ＝６４、Ｋ＝１６とした条件で、この発明のエコー消去装置４００において追加となるＦＦＴの演算量は、１４３３６回となる。このようにこの発明のエコー消去装置４００は、周波数領域への信号変換処理の演算量を約１/１０に削減することができる。

上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、DVD（Digital Versatile Disc）、DVD-RAM（Random Access Memory）、CD-ROM（Compact Disc Read Only Memory）、CD-R（Recordable）/RW（ReWritable）等を、光磁気記録媒体として、MO（Magneto Optical disc）等を、半導体メモリとしてEEP-ROM（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したDVD、CD-ROM等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (12)

1. 白色化相互相関を求める時間幅の遅延前信号から採取されるサンプル点数をＮとした時に、上記遅延前信号を上記Ｎ点よりも少ない２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理をＭ回単位で繰り返す第１の周波数分析過程と、
   上記遅延前信号を遅延させた遅延後信号を２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理を、上記Ｍフレームに遅延を推定する範囲であるＫフレームを加えたＭ＋Ｋ−１回の単位で繰り返す第２の周波数分析過程と、
   上記第１の周波数分析過程で周波数領域の信号に変換した遅延前信号を上記Ｍ個分保存する第１のバッファ過程と、
   上記第２の周波数分析過程で周波数領域の信号に変換した遅延後信号を、上記Ｍ＋Ｋ−１個分保存する第２のバッファ過程と、
   上記第１のバッファ過程で保存された遅延前信号の共役を取った信号と、上記第２のバッファ過程で保存された遅延後信号とを同一周波数毎に乗算し、当該乗算結果をフレームについて加算し上記Ｋ要素までのクロススペクトルを計算するクロススペクトル計算過程と、
   上記第１のバッファ過程で保存された遅延前信号のノルムの２乗の逆数を白色化係数として計算する白色化係数計算過程と、
   上記クロススペクトルに上記白色化係数を乗算する白色化過程と、
   上記白色化過程で白色化されたクロススペクトルを、時間領域に変換する時間領域変換過程と、
   上記時間領域変換過程の出力を時間シフトしながら第１要素から上記Ｋ要素まで加算して保存する時間シフト加算過程と、
   上記時間シフト加算過程で保存された信号の最大ピーク位置を推定遅延量として求める最大ピーク位置探索過程と、
   を備える遅延推定方法。
2. 請求項１に記載した遅延推定方法において、
   上記白色化係数計算過程は、
   上記第１のバッファ過程で保存された上記遅延前信号のノルムの２乗を計算するノルム計算ステップと、
   上記ノルムの２乗を周波数について加算平均してフルバンド平均パワーを計算するフルバンド平均パワー計算ステップと、
   上記フルバンド平均パワーに、予め設定された定数αを乗じる係数乗算ステップと、
   上記ノルムの２乗に、上記フルバンド平均パワーに定数αを乗じた値を加算する加算ステップと、
   上記加算ステップで出力される出力信号の逆数を計算し、白色化係数として出力する逆数計算ステップと、
   を含むことを特徴とする遅延推定方法。
3. 請求項１又は２に記載した遅延推定方法と、
   上記遅延推定方法に加えて、
   上記遅延前信号から遅延推定に適した区間を設定する推定区間設定過程を備え、
   上記推定区間設定過程で設定した区間の上記遅延前信号を用いて遅延推定を行うことを特徴とする遅延推定方法。
4. 請求項３に記載した遅延推定方法において、
   上記推定区間設定過程は、
   上記遅延前信号の振幅のディップホールドを行って雑音レベルを推定する雑音レベル推定ステップと、
   上記遅延前信号の振幅のピークホールドを行って残響音レベルを推定する残響レベル推定ステップと、
   上記雑音レベルに１以上の定数を乗じた値と、上記残響レベルに１未満の定数を乗じた値と、を比較し大きい方を閾値として設定する閾値設定ステップと、
   上記閾値と上記遅延前信号の振幅とを比較して当該閾値よりも上記遅延前信号の振幅が大きい区間を、遅延推定する区間として検出する区間検出ステップと、
   を含むことを特徴とする遅延推定方法。
5. 請求項１乃至４の何れかに記載した遅延推定方法と、
   上記遅延推定方法に加えて、
   上記第１の周波数分析過程で周波数領域の信号に変換された受話信号を、上記遅延推定方法で推定された推定遅延量分の時間遅延させる遅延挿入過程と、
   上記遅延挿入過程で遅延が挿入された受話信号を入力として擬似エコー信号を生成する周波数領域適応フィルタ過程と、
   上記周波数領域の遅延後信号から上記擬似エコー信号を除去する加算過程と、
   上記擬似エコー信号が除去された周波数領域の遅延後信号を、２Ｎ/Ｍ点の逆ＦＦＴをＫ回実行して時間領域の送話信号に変換する２Ｎ/Ｍ点ＩＦＦＴ過程と、
   を備えたエコー消去方法。
6. 白色化相互相関を求める時間幅の遅延前信号から採取されるサンプル点数をＮとした時に、上記遅延前信号を上記Ｎ点よりも少ない２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理をＭ回単位で繰り返す第１の周波数分析部と、
   上記遅延前信号を遅延させた遅延後信号を２Ｎ/Ｍ点の周波数分析で周波数領域に変換する処理を、上記Ｍフレームに遅延を推定する範囲であるＫフレームを加えたＭ＋Ｋ−１回の単位で繰り返す第２の周波数分析部と、
   上記第１の周波数分析部で周波数領域の信号に変換した遅延前信号を上記Ｍ個分保存する第１のバッファ部と、
   上記第２の周波数分析部で周波数領域の信号に変換した遅延後信号を、上記Ｍ＋Ｋ−１個分保存する第２のバッファ部と、
   上記第１のバッファ部で保存された遅延前信号の共役を取った信号と、上記第２のバッファ部で保存された遅延後信号とを同一周波数毎に乗算し、当該乗算結果をフレームについて加算し上記Ｋ要素までのクロススペクトルを計算するクロススペクトル計算部と、
   上記第１のバッファ部で保存された遅延前信号のノルムの２乗の逆数を白色化係数として計算する白色化係数計算部と、
   上記クロススペクトルに上記白色化係数を乗算する白色化部と、
   上記白色化部で白色化されたクロススペクトルを、時間領域に変換する時間領域変換部と、
   上記時間領域変換部の出力を時間シフトしながら第１要素から上記Ｋ要素まで加算して保存する時間シフト加算部と、
   上記時間シフト加算部で保存された信号の最大ピーク位置を推定遅延量として求める最大ピーク位置探索部と、
   を具備する遅延推定装置。
7. 請求項６に記載した遅延推定装置において、
   上記白色化係数計算部は、
   上記第１のバッファ部で保存された上記遅延前信号のノルムの２乗を計算するノルム計算手段と、
   上記ノルムの２乗を周波数について加算平均してフルバンド平均パワーを計算するフルバンド平均パワー計算手段と、
   上記フルバンド平均パワーに、予め設定された定数αを乗じる係数乗算手段と、
   上記ノルムの２乗に、上記フルバンド平均パワーに定数αを乗じた値を加算する加算手段と、
   上記加算手段が出力する出力信号の逆数を計算し、白色化係数として出力する逆数計算手段と、
   を備えることを特徴とする遅延推定装置。
8. 請求項６又は７に記載した遅延推定装置において、
   更に、上記遅延前信号から遅延推定に適した区間を設定する推定区間設定部を備え、
   上記推定区間設定部で設定した区間の上記遅延前信号を用いて遅延推定を行うことを特徴とする遅延推定装置。
9. 請求項８に記載した遅延推定装置において、
   上記推定区間設定部は、
   上記遅延前信号の振幅のディップホールドを行って雑音レベルを推定する雑音レベル推定手段と、
   上記遅延前信号の振幅のピークホールドを行って残響音レベルを推定する残響レベル推定手段と、
   上記雑音レベルに１以上の定数を乗じた値と、上記残響レベルに１未満の定数を乗じた値と、を比較し大きい方を閾値として設定する閾値設定手段と、
   上記閾値と上記遅延前信号の振幅とを比較して当該閾値よりも上記遅延前信号の振幅が大きい区間を、遅延推定する区間として検出する区間検出手段と、
   を備えることを特徴とする遅延推定装置。
10. 請求項６乃至９の何れかに記載した遅延推定装置と、
    更に、
    上記第１の周波数分析部で周波数領域の信号に変換された受話信号を、上記遅延推定装置で推定された推定遅延量分の時間遅延させる遅延挿入部と、
    上記遅延挿入部で遅延が挿入された受話信号を入力として擬似エコー信号を生成する周波数領域適応フィルタ部と、
    上記周波数領域の遅延後信号から上記擬似エコー信号を除去する加算部と、
    上記擬似エコー信号が除去された周波数領域の遅延後信号を、２Ｎ/Ｍ点の逆ＦＦＴをＫ回実行して時間領域の送話信号に変換する２Ｎ/Ｍ点ＩＦＦＴ部と、
    を具備するエコー消去装置。
11. 請求項６乃至９の何れかに記載した遅延推定装置、又は、請求項１０に記載したエコー消去装置としてコンピュータを機能させるためのプログラム。
12. 請求項１１に記載した何れかのプログラムを記録したコンピュータで読み取り可能な記録媒体。