反射音情報推定装置、反射音情報推定方法、プログラム
【課題】基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する技術を提供する。
【解決手段】Ｑを２以上の予め定められた整数として、Ｑ個の反射音（反射音は、所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（伝達特性関数）に複素振幅を乗じて表される）それぞれの複素振幅（到来振幅）を用いて、基準反射音以外の反射音（対象反射音）それぞれについて、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とする。
【解決手段】Ｑを２以上の予め定められた整数として、Ｑ個の反射音（反射音は、所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（伝達特性関数）に複素振幅を乗じて表される）それぞれの複素振幅（到来振幅）を用いて、基準反射音以外の反射音（対象反射音）それぞれについて、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、反射音に関する情報、特に基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する技術に関する。
【背景技術】
【０００２】
電話や音声会議といった音声情報をやりとりするシステムを一般に音声通信システムと呼ぶ。音声通信システムにおいて、反射音に関する情報（到来振幅、到来方向など）を得ることは非常に重要なことである。会議室のような残響環境下において、マイクロホンを通して収音される収音信号の中には発話者のような音源から直接到来する直接音だけではなく、床、壁や天井に反射して到来する反射音が混在する。したがって、このような残響環境下で或る話者の発言を収録すると、直接音から遅延して反射音が混入するため、聞き取りづらくなってしまう。収音信号から各反射音の到来情報を推定して、反射音を除去することができれば、聞き取りやすい音声に回復することができる。ここで、反射音情報を推定する従来研究として、非特許文献１が挙げられる。
【０００３】
非特許文献１に開示される技術を実現する機能構成を図１に示す。この技術における処理手順は次のとおりである。
【０００４】
１．インパルス音源１００から放射された音源信号を４ｃｈのマイクロホン１１０−１，１１０−２，１１０−３，１１０−４を用いて収音する。ＡＤ変換部１２０は、収音されたアナログ信号をデジタル信号x→(t)＝[x1(t),x2(t),x3(t),x4(t)]Tへ変換する。ここで、[・]Tは転置を表す。ｔは離散時間のインデックスを表す。４本のマイクロホンは正四面体の頂点に配置されていることとする。
【０００５】
２．インパルス応答算出部１３０は、デジタル信号x→(t)＝[x1(t),x2(t),x3(t),x4(t)]Tを入力とし、各マイクロホンのインパルス応答h→(t)＝[h1(t),h2(t),h3(t),h4(t)]Tを算出する。インパルス応答の算出方法には、ＴＳＰ法やＭ系列法等があり、いかなる方法を用いてインパルス応答を算出してもよい。
【０００６】
３．仮想音源算出部１４０は、４ｃｈのインパルス応答h→(t)＝[h1(t),h2(t),h3(t),h4(t)]Tを入力とし、仮想音源情報v→＝[v→1,…,v→D]Tを出力する。Ｄは仮想音源の数を表す。仮想音源とは、各反射音の到来振幅、到来方向、到来時間を表現するために仮想的に存在するとされる音源である。図２を参照して、仮想音源について説明する。図２には、右側の壁で反射した音源信号をマイクロホンで受音する経路が書かれている。右側の壁で反射して到来する音源信号（反射音）は、「仮想音源」と書かれた位置から直接到来する信号と等価である（ただし、壁での反射による減衰や距離減衰の影響は受ける）。
【０００７】
この従来技術の詳細について説明する。インパルス応答を近接した４つの受音点（マイクロホンの位置）で測定すると反射音の到来時刻にわずかな差が生じる。インパルス応答の短い区間の相互相関を利用して、各マイクロホンにおける反射音の対応付けを行うことで、図３のように、ｎ番目の反射波に関するそれぞれの受音点での到来時刻ｔ1n，ｔ2n，ｔ3n，ｔ4n（１≦ｎ≦Ｄ）が求まる。正四面体マイクロホンアレーの辺の長さをｄ、音速をｃとすると、各仮想音源情報vn→＝[Xn,Yn,Zn,Sn]Tが求まる。ここで、Xn,Yn,Znはｎ番目の仮想音源の位置を表し（式（１）−（３）参照）、これは各反射音の到来方向と到来時間に対応する情報を持つ。また、Snはｎ番目の仮想音源の強さを表し、４ｃｈのインパルスで対応付けされたｎ番目の反射音の振幅の平均で求まる。
【数１】
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】山崎芳男ら、「近接する４点のインパルス応答により求めたホールの空間情報」、日本音響学会講演論文集、１９８１年５年、pp.759-760.
【発明の概要】
【発明が解決しようとする課題】
【０００９】
従来技術によると、仮想音源情報と呼んでいた反射音の「到来振幅」、「到来方向」、「到来時刻」を推定するためには、インパルス応答をあらかじめ用意することが必要であった。しかし、インパルス応答を用意するためには特殊な信号を用いて観測する必要があるため、あらゆる位置でのインパルス応答が事前に用意されているという条件は現実的ではない。
【００１０】
そこで本発明は、特殊な信号を用いることなく、反射音情報、特に基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する技術を提供することを目的とする。
【課題を解決するための手段】
【００１１】
反射音を、所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（伝達特性関数）に複素振幅を乗じて表されるとする。Ｑを２以上の予め定められた整数として、Ｑ個の反射音それぞれの複素振幅（到来振幅）を用いて、基準反射音以外の反射音（対象反射音）それぞれについて、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とする。
【００１２】
基準反射音として、Ｑ個の反射音の中から各対象反射音に共通の基準反射音を決定してもよいし、対象反射音ごとにＱ個の反射音の中から基準反射音を決定してもよい。
【００１３】
Ｑ個の反射音の到来振幅はそれぞれ、音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（観測信号）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものであるとしてもよい。
【００１４】
あるいは、Ｑ個の反射音の到来振幅は、相異なるＱ個の位置に対応するＱ個の反射音を観測信号から減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅であるとしてもよい。
【発明の効果】
【００１５】
本発明に拠ると、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とすることから、インパルス応答を求めるために音源信号に特殊な信号を用いることなく、反射音情報（基準反射音に対する他の反射音の相対的な到来時刻差）を推定することが可能である。反射音情報が得られると、従来の音声情報処理技術では実現できなかった音源距離の推定や、音声強調（遠方音の収音や距離別の収音）といった用途に応用できる。
【図面の簡単な説明】
【００１６】
【図１】従来技術における反射音情報推定技術の機能構成を示す図。
【図２】仮想音源を説明するための図。
【図３】従来技術における反射音の対応付けを説明するための図。
【図４】第１実施形態に係る反射音情報推定装置の機能構成を示す図。
【図５】第１実施形態に係る反射音情報推定方法の処理手順を示す図。
【図６】２次元マイクロホンアレーの構成例を示す図。
【図７】ｐ番目の点[xp,yp,zp]とｍ番目の受音点[um,vm,wm]との間の伝達特性を説明するための図。
【図８】図６に示すマイクロホンアレーを用いて観測した或る平面での音圧分布が例えば直接音と反射音１と反射音２との重畳で得られていることを説明するための図。
【図９】本発明の原理を説明するための図。
【図１０】実験例においてテンプレート情報作成のための区間座標とマイクロホンアレーとの位置関係を示す図。
【図１１】実験結果を示す図。
【図１２】実験例により本発明の有効性を説明するための図。
【図１３】第２実施形態に係る反射音情報推定方法の処理手順を示す図。
【図１４】（ａ）理想的には、推定到来方向に関する情報だけが抽出されるべきことを説明するための図。（ｂ）実際には、推定到来方向以外の方向に関する情報も混在してしまうことを説明するための図。
【図１５】残差信号のパワーを全周波数に亘って総括することにより、推定到来方向以外の方向の影響を減らすことを説明するための図。
【図１６】実用レベルの２次元マトリクスマイクロホンアレーを用いた場合における音圧分布とその分解を示す図。
【発明を実施するための形態】
【００１７】
《第１実施形態》
本発明は、発話信号のような音源から放射された音声信号（音源信号）を複数のマイクロホンで構成されるマイクロホンアレーで収音した信号（収音信号）から反射音の「到来振幅」と「到来方向」の少なくともいずれか一つを推定する。第１実施形態の機能構成および処理フローを図４と図５に示す。
【００１８】
音源２００から放射された音源信号をＭｃｈのマイクロホン２１０−１，…，２１０−Ｍを用いて収音する（ステップＳ１）。Ｍは、４より大きい値が望ましい。ＡＤ変換部２２０が、収音されたアナログ信号をデジタル信号xx→(t)＝[xx1(t),…,xxM(t)]Tへ変換する（ステップＳ２）。ここで、[・]Tは転置を表す。ｔは離散時間のインデックスを表す。
【００１９】
Ｍ本のマイクロホンの並べ方は、２次元または３次元的に等間隔で配置することが望ましい。これは、反射音の到来方向とテンプレート（後で説明するが、反射音の伝達特性を模擬したものである）の対応を一意に定めるためである。なお、原理的に、１次元的にマイクロホンを配置しても、あるいは等間隔に配置しなくても、本発明を実施できるが、反射音の伝達特性と反射音の到来方向が一対一の関係とならないため、２次元または３次元的に等間隔で配置することが望ましい。２次元平面上に等間隔にマイクロホンを並べた場合の一例を図６に示す。マイクロホン間隔ｄは、空間サンプリング定理を満たすように設定されていることが望ましい。空間サンプリング定理を満たす場合、マイクロホン間隔ｄは、式（４）を満たす数値となる。ｃは音速であり、ｆは解析対象とする周波数である。例えば、４ｋＨｚの周波数を解析する場合、マイクロホン間隔を４ｃｍ程度に設定するのがよい。
【数２】
【００２０】
フレーム分割部２３０は、ＡＤ変換部２００が出力したデジタル信号xx→(t)＝[xx1(t),…,xxM(t)]Tを入力とし、チャネルごとに複数サンプルから成るデジタル信号の組（フレーム）に分割された信号x→(k)＝[x1(k),…,xM(k)]Tを出力する（ステップＳ３）。ｋはフレーム番号を表すインデックスである。フレーム分割は、各チャネルのデジタル信号xxi(t)（１≦ｉ≦Ｍ）ごとにＷ点分をバッファリングして出力する処理である。Ｗはサンプリング周波数にもよるが、１６ｋＨｚサンプリングの場合には５１２点あたりが妥当である。
【００２１】
周波数領域変換部２４０は、各フレームのデジタル信号x→(k)を入力として、周波数領域の信号X→(ω,k)＝[X1(ω,k),…,XM(ω,k)]Tに変換して出力する（ステップＳ４）。この信号X→(ω,k)を観測信号と呼ぶことにする。ここで、ωは離散周波数のインデックスを指し（周波数ｆと角周波数ωとの間にはω=2πfの関係があるから、周波数のインデックスωをこの角周波数ωと同一視してもかまわない。以下、ωに関して「周波数のインデックス」を単に「周波数」ともいう）、ｋはフレームのインデックスを指す。周波数領域に変換する方法の一つに、離散フーリエ変換があるが、周波数領域に変換するのであれば、他の方法を用いてもよい。周波数領域の観測信号X→(ω,k)は、各周波数ω、フレームｋごとに出力される。
【００２２】
テンプレート生成部２５０は、Ｐ個のテンプレートSp→(ω)の集合（ただし、計算の都合、ベクトル表記している）であるテンプレート情報S→(ω)＝[S1→(ω),…,SP→(ω)]（∀ω∈Ω；Ωは周波数のインデックスωの集合）を周波数ωごとに生成する（ステップＳｐ）。この処理は通常、ステップＳ１−Ｓ４の各処理に先立ち実施される。Ｐはテンプレートの総数を表し、予め２以上の整数値に設定されている。テンプレート総数Ｐは多ければ多いほど高精度な反射音情報の推定に繋がるが、計算量が多くなるので、例えばＰ＝１０００くらいに設定するのが良い。この処理は、マイクロホンで信号を観測する以前にあらかじめ行う処理である。また、マイクロホンの位置（例えばマイクロホンの間隔ｄ）を変更したり、テンプレート総数Ｐを変更したりしない限り、通常、テンプレートを毎回作り直す必要はない。ここで言う“テンプレート”とは、反射音の到来方向に対応する伝達特性（音響伝播特性）を模擬したものである。ｐ番目（１≦ｐ≦Ｐ）のテンプレートSp→(ω)＝[Sp1(ω),…,SpM(ω)]T（ω∈Ω）は、予め定められたｐ番目の点[xp,yp,zp]とＭ個の受音点（ここで受音点はマイクロホンが配置される位置であり、ｍ番目（１≦ｍ≦Ｍ）の受音点を[um,vm,wm]とする）の間の周波数ごとの伝達特性を表す（図７を参照）。ｐ番目のテンプレートSp→(ω)の各要素Spm(ω)の算出式の一例を式（５）に示す。記号ｉは虚数単位を表す。
【数３】
【００２３】
ｐ番目のテンプレートSp→(ω)には方向情報θp→(ω)が対応付けられている。方向情報θp→(ω)は、ｐ番目の点[xp,yp,zp]および受音点[um,vm,wm]の位置座標の基準となる３次元直交座標系の原点からｐ番目の点[xp,yp,zp]を見た方向であり、例えば（当該３次元直交座標系の原点と共通の原点を持つ）球座標系における二つの偏角（極角θp,polと方位角θp,azi）として表される。つまり、θp→(ω)＝[θp,pol(ω),θp,azi(ω)]である。なお、ｐ番目のテンプレートSp→(ω)にｐ番目の点[xp,yp,zp]が関連付けられていれば方向情報θp→(ω)は位置[xp,yp,zp]から計算可能であるから、ｐ番目のテンプレートSp→(ω)に方向情報θp→(ω)が対応付けられていることは必須要件ではない。なお、３次元直交座標系と球座標系とは相互に転換可能であるから（座標変換）、式（５）の右辺は位置[x,y,z]でなく方向情報θp→(ω)＝[θp,pol(ω),θp,azi(ω)]を用いて、例えば式（５ａ）のように表すこともできる。ここで、ｄはマイクロホン間隔であり、マイクロホンアレーをΦ行Ξ列（Φ×Ξ＝Ｍ）の２次元マイクロホンアレーとし、ｍ番目のマイクロホンの位置をφ行ξ列（１≦φ≦Φ，１≦ξ≦Ξ）にあるとする。
【数４】
【００２４】
また、第１実施形態のようにテンプレートが方向に対応している場合、Ｐ個の点[xp,yp,zp]（１≦ｐ≦Ｐ）の位置は互いに方向の異なる位置であることが好ましく、例えば各点[xp,yp,zp]が原点から十分に離れた等距離にあるとして、上記原点を中心とする球面上の異なるＰ個の点とすればよい。各点[xp,yp,zp]を原点から十分に離れた位置とする理由は、音源ないし仮想音源から放射された信号は球面的に伝達するが音源ないし仮想音源から十分に離れた位置（原点）での局所領域では直接音ないし反射音を平面波として模擬できるからである。ただし、テンプレート情報が同じ方向の位置に対応するテンプレートを含むことを排除する趣旨ではない。なお、マイクロホンアレーは上記座標系の原点の近傍（局所領域）に配置されているとする。
【００２５】
テンプレート記憶部２６０は、テンプレート生成部２５０が出力したテンプレート情報S→(ω)を記憶し、解析時に反射音情報推定部２７０にテンプレート情報S→(ω)を提供する役割を果たす。
【００２６】
反射音情報推定部２７０は、周波数領域の観測信号X→(ω,k)とテンプレート情報S→(ω)を入力として、Ｑ個の反射音情報成分rsq→(ω,k)の集合（ただし、計算の都合、ベクトル表記している）である反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを各フレームｋについて周波数ωごとに出力する（ステップＳ５）。ここで、Ｑは推定される反射音の総数を表し、予め１以上の整数値に設定されている。ｑ番目（１≦ｑ≦Ｑ）の反射音情報成分rsq→(ω,k)は、rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]の２要素から成り、rsAq(ω,k)はｑ番目の反射音の到来振幅であり、rsBq(ω,k)はｑ番目の反射音の到来方向である。
【００２７】
反射音情報を推定する原理について説明する。図６に示すような２次元マイクロホンアレーを用いて観測した或る平面での音圧分布の一例を図８の左端の濃淡図として示す。濃淡図として示された音圧分布の見方について、黒い部分は音圧が小さく、白い部分は音圧が大きいことを示す。観測した音圧分布には直接音の音圧分布だけではなく、反射音の音圧分布も混入している。直接音や反射音が十分に遠方より到来する場合において、２次元平面上でのそれぞれの音圧分布は、図８の右側の３つの濃淡図のように縞模様となる。縞模様の「濃淡」が直接音ないし反射音の到来振幅、「回転・周期」が直接音ないし反射音の到来方向にそれぞれ対応する。図８の例では、到来振幅や到来方向が異なる直接音、反射音１、反射音２の各音圧分布の重畳で観測信号の音圧分布が構成されることを示している。周波数領域で考えると、直接音や各反射音は到来方向に応じて周波数の変化する複素正弦波で表され、観測信号は直接音と各反射音に対応する複数の複素正弦波が重畳したものとして表される。ところで、本発明で解決する問題は、観測信号のみを用いて、反射音の到来振幅および／または到来方向を推定することである。この課題解決は、図８の左端に描かれた音圧分布から図８の右側の３つの濃淡図の直接音や各反射音に対応する縞模様の「濃淡」および／または「回転・周期」を推定することに対応する。
【００２８】
図９を参照して、反射音情報rs→(ω,k)を推定する手法の概略について説明する。ある２次元平面で観測した観測信号に含まれているパワーの最も強い反射音０（q=1に相当し、最も強いパワーを持つことから、通常、この反射音０は「直接音」として理解される）を推定し、観測信号から当該反射音０を減算して残差信号E2を得る。次に当該残差信号E2に含まれているパワーの最も強い反射音１（q=2に相当する）を推定し、当該残差信号E2から当該反射音１を減算して新たな残差信号E3を得る。次に、当該残差信号E3に含まれているパワーの最も強い反射音２（q=Q=3に相当する）を推定する。ここでは、Ｑ＝３の場合を説明したが、一般的に、ｑ番目の残差信号Eq（ただし１番目の残差信号は観測信号とする）に含まれているパワーの最も強いｑ番目の反射音q-1（ただし反射音０は直接音である）を減算する操作をｑ＝Ｑまで逐次実行することでＱ個の反射音情報成分（rs1→(ω,k),…,rsQ→(ω,k)）を得る。１番目の反射音情報成分rs1→(ω,k)は反射音０（直接音）に対応し、２番目の反射音情報成分rs2→(ω,k)は反射音１に対応し、３番目の反射音情報成分rs3→(ω,k)は反射音２に対応し、・・・、Ｑ番目の反射音情報成分rsQ→(ω,k)は反射音Q-1に対応する。Ｑは、計算パワーや反射音情報を用いるアプリケーションにも依存するが、３０くらいに設定するのが良い。
【００２９】
なお、図８および図９の音圧分布はそれぞれ高解像度の濃淡図として示されているが、このような高解像度の濃淡図として音圧分布を示すためには極めて多くのマイクロホンを必要とし、実用的ではない。他方、実用レベルの２次元マトリクスマイクロホンアレーとして例えば１００個のマイクロホンを１０×１０の２次元マトリクスマイクロホンアレーとして用いた場合でさえ、粗い（低解像度）濃淡図（図１６参照）として示される音圧分布しか得られない。そこで、実用の観点から、低解像度の音圧分布しか得られないような状況の下で、精度良く反射音の到来振幅や到来方向を推定することが求められる。本発明では、空間分解能の向上のために任意の位置から到来する平面波を具体的に表現することとし（定式化）、パワーが大きな反射音の影響を受けてパワーの小さいな反射音を推定できなくなることを防止するために、既に推定された反射音を観測信号から除去して次の反射音を推定する（分解）。定式化についてはテンプレート情報として説明したとおりであり、分解については反射音情報rs→(ω,k)の推定手法の概略で説明したとおりである。
【００３０】
上で述べた反射音情報rs→(ω,k)を推定する手法について詳細を述べる。説明に先立ち、記号の定義を行う。ｑ番目の残差信号をEq→(ω,k)=[Eq1(ω,k),…,EqM(ω,k)]T、推定されたｑ番目の反射音（q=1の場合は直接音を表す）をAq(ω,k,g(q,k))Sg(q,k)→(ω)とする。ここで、g(q,k)は、フレームｋにおけるｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(q,k))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(q,k)→(ω)と反射音との相違を表す。ｑの昇順に残差信号から反射音を減算する上述の方法を式で表すと式（６）のようになる。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数５】
【００３１】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（７）、式（８）で与えられる。方向情報θg(q,k)→(ω,k)はテンプレートSg(q,k)→(ω)に対応する方向情報である。
【数６】
【００３２】
次に、ｑ番目の反射音Aq(ω,k,g(q,k))Sg(q,k)→(ω)を推定する手法について述べる。ｑ番目の反射音Aq(ω,k,g(q,k))Sg(q,k)→(ω)は、ｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするものとして推定される。記号Ｈは共役転置を表す。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(q,k))とSg(q,k)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の＜処理１＞と＜処理２＞はｑの昇順に各ｑについて行われる。
【００３３】
＜処理１＞
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（１０）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(1,k),…,g(q-1,k)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目のテンプレートSp→(ω)が残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（９）により求められる。なお、この段階では、式（９）左辺のｑは意味を持たないことに留意されたい。
【数７】
【００３４】
＜処理２＞
集合Λの要素の個数（濃度）を|Λ|とすると、式（９）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）を用いて、テンプレートSg(q,k)→(ω)のインデックスを表すg(q,k)は、残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするインデックスとして式（１０）により得られる。
【数８】
【００３５】
ｑ番目の反射音を構成する係数Aq(ω,k,g(q,k))は、式（１０）で得られたg(q,k)を式（９）のｐに代入して得られる。式（１０）の計算の際に係数Aq(ω,k,g(q,k))は計算されているので、メモリに余裕がある場合には、式（１０）の際に得られた係数Aq(ω,k,g(q,k))をメモリに格納しておき、式（１０）でg(q,k)が得られた後に、メモリからp=g(q,k)に対応する係数Aq(ω,k,g(q,k))を呼び出せばよい。
【００３６】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。
【００３７】
到来時刻差推定部２８０は、反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを入力として（正確には、複素振幅{rsA1(ω,k),…,rsAQ(ω,k)}のみで十分である）、基準となる反射音（基準反射音）に対する各反射音の相対的な到来時刻差rsDq(k)（１≦ｑ≦Ｑ）の集合である到来時刻差情報rsD→(k)＝[rsD1(k),…,rsDQ(k)]T（ただし、到来時刻差を用いる任意のアプリケーションの便宜上、ベクトル表記している）を、各フレームｋについて出力する（ステップＳ６）。具体的には、後述の((処理１))と((処理２))を各ｑ（１≦ｑ≦Ｑ）について行うことにより、rsDq(k)（１≦ｑ≦Ｑ）が求められる。ここでは処理対象のｑを特にｑｃと記す。
【００３８】
((処理１))
推定されたＱ個の反射音から基準反射音のインデックスＮを決める。((処理１))を各ｑ（１≦ｑ≦Ｑ）について行う場合、Ｎ＝ｑｃならば基準反射音に対するｑｃ番目の反射音の到来時刻差が０であることは自明であるから、インデックスＮは反射音を識別するインデックスの差集合{１，…，Ｑ}−{ｑｃ}の中から選択されることが好ましい。基準反射音の決定方法は種々考えられるが、例えば、差集合{１，…，Ｑ}−{ｑｃ}に対応する反射音のうちパワーが強い反射音を選択することが好適である。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して到来振幅の大きさの総和を計算することにより、周波数による影響をキャンセルすることができる。この場合、Ｎは式（１０ａ）により決定される。なお、解析する周波数帯域に含まれる周波数のインデックスωの集合をΩとし、例えば音声信号を扱うのであれば、1.0〜3.0kHz帯域に対応するインデックスの集合をΩとすればよい。
【数９】
【００３９】
((処理２))
Ｎ番目の反射音の到来振幅rsAN(ω,k)とｑｃ番目（ｑｃ≠Ｎ）の反射音の到来振幅rsAqc(ω,k)との位相差を用いて相対的な到来時刻差rsDqc(k)を算出する。一般に、到来振幅rsAq(ω,k)は複素振幅であるから振幅λq(ω,k)と時刻τq(ω,k)を用いてrsAq(ω,k)=λq(ω,k)exp[iωτq(ω,k)]と表すことができる。従って、Ｎ番目の反射音に対するｑ番目（ｑ≠Ｎ）の反射音の位相差は、Ｎ番目の反射音の到来振幅rsAN(ω,k)に対するｑ番目（ｑ≠Ｎ）の反射音の到来振幅rsAq(ω,k)の偏角arg(rsAq(ω,k)／rsAN(ω,k))として与えられる。偏角arg(rsAq(ω,k)／rsAN(ω,k))を周波数のインデックスωで除することによりＮ番目の反射音とｑ番目（ｑ≠Ｎ）の反射音との時刻差λq(ω,k)-λN(ω,k)を求めることができる。つまり、λq(ω,k)-λN(ω,k)=arg(rsAq(ω,k)／rsAN(ω,k))／ωである。時刻λN(ω,k)に対する時刻λq(ω,k)が「早い」か「遅い」かの区別を行うため負号を導入する。つまり、-arg(rsAq(ω,k)／rsAN(ω,k))／ω＜０であれば時刻λN(ω,k)に対する時刻λq(ω,k)が「早い」ことを表し、-arg(rsAq(ω,k)／rsAN(ω,k))／ω＞０であれば時刻λN(ω,k)に対する時刻λq(ω,k)が「遅い」ことを表す。-arg(rsAq(ω,k)／rsAN(ω,k))／ω=０の場合、時刻差が無いことを表す。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して-arg(rsAq(ω,k)／rsAN(ω,k))／ωの総和を計算することにより、周波数による影響をキャンセルすることができる。即ち、相対的な到来時刻差rsDqc(k)は、式（１０ｂ）によって与えられる。集合Ωの要素の個数（濃度）を|Ω|とする。
【数１０】
【００４０】
なお、（(処理１)）にて基準反射音のインデックスを決定するが、基準反射音を各ｑについて共通とする場合には、（(処理１)）はひとたび行われれば十分である。この場合、先ず((処理１ａ))を行い、次いで差集合{１，…，Ｑ}−{Ｎ}の各要素ｑｃについて上述の((処理２))を行えばよい。
【００４１】
((処理１ａ))
推定されたＱ個の反射音から基準反射音のインデックスＮを決める。基準反射音の決定方法は種々考えられるが、例えば、Ｑ個の反射音のうちパワーが強い反射音（直接音）を選択することが好適である。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して到来振幅の大きさの総和を計算することにより、周波数による影響をキャンセルすることができる。この場合、Ｎは式（１０ｃ）により決定される。
【数１１】
【００４２】
なお、マイクロホンアレーから音源ないし仮想音源までの音源距離disq(ω,k)は、相対的な到来時刻差を定数倍した量であるから、式（１０ｄ）によって与えられる。ｃは音速である。
【数１２】
【００４３】
以下、相対的な到来時刻差を求めるに必要な反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを求めるいくつかの実施形態を説明する。いずれの実施形態であっても、第１実施形態のステップＳ６の処理を行うことで、相対的な到来時刻差rsDq(k)（１≦ｑ≦Ｑ）を求めることができる。
【００４４】
《第２実施形態》
第１実施形態ではテンプレート情報S→(ω)を用いて反射音情報rs→(ω,k)を求めたが、Ｐ個のテンプレートSp→(ω)の集合であるテンプレート情報S→(ω)を事前に求めておくことは必ずしも必須ではない。テンプレート情報S→(ω)を事前に求めておかない実施形態を第２実施形態として説明する。
【００４５】
第２実施形態では、第１実施形態におけるステップＳ１−Ｓ４の各処理が実施されるが、第１実施形態におけるステップＳｐの処理が不要であり、さらに第１実施形態のステップＳ５の処理に替えてステップＳ５ａの処理が行われる（図１３参照）。そこで、第１実施形態と同じ事項については重複説明を省略し、第１実施形態と異なる事項について説明する。
【００４６】
記号の定義を行う。ｑ番目の残差信号をEq→(ω,k)=[Eq1(ω,k),…,EqM(ω,k)]T、ｑ番目の反射音（q=1の場合は直接音を表す）をAq(ω,k)Rq→(ω,θq→(ω,k))とする。反射音を構成するRq→(ω,θq→(ω,k))=[R1(ω,θq→(ω,k)),…,RM(ω,θq→(ω,k))]Tは、空間中の任意の位置[x,y,z]と各マイクロホンとの間の周波数ごとの伝達特性を模擬した関数（以下、伝達特性関数という）であり、各マイクロホンに対する伝達特性を模擬した関数であれば何でもよい。通常、伝達特性関数を構成する各伝達特性Rm(ω,θq→(ω,k))とテンプレートの各要素Spm(ω)の算出式とは同じである。この場合、方向情報θq→(ω,k)で表される方向に在る位置[x,y,z]とｍ番目の受音点[um,vm,wm]との間の周波数ごとの伝達特性Rm(ω,θq→(ω,k))は式（１１）で表される。なお、方向情報θq→(ω,k)で表される方向に在る位置[x,y,z]は、例えば、上記座標系原点から十分に離れた球面上の位置とすればよい。位置[x,y,z]を原点から十分に離れた位置とする理由は既述のとおりであり、詳しくは位置[x,y,z]はマイクロホンアレーが配置されている局所領域にて音源ないし仮想音源からの直接音ないし反射音を平面波として模擬できる距離にある空間中の任意の位置であることが好ましい。なお、３次元直交座標系と球座標系とは相互に転換可能であるから（座標変換）、式（１１）の右辺は位置[x,y,z]でなく方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]を用いて、例えば式（１１ａ）のように表すこともできる。ここで、ｄはマイクロホン間隔であり、マイクロホンアレーをΦ行Ξ列（Φ×Ξ＝Ｍ）の２次元マイクロホンアレーとし、ｍ番目のマイクロホンの位置をφ行ξ列（１≦φ≦Φ，１≦ξ≦Ξ）にあるとする。
【数１３】
【００４７】
反射音を構成するAq(ω,k)は、音源２００自身が持つ位相や壁での反射、距離による減衰といったテンプレートRq→(ω,θq→(ω,k))と反射音との相違を表し、到来振幅に相当する。ｑの昇順に残差信号から反射音を減算する上述の方法を式で表すと式（１２）のようになる。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数１４】
【００４８】
次に、反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する方法について説明する。
ｑ番目の最適化された反射音Aq(ω,k)Rq→(ω,θq→(ω,k))は、式（１２）で表されるｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って決定される。具体的には、伝達特性関数Rq→(ω,θq→(ω,k))が方向情報θq→(ω,k)で決定されることに注意すると、ｑ番目の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を表現するパラメータAq(ω,k)，θq→(ω,k)の最適値Aq,opt(ω,k)，θq,opt→(ω,k)は式（１３）によって得られる。なお、記号Ｈは共役転置を表す。
【数１５】
【００４９】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（１４）、式（１５）で与えられる。
【数１６】
【００５０】
式（１３）の具体的な計算方法は種々考えられるが、ここではその一例を示す。方向情報θq→(ω,k)を、例えばビームフォーマ法などの到来方位推定方法によって定める。ビームフォーマ法は、指向性ビームを空間走査し、得られた電力スペクトルから電力が大きくなる方向を探索する方法である。ここでは、ビームフォーマ法によりＰ個の到来方向が推定できたとする。
【００５１】
実際には、ビームフォーマ法によって得られる電力スペクトルは到来方向に対して急峻でないことがあり、このような場合、例えば、予め定めたスペクトル強度以上のスペクトル強度を示す電力スペクトルに対応する方向の範囲にて予め定めた間隔で到来方向を定めればよい。具体例として、極角５°、方位角１０°から２０°の範囲で予め定めたスペクトル強度以上のスペクトル強度を示す電力スペクトルが得られたとすると、予め定めた間隔２°ごとに到来方向を定めるとして、（極角５°，方位角１０°），（極角５°，方位角１２°），（極角５°，方位角１４°），（極角５°，方位角１６°），（極角５°，方位角１８°），（極角５°，方位角２０°）を到来方向とすればよい。
【００５２】
また、電力スペクトルが或る方向にて急峻なピークを示したとしても、単純に当該方向を到来方向の一つとして定めるのではなく、当該方向の所定の範囲で到来方向を定めてもよい。具体例として、極角３０°、方位角５０°で急峻なピークを示す電力スペクトルが得られたとすると、所定の範囲（極角±４°，方位角±４°，間隔２°）で到来方向を定めるとして、（極角２６°，方位角４６°），（極角２８°，方位角４６°），（極角３０°，方位角４６°），（極角３２°，方位角４６°），（極角３４°，方位角４６°），（極角２６°，方位角４８°），（極角２８°，方位角４８°），（極角３０°，方位角４８°），（極角３２°，方位角４８°），（極角３４°，方位角４８°），（極角２６°，方位角５０°），（極角２８°，方位角５０°），（極角３０°，方位角５０°），（極角３２°，方位角５０°），（極角３４°，方位角５０°），（極角２６°，方位角５２°），（極角２８°，方位角５２°），（極角３０°，方位角５２°），（極角３２°，方位角５２°），（極角３４°，方位角５２°），（極角２６°，方位角５４°），（極角２８°，方位角５４°），（極角３０°，方位角５４°），（極角３２°，方位角５４°），（極角３４°，方位角５４°）を到来方向とすればよい。なお、第１実施形態ではＰは固定値であったが、第２実施形態ではＰはビームフォーマ法などの到来方位推定方法による推定成果に依存する値であることに留意されたい。
【００５３】
ビームフォーマ法によって得られたＰ個の到来方向についてテンプレートを生成する。テンプレートの各要素の算出式は例えば式（６）である。このＰ個のテンプレート（テンプレート情報S→(ω)）を用いて第１実施形態で説明した処理を行えばよい。
【００５４】
《第３実施形態》
第３実施形態は、第２実施形態と異なり、テンプレート情報を利活用して反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する。式（１３）の具体的な計算方法の一例を示す。下記に説明する最適化方法はｑの昇順に各ｑに対して適用される。
【００５５】
§１ 方向情報の初期値設定
最初に、方向情報θq→(ω,k)の初期値θini,q→(ω,k)をテンプレート情報S→(ω)を用いて決定する。このために、推定されるべき到来方向に最も近いと考えられる方向に対応するテンプレートを決定し、この決定されたテンプレートに対応する方向情報を方向情報θq→(ω,k)の初期値θini,q→(ω,k)とすればよい。
【００５６】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、ｑ＋１番目の残差信号Eq+1→(ω,k)は式（１６）のように表される。ただしE1→(ω,k)=X→(ω,k)である。
【数１７】
【００５７】
反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（１６）に基づくｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の=処理１=と=処理２=はｑの昇順に各ｑについて行われる。
【００５８】
=処理１=
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（１８）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目のテンプレートSp→(ω)が残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（１７）により求められる。なお、この段階では、式（１７）左辺のｑは意味を持たないことに留意されたい。
【数１８】
【００５９】
=処理２=
集合Λの要素の個数（濃度）を|Λ|とすると、式（１７）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）を用いて、テンプレートSg(ω,q)→(ω)のインデックスを表すg(ω,q)は、残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするインデックスとして式（１８）により得られる。
【数１９】
【００６０】
従って、方向情報θq→(ω,k)の初期値θini,q→(ω,k)は、式（１８）により得られるg(ω,q)をインデックスに持つテンプレートSg(ω,q)→(ω)に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００６１】
§２ 反射音の最適化
次に、方向情報θq→(ω,k)の初期値θini,q→(ω,k)を起点として、式（１２）で表されるｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするように、反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する。反射音は、係数Aq(ω,k)とRq→(ω,θq→(ω,k))の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。この最適化方法は様々あるが、そのうちの一つの方法（勾配法）について述べる。例示する方法では、方向情報θq→(ω,k)の補正と係数Aq(ω,k)の補正が交互に所定回数（δ回）反復して行われることにより反射音Aq(ω,k)Rq→(ω,θq→(ω,k))が最適化される。δは例えば５０程度の値とされるが１でもよい。
【００６２】
§２．１ 方向情報の補正
方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（１９）による更新によって行われる。初めて§２．１の処理を行う場合、式（１９）右辺の方向情報θq→(ω,k)は§１の処理で得られた初期値θini,q→(ω,k)であり、§２．１の処理が初めてではない場合、式（１９）右辺の方向情報θq→(ω,k)は直前の§２．１の処理で得られた方向情報とする。また、初めて§２．１の処理を行う場合、パワー(Eq+1→(ω,k))HEq+1→(ω,k)の計算に用いられる係数Aq(ω,k)は式（１７）で得られたAq(ω,k,p)とし、§２．１の処理が初めてではない場合、パワー(Eq+1→(ω,k))HEq+1→(ω,k)の計算に用いられる係数Aq(ω,k)は直前の§２．２の処理（後述する）で得られた係数Aq(ω,k)とする。ステップ幅α1およびα2は小さい正の定数であり、収束速度などを考慮して決定されるが、例えばそれぞれ0.1程度の値とされる。
【数２０】
【００６３】
§２．２ 係数の補正
係数Aq(ω,k)の補正は、最小二乗法に基づき、式（２０）に従って新たな係数Aq(ω,k)を求めることにより行われる。式（２０）で用いるRq→(ω,θq→(ω,k))は§２．１の処理で得られた方向情報θq→(ω,k)と式（６）から得られる。
【数２１】
【００６４】
δ回の反復処理が終了した時点で得られている係数Aq(ω,k)と方向情報θq→(ω,k)がAq,opt(ω,k)とθq,opt→(ω,k)であり、ｑ番目の反射音情報成分rsq→(ω,k)となる。すなわち、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（２１）、式（２２）で与えられる。
【数２２】
【００６５】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。なお、δ＝１に設定されている場合、係数の補正を行わないことにより、反射音情報として到来方向のみを求めることができる。
【００６６】
《第４実施形態》
次に、第１、第２、第３実施形態と異なり、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を一括して最適化する方法について説明する。
【００６７】
Ｑ個の最適化された反射音Aq(ω,k)Rq→(ω,θq→(ω,k))は、観測信号からＱ個の反射音を除去して得られる残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)（式（２３ａ）参照）を最小とする基準に従って決定される。具体的には、伝達特性関数Rq→(ω,θq→(ω,k))が方向情報θq→(ω,k)で決定されることに注意すると、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）を表現する各パラメータAq(ω,k)，θq→(ω,k)（１≦ｑ≦Ｑ）の最適値Aq,opt(ω,k)，θq,opt→(ω,k)（１≦ｑ≦Ｑ）は式（２３ｂ）によって得られる。なお、記号Ｈは共役転置を表す。式（２３ｂ）にて、｛(Aq,opt(ω,k)，θq,opt→(ω,k))｝q∈{1,…,Q}は、｛(A1,opt(ω,k)，θ1,opt→(ω,k)),…,(Aq,opt(ω,k)，θq,opt→(ω,k)),…,(AQ,opt(ω,k)，θQ,opt→(ω,k))｝を表し、｛(Aq(ω,k)，θq→(ω,k))｝q∈{1,…,Q}は、｛(A1(ω,k)，θ1→(ω,k)),…,(Aq(ω,k)，θq→(ω,k)),…,(AQ(ω,k)，θQ→(ω,k))｝を表す。
【数２３】
【００６８】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（２４）、式（２５）で与えられる。
【数２４】
【００６９】
式（２３ｂ）の具体的な計算方法は種々考えられるが、ここではその一例を示す。
【００７０】
§１．１ 方向情報の初期値設定
最初に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を設定する。Ｑ個の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定方法として、この実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を、観測信号X→(ω,k)とテンプレート情報S→(ω)を用いて決定するビーム探索的決定方法を説明する。この方法によると、Ｑ個の推定されるべき到来方向それぞれに最も近いと考えられるＱ個の方向に対応するＱ個のテンプレートを決定し、この決定されたＱ個のテンプレートに対応する方向情報をＱ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）とすればよい。なお、この場合、ＰとＱとの間にＱ＜Ｐなる関係がある。
【００７１】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、ｑ番目の反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、観測信号X→(ω,k)からｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)を除去して得られる残差信号Eq→(ω,k)は式（２６）のように表される。
【数２５】
【００７２】
ｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（２６）に基づく残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。
【００７３】
反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。まず、ｐ番目のテンプレートSp→(ω)が残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（２７）により求められる。なお、この段階では、式（２７）左辺のｑは意味を持たないことに留意されたい。
【数２６】
【００７４】
次に、式（２７）に基づき得られたＰ個の係数Aq(ω,k,p)（１≦ｐ≦Ｐ）の中から、その絶対値の大きい方から順にＱ個の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（２８）参照）。記号Λは、インデックスｐの全体の集合{1,…,p,…,P}から式（２８）により決定されたインデックスの集合を除いた集合であり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}である。
【数２７】
【００７５】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（２８）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００７６】
§１．２ 係数Aq(ω,k)の初期値設定
次に、Ｑ個の係数Aq(ω,k)の初期値Aini,q(ω,k)を設定する。Ｑ個の初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）の決定方法として種々のものが考えられるが、ここでは一例として、パワー最小化基準でＱ個の初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）を決定する方法を説明する。まず、Aq(ω,k,p)=0（１≦ｐ≦Ｐ）とする。そして、初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）は、残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)を最小化するように、最小二乗法に基づき、式（２９）により求められる。式（２９）にてFq→(ω,k)は式（３０）で与えられる。式（３０）にてΥ={1,…,q-1,q+1,…,Q}であり、Fq→(ω,k)は観測信号からｑ番目の反射音を除去した残差信号である。なお、Ｑ個の方向情報θq→(ω,k)として§１．１で決定されたＱ個の初期値θini,q→(ω,k)を用いる。式（２９）で用いるRq→(ω,θq→(ω,k))は方向情報の初期値θini,q→(ω,k)と式（６）から得られる。
【数２８】
【００７７】
§２ 反射音の最適化
次に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を起点として、式（２３ａ）で表される残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)を最小とするように、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）を一括して最適化する。各反射音は、係数Aq(ω,k)とRq→(ω,θq→(ω,k))の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。この最適化方法は様々あるが、そのうちの一つの方法（勾配法）について述べる。例示する方法では、方向情報θq→(ω,k)の補正と係数Aq(ω,k)の補正が交互に所定回数（δ回）反復して行われることによりＱ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）が最適化される。δは例えば１００程度の値とされるが１でもよい。
【００７８】
§２．１ 方向情報の補正
Ｑ個の方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]（１≦ｑ≦Ｑ）の補正は、式（３１）による更新によって行われる。各ｑ（１≦ｑ≦Ｑ）について、初めて§２．１の処理を行う場合、式（３１）右辺の方向情報θq→(ω,k)は§１．１の処理で得られた初期値θini,q→(ω,k)であり、§２．１の処理が初めてではない場合、式（３１）右辺の方向情報θq→(ω,k)は直前の§２．１の処理で得られた方向情報とする。また、初めて§２．１の処理を行う場合、パワー(Fq→(ω,k))HFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として式（２９）で得られたAini,q(ω,k)（１≦ｑ≦Ｑ）を用い、§２．１の処理が初めてではない場合、パワー(Fq→(ω,k))HFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として直前の§２．２の処理（後述する）で得られた係数Aq(ω,k)（１≦ｑ≦Ｑ）を用いる。ステップ幅α1およびα2は小さい正の定数であり、収束速度などを考慮して決定されるが、例えばそれぞれ0.1程度の値とされる。
【数２９】
【００７９】
§２．２ 係数の補正
Ｑ個の係数Aq(ω,k)（１≦ｑ≦Ｑ）の補正は、最小二乗法に基づき、式（３２）に従って新たな係数Aq(ω,k)（１≦ｑ≦Ｑ）を求めることにより行われる。式（３２）で用いるRq→(ω,θq→(ω,k))は§２．１の処理で得られた方向情報θq→(ω,k)と式（６）から得られる。式（３２）にてFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として、初めて§２．２の処理を行う場合、§１．１の処理で得られた初期値Aini,q(ω,k)を用い、§２．２の処理が初めてではない場合、直前の§２．２の処理で得られた係数Aq(ω,k)（１≦ｑ≦Ｑ）を用いる。
【数３０】
【００８０】
δ回の反復処理が終了した時点で得られている係数Aq(ω,k)と方向情報θq→(ω,k)のＱ個の組み合わせ(Aq(ω,k)，θq→(ω,k))（１≦ｑ≦Ｑ）が｛(Aq,opt(ω,k)，θq,opt→(ω,k))｝q∈{1,…,Q}であり、Ｑ個の反射音情報成分rsq→(ω,k)（１≦ｑ≦Ｑ）となる。すなわち、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（３３）、式（３４）で与えられる。
【数３１】
【００８１】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。なお、δ＝１に設定されている場合、係数の補正を行わないことにより、反射音情報として到来方向のみを求めることができる。
【００８２】
《第５実施形態》
第５実施形態は、「§１．１ 方向情報の初期値設定」が第４実施形態と異なる。そこで、第４実施形態と同じ事項については重複説明を省略し、第４実施形態と異なる事項について説明する。
【００８３】
§１．１ 方向情報の初期値設定
最初に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を設定する。Ｑ個の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定方法として、この実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を、観測信号X→(ω,k)とテンプレート情報S→(ω)を用いて決定する一般化調和解析的決定方法を説明する。この方法によると、Ｑ個の推定されるべき到来方向それぞれに最も近いと考えられるＱ個の方向に対応するＱ個のテンプレートを決定し、この決定されたＱ個のテンプレートに対応する方向情報をＱ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）とすればよい。なお、この場合、ＰとＱとの間にＱ＜Ｐなる関係がある。
【００８４】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、ｑ番目の反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、観測信号から１番目からｑ番目までのｑ個の反射音を除去して得られる残差信号Eq+1→(ω,k)は式（３５）のように表される。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数３２】
【００８５】
ｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（３５）に基づく残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の＜処理１＞と＜処理２＞はｑの昇順に各ｑについて行われる。
【００８６】
＜処理１＞
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（３７）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目（ｐ∈Λ）のテンプレートSp→(ω)が式（３５）に基づく残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（３６）により求められる。なお、この段階では、式（３６）左辺のｑは意味を持たないことに留意されたい。
【数３３】
【００８７】
＜処理２＞
集合Λの要素の個数（濃度）を|Λ|とすると、式（３６）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）の中から、その絶対値が最大の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（３７）参照）。
【数３４】
【００８８】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（３７）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００８９】
第５実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定に一般化調和解析的決定方法を用いた。一般化調和解析的決定方法によると、初期値計算量がビーム探索的決定方法に比べて増大するが、推定されるべき反射音（正解）に近い初期値を設定できる可能性が高く、この場合、推定精度の向上を望めるだけでなく、反復処理の回数を減らすことができる。
【００９０】
《第６実施形態》
第６実施形態では、「§１．１ 方向情報の初期値設定」が第４実施形態と異なる。そこで、第４実施形態と同じ事項については重複説明を省略し、第４実施形態と異なる事項について説明する。第６実施形態のコンセプトは、第４実施形態よりも簡便でありながら、上述の初期値の偏りを防止することにある。
【００９１】
§１．１ 方向情報の初期値設定
まず、方向情報を構成するθpol(ω)とθazi(ω)のいずれか一方について、偏りが無いように複数の方向を決定する。この例では、極角θpol(ω)について偏りが無いようにβ個（β≧2）の方向{θ1,pol(ω),…,θβ,pol(ω)}を決定する。通常、極角は0°≦θpol(ω)≦180°を満たすから、例えば等間隔10°ごとに方向を定めることにより{θ1,pol(ω),…,θβ,pol(ω)}={0,10,20,…,180}となる（β=19）。テンプレート情報に含まれるテンプレートのうち集合{θ1,pol(ω),…,θβ,pol(ω)}の要素のいずれかを極角θpol(ω)として持つ方向情報（位置）に対応するテンプレートのインデックスの集合をΨとする。この際、集合Ψの要素の個数（濃度）|Ψ|が、Ｑ≦|Ψ|＜Ｐを満たすようになることが好ましい。集合Ψは、テンプレート情報に含まれるテンプレートに対応する方向情報（位置）の集合の真部分集合である。
【００９２】
そして、ｐ番目のテンプレートSp→(ω)が残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)を、最小二乗法に基づき、式（３８ａ）により求める。ただし、式（３８ａ）の右辺で用いるテンプレートのインデックスｐはｐ∈Ψである。なお、この段階では、式（３８ａ）左辺のｑは意味を持たないことに留意されたい。
【数３５】
【００９３】
次に、ｐ∈Ψなる条件の下で式（３８ａ）に基づき得られた|Ψ|個の係数Aq(ω,k,p)（１≦ｐ≦|Ψ|）の中から、その絶対値の大きい方から順にＱ個の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（３８ｂ）参照）。記号Γは、集合Ψから式（３８ｂ）により決定されたインデックスの集合を除いた集合であり、Γ＝Ψ-{g(ω,1),…,g(ω,q-1)}である。
【数３６】
【００９４】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（３８）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００９５】
これらの実施形態のほか、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）をランダムに設定する実施形態も許容される。
【００９６】
＜変形例＞
上述の実施形態では周波数ごとに観測信号X→(ω,k)を用いて反射音情報rs→(ω,k)を推定したが、周波数ごとに反射音情報を推定すると、一意に推定されるべき仮想音源の方向（推定到来方向）以外の方向に関する情報も含んでしまうことがあり、この結果、反射音情報に誤差が生じることがありうる。例えば、図１４（ａ）に示すように推定到来方向に関する情報だけを抽出できることが望ましいが、実際には図１４（ｂ）に示すように推定到来方向以外の方向に関する情報が混在してしまうことがありえる。
【００９７】
そこで変形例では、全周波数に亘り一括してパワーを算出することで、反射音情報の推定誤差を小さくする。つまり、図１５に示すように、残差信号のパワーを全周波数に亘り統括することによって、推定到来方向以外の方向の影響を極力減らすことができる。一般的に推定到来方向以外の方向では各周波数でのパワーにバラつきが生じるので、残差信号のパワーを全周波数に亘り統括することにより、推定到来方向のパワーに比してそれ以外の方向のパワーの相対的な影響を低減することができる。なお、図１５では、縦軸のパワーは相対値を示しているので各グラフのスケールが同じであるわけではないことに留意されたい。
【００９８】
この変形例での処理は次のとおりである。解析する周波数帯域に含まれる周波数のインデックスωの集合をΩとする。例えば、音声信号を扱うのであれば、1.0〜3.0kHz帯域に対応するインデックスの集合をΩとすればよい。そして、テンプレートSg(q,k)→(ω)のインデックスg(q,k)を式（１０）の替わりに式（３９）によって求める。
【数３７】
【００９９】
テンプレートSg(ω,q)→(ω)のインデックスg(ω,q)を式（１８）の替わりに式（４０）によって求める。また、方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（１９）の替わりに式（４１）による更新によって行われる。
【数３８】
【０１００】
テンプレートSg(ω,q)→(ω)のインデックスg(ω,q)を式（２８）や式（３７）の替わりに式（４２）によって求める。同様に、式（３８）の替わりに式（４３）によって求める。また、方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（３１）の替わりに式（４４）による更新によって行われる。
【数３９】
【実施例】
【０１０１】
本発明を用いて相対的な到来時刻差を推定した実験結果を示す。直方体の部屋で１００本のマイクロホンを１０行１０列に等間隔に２次元的に並べ、壁際に配置した環境をシミュレートした。マイクロホン間隔ｄは４ｃｍである。詳細な実験条件は図１１に示してある。テンプレートの座標として、図１０に示したように、マイクロホンアレーを取り囲む半球上に等間隔となるように配置した。1.0-3.0kHz帯域を解析した場合の実験結果を図１２に示す。基準反射音として直接音を選択した。図が煩雑になることを避けるため、ここでは反射音１−３についてのみ図示している。反射音１について、直接音からの到来時刻差は周波数平均で1.4msecであり、反射音２について、直接音からの到来時刻差は周波数平均で2.5msecであり、反射音３について、直接音からの到来時刻差は周波数平均で4.7msecであった。低周波数帯域では距離に対する位相の変化が小さいために到来時刻差の推定値の誤差が大きくなる傾向があるため、周波数帯域全体に亘って平均をとることによって到来時刻差を求めることがよい。
【０１０２】
＜応用例＞
反射音情報は人間が生活する上で、非常に重要な音声情報である。例えば、視覚障害者は、タッピングによって発した音源信号が壁や天井等で反射して耳で観測することにより、環境を把握している。また、日常会話でも、適度な反射が生じる部屋で会話することと、反射音が比較的少ない環境で会話することでは会話のしやすさに相違が生じる。以下、本発明により推定された反射音情報を用いたサービス例について述べる。
１つ目は、会議システムに本発明を組み込んだ例である。指向性音源の向きに応じて反射音の振幅は変化するので、反射音情報が分かると、どの方向に音源が向いているのかを推定することができる。会議システムに音源向きの推定装置を組み込めば、誰に向かって発言したのかを提示することに応用できる。
２つ目は、自由な位置で映像や音声を鑑賞できるシステムである。遠方にある音は直接到来する音源のパワーが小さいので収音することが困難である。反射音情報が分かると、直接音だけでなく、反射音も強調収音できるので、遠方の音を強調することが可能となる。また、音声処理の分野では、方向別に音源の強調収音は可能であるが、距離別に音声を強調収音することは非常に難しいとされている。反射音情報が分かると、距離に対応する物理的な特徴量が得られるので、距離別に収音することが可能となる。遠方の音を収音したり、方向別、距離別に収音することができれば、視聴者の選択した位置に対応した音場を擬似的に生成することが可能となる。
【０１０３】
音声通信システムにおいて、反射音情報を推定することは、直接音だけでは得られなかった音場の情報を得ることにつながる。反射音情報が分かれば、これまでの音声強調技術ではできなかったような遠方音の収音や距離別の収音に結びついたり、従来の収音技術では推定できなかった音場の情報（例えば音源の向き）を推定できる。こういった音場の情報の推定は、これまでの技術では実現できなかった音声処理装置の開発に繋がる。反射音情報の推定に関する従来技術は、インパルス応答を求めるために特殊な信号を観測する必要があったが、本発明は音声信号のような一般的な観測信号で反射音情報を得られるという利点を持つ。
【０１０４】
＜反射音情報推定装置のハードウェア構成例＞
上述の実施形態に関わる反射音情報推定装置は、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ＣＰＵ（Central Processing Unit）〔キャッシュメモリなどを備えていてもよい。〕、メモリであるＲＡＭ（Random Access Memory）やＲＯＭ（Read Only Memory）と、ハードディスクである外部記憶装置、並びにこれらの入力部、出力部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置間のデータのやり取りが可能なように接続するバスなどを備えている。また必要に応じて、反射音情報推定装置に、ＣＤ−ＲＯＭなどの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【０１０５】
反射音情報推定装置の外部記憶装置には、反射音情報を推定するためのプログラム並びにこのプログラムの処理において必要となるデータなどが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。以下、データやその格納領域のアドレスなどを記憶する記憶装置を単に「記憶部」と呼ぶことにする。
【０１０６】
反射音情報推定装置の記憶部には、アナログ信号に対してＡＤ変換を行うためのプログラム、フレーム分割処理を行うためのプログラム、フレームごとのデジタル信号を周波数領域の観測信号に変換するためのプログラム、テンプレート情報を生成するためのプログラム、周波数領域の観測信号とテンプレート情報を用いて反射音情報を推定するためのプログラムが記憶されている。
【０１０７】
反射音情報推定装置では、記憶部に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭに読み込まれて、ＣＰＵで解釈実行・処理される。この結果、ＣＰＵが所定の機能（ＡＤ変換部、フレーム分割部、周波数領域変換部、テンプレート生成部、反射音情報推定部、到来時刻差推定部）を実現することで反射音情報の推定が実現される。
【０１０８】
＜補記＞
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１０９】
また、上記実施形態において説明したハードウェアエンティティ（反射音情報推定装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１１０】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１１３】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【技術分野】
【０００１】
本発明は、反射音に関する情報、特に基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する技術に関する。
【背景技術】
【０００２】
電話や音声会議といった音声情報をやりとりするシステムを一般に音声通信システムと呼ぶ。音声通信システムにおいて、反射音に関する情報（到来振幅、到来方向など）を得ることは非常に重要なことである。会議室のような残響環境下において、マイクロホンを通して収音される収音信号の中には発話者のような音源から直接到来する直接音だけではなく、床、壁や天井に反射して到来する反射音が混在する。したがって、このような残響環境下で或る話者の発言を収録すると、直接音から遅延して反射音が混入するため、聞き取りづらくなってしまう。収音信号から各反射音の到来情報を推定して、反射音を除去することができれば、聞き取りやすい音声に回復することができる。ここで、反射音情報を推定する従来研究として、非特許文献１が挙げられる。
【０００３】
非特許文献１に開示される技術を実現する機能構成を図１に示す。この技術における処理手順は次のとおりである。
【０００４】
１．インパルス音源１００から放射された音源信号を４ｃｈのマイクロホン１１０−１，１１０−２，１１０−３，１１０−４を用いて収音する。ＡＤ変換部１２０は、収音されたアナログ信号をデジタル信号x→(t)＝[x1(t),x2(t),x3(t),x4(t)]Tへ変換する。ここで、[・]Tは転置を表す。ｔは離散時間のインデックスを表す。４本のマイクロホンは正四面体の頂点に配置されていることとする。
【０００５】
２．インパルス応答算出部１３０は、デジタル信号x→(t)＝[x1(t),x2(t),x3(t),x4(t)]Tを入力とし、各マイクロホンのインパルス応答h→(t)＝[h1(t),h2(t),h3(t),h4(t)]Tを算出する。インパルス応答の算出方法には、ＴＳＰ法やＭ系列法等があり、いかなる方法を用いてインパルス応答を算出してもよい。
【０００６】
３．仮想音源算出部１４０は、４ｃｈのインパルス応答h→(t)＝[h1(t),h2(t),h3(t),h4(t)]Tを入力とし、仮想音源情報v→＝[v→1,…,v→D]Tを出力する。Ｄは仮想音源の数を表す。仮想音源とは、各反射音の到来振幅、到来方向、到来時間を表現するために仮想的に存在するとされる音源である。図２を参照して、仮想音源について説明する。図２には、右側の壁で反射した音源信号をマイクロホンで受音する経路が書かれている。右側の壁で反射して到来する音源信号（反射音）は、「仮想音源」と書かれた位置から直接到来する信号と等価である（ただし、壁での反射による減衰や距離減衰の影響は受ける）。
【０００７】
この従来技術の詳細について説明する。インパルス応答を近接した４つの受音点（マイクロホンの位置）で測定すると反射音の到来時刻にわずかな差が生じる。インパルス応答の短い区間の相互相関を利用して、各マイクロホンにおける反射音の対応付けを行うことで、図３のように、ｎ番目の反射波に関するそれぞれの受音点での到来時刻ｔ1n，ｔ2n，ｔ3n，ｔ4n（１≦ｎ≦Ｄ）が求まる。正四面体マイクロホンアレーの辺の長さをｄ、音速をｃとすると、各仮想音源情報vn→＝[Xn,Yn,Zn,Sn]Tが求まる。ここで、Xn,Yn,Znはｎ番目の仮想音源の位置を表し（式（１）−（３）参照）、これは各反射音の到来方向と到来時間に対応する情報を持つ。また、Snはｎ番目の仮想音源の強さを表し、４ｃｈのインパルスで対応付けされたｎ番目の反射音の振幅の平均で求まる。
【数１】
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】山崎芳男ら、「近接する４点のインパルス応答により求めたホールの空間情報」、日本音響学会講演論文集、１９８１年５年、pp.759-760.
【発明の概要】
【発明が解決しようとする課題】
【０００９】
従来技術によると、仮想音源情報と呼んでいた反射音の「到来振幅」、「到来方向」、「到来時刻」を推定するためには、インパルス応答をあらかじめ用意することが必要であった。しかし、インパルス応答を用意するためには特殊な信号を用いて観測する必要があるため、あらゆる位置でのインパルス応答が事前に用意されているという条件は現実的ではない。
【００１０】
そこで本発明は、特殊な信号を用いることなく、反射音情報、特に基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する技術を提供することを目的とする。
【課題を解決するための手段】
【００１１】
反射音を、所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（伝達特性関数）に複素振幅を乗じて表されるとする。Ｑを２以上の予め定められた整数として、Ｑ個の反射音それぞれの複素振幅（到来振幅）を用いて、基準反射音以外の反射音（対象反射音）それぞれについて、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とする。
【００１２】
基準反射音として、Ｑ個の反射音の中から各対象反射音に共通の基準反射音を決定してもよいし、対象反射音ごとにＱ個の反射音の中から基準反射音を決定してもよい。
【００１３】
Ｑ個の反射音の到来振幅はそれぞれ、音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（観測信号）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものであるとしてもよい。
【００１４】
あるいは、Ｑ個の反射音の到来振幅は、相異なるＱ個の位置に対応するＱ個の反射音を観測信号から減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅であるとしてもよい。
【発明の効果】
【００１５】
本発明に拠ると、基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、基準反射音に対する対象反射音の相対的な到来時刻差とすることから、インパルス応答を求めるために音源信号に特殊な信号を用いることなく、反射音情報（基準反射音に対する他の反射音の相対的な到来時刻差）を推定することが可能である。反射音情報が得られると、従来の音声情報処理技術では実現できなかった音源距離の推定や、音声強調（遠方音の収音や距離別の収音）といった用途に応用できる。
【図面の簡単な説明】
【００１６】
【図１】従来技術における反射音情報推定技術の機能構成を示す図。
【図２】仮想音源を説明するための図。
【図３】従来技術における反射音の対応付けを説明するための図。
【図４】第１実施形態に係る反射音情報推定装置の機能構成を示す図。
【図５】第１実施形態に係る反射音情報推定方法の処理手順を示す図。
【図６】２次元マイクロホンアレーの構成例を示す図。
【図７】ｐ番目の点[xp,yp,zp]とｍ番目の受音点[um,vm,wm]との間の伝達特性を説明するための図。
【図８】図６に示すマイクロホンアレーを用いて観測した或る平面での音圧分布が例えば直接音と反射音１と反射音２との重畳で得られていることを説明するための図。
【図９】本発明の原理を説明するための図。
【図１０】実験例においてテンプレート情報作成のための区間座標とマイクロホンアレーとの位置関係を示す図。
【図１１】実験結果を示す図。
【図１２】実験例により本発明の有効性を説明するための図。
【図１３】第２実施形態に係る反射音情報推定方法の処理手順を示す図。
【図１４】（ａ）理想的には、推定到来方向に関する情報だけが抽出されるべきことを説明するための図。（ｂ）実際には、推定到来方向以外の方向に関する情報も混在してしまうことを説明するための図。
【図１５】残差信号のパワーを全周波数に亘って総括することにより、推定到来方向以外の方向の影響を減らすことを説明するための図。
【図１６】実用レベルの２次元マトリクスマイクロホンアレーを用いた場合における音圧分布とその分解を示す図。
【発明を実施するための形態】
【００１７】
《第１実施形態》
本発明は、発話信号のような音源から放射された音声信号（音源信号）を複数のマイクロホンで構成されるマイクロホンアレーで収音した信号（収音信号）から反射音の「到来振幅」と「到来方向」の少なくともいずれか一つを推定する。第１実施形態の機能構成および処理フローを図４と図５に示す。
【００１８】
音源２００から放射された音源信号をＭｃｈのマイクロホン２１０−１，…，２１０−Ｍを用いて収音する（ステップＳ１）。Ｍは、４より大きい値が望ましい。ＡＤ変換部２２０が、収音されたアナログ信号をデジタル信号xx→(t)＝[xx1(t),…,xxM(t)]Tへ変換する（ステップＳ２）。ここで、[・]Tは転置を表す。ｔは離散時間のインデックスを表す。
【００１９】
Ｍ本のマイクロホンの並べ方は、２次元または３次元的に等間隔で配置することが望ましい。これは、反射音の到来方向とテンプレート（後で説明するが、反射音の伝達特性を模擬したものである）の対応を一意に定めるためである。なお、原理的に、１次元的にマイクロホンを配置しても、あるいは等間隔に配置しなくても、本発明を実施できるが、反射音の伝達特性と反射音の到来方向が一対一の関係とならないため、２次元または３次元的に等間隔で配置することが望ましい。２次元平面上に等間隔にマイクロホンを並べた場合の一例を図６に示す。マイクロホン間隔ｄは、空間サンプリング定理を満たすように設定されていることが望ましい。空間サンプリング定理を満たす場合、マイクロホン間隔ｄは、式（４）を満たす数値となる。ｃは音速であり、ｆは解析対象とする周波数である。例えば、４ｋＨｚの周波数を解析する場合、マイクロホン間隔を４ｃｍ程度に設定するのがよい。
【数２】
【００２０】
フレーム分割部２３０は、ＡＤ変換部２００が出力したデジタル信号xx→(t)＝[xx1(t),…,xxM(t)]Tを入力とし、チャネルごとに複数サンプルから成るデジタル信号の組（フレーム）に分割された信号x→(k)＝[x1(k),…,xM(k)]Tを出力する（ステップＳ３）。ｋはフレーム番号を表すインデックスである。フレーム分割は、各チャネルのデジタル信号xxi(t)（１≦ｉ≦Ｍ）ごとにＷ点分をバッファリングして出力する処理である。Ｗはサンプリング周波数にもよるが、１６ｋＨｚサンプリングの場合には５１２点あたりが妥当である。
【００２１】
周波数領域変換部２４０は、各フレームのデジタル信号x→(k)を入力として、周波数領域の信号X→(ω,k)＝[X1(ω,k),…,XM(ω,k)]Tに変換して出力する（ステップＳ４）。この信号X→(ω,k)を観測信号と呼ぶことにする。ここで、ωは離散周波数のインデックスを指し（周波数ｆと角周波数ωとの間にはω=2πfの関係があるから、周波数のインデックスωをこの角周波数ωと同一視してもかまわない。以下、ωに関して「周波数のインデックス」を単に「周波数」ともいう）、ｋはフレームのインデックスを指す。周波数領域に変換する方法の一つに、離散フーリエ変換があるが、周波数領域に変換するのであれば、他の方法を用いてもよい。周波数領域の観測信号X→(ω,k)は、各周波数ω、フレームｋごとに出力される。
【００２２】
テンプレート生成部２５０は、Ｐ個のテンプレートSp→(ω)の集合（ただし、計算の都合、ベクトル表記している）であるテンプレート情報S→(ω)＝[S1→(ω),…,SP→(ω)]（∀ω∈Ω；Ωは周波数のインデックスωの集合）を周波数ωごとに生成する（ステップＳｐ）。この処理は通常、ステップＳ１−Ｓ４の各処理に先立ち実施される。Ｐはテンプレートの総数を表し、予め２以上の整数値に設定されている。テンプレート総数Ｐは多ければ多いほど高精度な反射音情報の推定に繋がるが、計算量が多くなるので、例えばＰ＝１０００くらいに設定するのが良い。この処理は、マイクロホンで信号を観測する以前にあらかじめ行う処理である。また、マイクロホンの位置（例えばマイクロホンの間隔ｄ）を変更したり、テンプレート総数Ｐを変更したりしない限り、通常、テンプレートを毎回作り直す必要はない。ここで言う“テンプレート”とは、反射音の到来方向に対応する伝達特性（音響伝播特性）を模擬したものである。ｐ番目（１≦ｐ≦Ｐ）のテンプレートSp→(ω)＝[Sp1(ω),…,SpM(ω)]T（ω∈Ω）は、予め定められたｐ番目の点[xp,yp,zp]とＭ個の受音点（ここで受音点はマイクロホンが配置される位置であり、ｍ番目（１≦ｍ≦Ｍ）の受音点を[um,vm,wm]とする）の間の周波数ごとの伝達特性を表す（図７を参照）。ｐ番目のテンプレートSp→(ω)の各要素Spm(ω)の算出式の一例を式（５）に示す。記号ｉは虚数単位を表す。
【数３】
【００２３】
ｐ番目のテンプレートSp→(ω)には方向情報θp→(ω)が対応付けられている。方向情報θp→(ω)は、ｐ番目の点[xp,yp,zp]および受音点[um,vm,wm]の位置座標の基準となる３次元直交座標系の原点からｐ番目の点[xp,yp,zp]を見た方向であり、例えば（当該３次元直交座標系の原点と共通の原点を持つ）球座標系における二つの偏角（極角θp,polと方位角θp,azi）として表される。つまり、θp→(ω)＝[θp,pol(ω),θp,azi(ω)]である。なお、ｐ番目のテンプレートSp→(ω)にｐ番目の点[xp,yp,zp]が関連付けられていれば方向情報θp→(ω)は位置[xp,yp,zp]から計算可能であるから、ｐ番目のテンプレートSp→(ω)に方向情報θp→(ω)が対応付けられていることは必須要件ではない。なお、３次元直交座標系と球座標系とは相互に転換可能であるから（座標変換）、式（５）の右辺は位置[x,y,z]でなく方向情報θp→(ω)＝[θp,pol(ω),θp,azi(ω)]を用いて、例えば式（５ａ）のように表すこともできる。ここで、ｄはマイクロホン間隔であり、マイクロホンアレーをΦ行Ξ列（Φ×Ξ＝Ｍ）の２次元マイクロホンアレーとし、ｍ番目のマイクロホンの位置をφ行ξ列（１≦φ≦Φ，１≦ξ≦Ξ）にあるとする。
【数４】
【００２４】
また、第１実施形態のようにテンプレートが方向に対応している場合、Ｐ個の点[xp,yp,zp]（１≦ｐ≦Ｐ）の位置は互いに方向の異なる位置であることが好ましく、例えば各点[xp,yp,zp]が原点から十分に離れた等距離にあるとして、上記原点を中心とする球面上の異なるＰ個の点とすればよい。各点[xp,yp,zp]を原点から十分に離れた位置とする理由は、音源ないし仮想音源から放射された信号は球面的に伝達するが音源ないし仮想音源から十分に離れた位置（原点）での局所領域では直接音ないし反射音を平面波として模擬できるからである。ただし、テンプレート情報が同じ方向の位置に対応するテンプレートを含むことを排除する趣旨ではない。なお、マイクロホンアレーは上記座標系の原点の近傍（局所領域）に配置されているとする。
【００２５】
テンプレート記憶部２６０は、テンプレート生成部２５０が出力したテンプレート情報S→(ω)を記憶し、解析時に反射音情報推定部２７０にテンプレート情報S→(ω)を提供する役割を果たす。
【００２６】
反射音情報推定部２７０は、周波数領域の観測信号X→(ω,k)とテンプレート情報S→(ω)を入力として、Ｑ個の反射音情報成分rsq→(ω,k)の集合（ただし、計算の都合、ベクトル表記している）である反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを各フレームｋについて周波数ωごとに出力する（ステップＳ５）。ここで、Ｑは推定される反射音の総数を表し、予め１以上の整数値に設定されている。ｑ番目（１≦ｑ≦Ｑ）の反射音情報成分rsq→(ω,k)は、rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]の２要素から成り、rsAq(ω,k)はｑ番目の反射音の到来振幅であり、rsBq(ω,k)はｑ番目の反射音の到来方向である。
【００２７】
反射音情報を推定する原理について説明する。図６に示すような２次元マイクロホンアレーを用いて観測した或る平面での音圧分布の一例を図８の左端の濃淡図として示す。濃淡図として示された音圧分布の見方について、黒い部分は音圧が小さく、白い部分は音圧が大きいことを示す。観測した音圧分布には直接音の音圧分布だけではなく、反射音の音圧分布も混入している。直接音や反射音が十分に遠方より到来する場合において、２次元平面上でのそれぞれの音圧分布は、図８の右側の３つの濃淡図のように縞模様となる。縞模様の「濃淡」が直接音ないし反射音の到来振幅、「回転・周期」が直接音ないし反射音の到来方向にそれぞれ対応する。図８の例では、到来振幅や到来方向が異なる直接音、反射音１、反射音２の各音圧分布の重畳で観測信号の音圧分布が構成されることを示している。周波数領域で考えると、直接音や各反射音は到来方向に応じて周波数の変化する複素正弦波で表され、観測信号は直接音と各反射音に対応する複数の複素正弦波が重畳したものとして表される。ところで、本発明で解決する問題は、観測信号のみを用いて、反射音の到来振幅および／または到来方向を推定することである。この課題解決は、図８の左端に描かれた音圧分布から図８の右側の３つの濃淡図の直接音や各反射音に対応する縞模様の「濃淡」および／または「回転・周期」を推定することに対応する。
【００２８】
図９を参照して、反射音情報rs→(ω,k)を推定する手法の概略について説明する。ある２次元平面で観測した観測信号に含まれているパワーの最も強い反射音０（q=1に相当し、最も強いパワーを持つことから、通常、この反射音０は「直接音」として理解される）を推定し、観測信号から当該反射音０を減算して残差信号E2を得る。次に当該残差信号E2に含まれているパワーの最も強い反射音１（q=2に相当する）を推定し、当該残差信号E2から当該反射音１を減算して新たな残差信号E3を得る。次に、当該残差信号E3に含まれているパワーの最も強い反射音２（q=Q=3に相当する）を推定する。ここでは、Ｑ＝３の場合を説明したが、一般的に、ｑ番目の残差信号Eq（ただし１番目の残差信号は観測信号とする）に含まれているパワーの最も強いｑ番目の反射音q-1（ただし反射音０は直接音である）を減算する操作をｑ＝Ｑまで逐次実行することでＱ個の反射音情報成分（rs1→(ω,k),…,rsQ→(ω,k)）を得る。１番目の反射音情報成分rs1→(ω,k)は反射音０（直接音）に対応し、２番目の反射音情報成分rs2→(ω,k)は反射音１に対応し、３番目の反射音情報成分rs3→(ω,k)は反射音２に対応し、・・・、Ｑ番目の反射音情報成分rsQ→(ω,k)は反射音Q-1に対応する。Ｑは、計算パワーや反射音情報を用いるアプリケーションにも依存するが、３０くらいに設定するのが良い。
【００２９】
なお、図８および図９の音圧分布はそれぞれ高解像度の濃淡図として示されているが、このような高解像度の濃淡図として音圧分布を示すためには極めて多くのマイクロホンを必要とし、実用的ではない。他方、実用レベルの２次元マトリクスマイクロホンアレーとして例えば１００個のマイクロホンを１０×１０の２次元マトリクスマイクロホンアレーとして用いた場合でさえ、粗い（低解像度）濃淡図（図１６参照）として示される音圧分布しか得られない。そこで、実用の観点から、低解像度の音圧分布しか得られないような状況の下で、精度良く反射音の到来振幅や到来方向を推定することが求められる。本発明では、空間分解能の向上のために任意の位置から到来する平面波を具体的に表現することとし（定式化）、パワーが大きな反射音の影響を受けてパワーの小さいな反射音を推定できなくなることを防止するために、既に推定された反射音を観測信号から除去して次の反射音を推定する（分解）。定式化についてはテンプレート情報として説明したとおりであり、分解については反射音情報rs→(ω,k)の推定手法の概略で説明したとおりである。
【００３０】
上で述べた反射音情報rs→(ω,k)を推定する手法について詳細を述べる。説明に先立ち、記号の定義を行う。ｑ番目の残差信号をEq→(ω,k)=[Eq1(ω,k),…,EqM(ω,k)]T、推定されたｑ番目の反射音（q=1の場合は直接音を表す）をAq(ω,k,g(q,k))Sg(q,k)→(ω)とする。ここで、g(q,k)は、フレームｋにおけるｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(q,k))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(q,k)→(ω)と反射音との相違を表す。ｑの昇順に残差信号から反射音を減算する上述の方法を式で表すと式（６）のようになる。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数５】
【００３１】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（７）、式（８）で与えられる。方向情報θg(q,k)→(ω,k)はテンプレートSg(q,k)→(ω)に対応する方向情報である。
【数６】
【００３２】
次に、ｑ番目の反射音Aq(ω,k,g(q,k))Sg(q,k)→(ω)を推定する手法について述べる。ｑ番目の反射音Aq(ω,k,g(q,k))Sg(q,k)→(ω)は、ｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするものとして推定される。記号Ｈは共役転置を表す。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(q,k))とSg(q,k)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の＜処理１＞と＜処理２＞はｑの昇順に各ｑについて行われる。
【００３３】
＜処理１＞
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（１０）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(1,k),…,g(q-1,k)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目のテンプレートSp→(ω)が残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（９）により求められる。なお、この段階では、式（９）左辺のｑは意味を持たないことに留意されたい。
【数７】
【００３４】
＜処理２＞
集合Λの要素の個数（濃度）を|Λ|とすると、式（９）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）を用いて、テンプレートSg(q,k)→(ω)のインデックスを表すg(q,k)は、残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするインデックスとして式（１０）により得られる。
【数８】
【００３５】
ｑ番目の反射音を構成する係数Aq(ω,k,g(q,k))は、式（１０）で得られたg(q,k)を式（９）のｐに代入して得られる。式（１０）の計算の際に係数Aq(ω,k,g(q,k))は計算されているので、メモリに余裕がある場合には、式（１０）の際に得られた係数Aq(ω,k,g(q,k))をメモリに格納しておき、式（１０）でg(q,k)が得られた後に、メモリからp=g(q,k)に対応する係数Aq(ω,k,g(q,k))を呼び出せばよい。
【００３６】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。
【００３７】
到来時刻差推定部２８０は、反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを入力として（正確には、複素振幅{rsA1(ω,k),…,rsAQ(ω,k)}のみで十分である）、基準となる反射音（基準反射音）に対する各反射音の相対的な到来時刻差rsDq(k)（１≦ｑ≦Ｑ）の集合である到来時刻差情報rsD→(k)＝[rsD1(k),…,rsDQ(k)]T（ただし、到来時刻差を用いる任意のアプリケーションの便宜上、ベクトル表記している）を、各フレームｋについて出力する（ステップＳ６）。具体的には、後述の((処理１))と((処理２))を各ｑ（１≦ｑ≦Ｑ）について行うことにより、rsDq(k)（１≦ｑ≦Ｑ）が求められる。ここでは処理対象のｑを特にｑｃと記す。
【００３８】
((処理１))
推定されたＱ個の反射音から基準反射音のインデックスＮを決める。((処理１))を各ｑ（１≦ｑ≦Ｑ）について行う場合、Ｎ＝ｑｃならば基準反射音に対するｑｃ番目の反射音の到来時刻差が０であることは自明であるから、インデックスＮは反射音を識別するインデックスの差集合{１，…，Ｑ}−{ｑｃ}の中から選択されることが好ましい。基準反射音の決定方法は種々考えられるが、例えば、差集合{１，…，Ｑ}−{ｑｃ}に対応する反射音のうちパワーが強い反射音を選択することが好適である。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して到来振幅の大きさの総和を計算することにより、周波数による影響をキャンセルすることができる。この場合、Ｎは式（１０ａ）により決定される。なお、解析する周波数帯域に含まれる周波数のインデックスωの集合をΩとし、例えば音声信号を扱うのであれば、1.0〜3.0kHz帯域に対応するインデックスの集合をΩとすればよい。
【数９】
【００３９】
((処理２))
Ｎ番目の反射音の到来振幅rsAN(ω,k)とｑｃ番目（ｑｃ≠Ｎ）の反射音の到来振幅rsAqc(ω,k)との位相差を用いて相対的な到来時刻差rsDqc(k)を算出する。一般に、到来振幅rsAq(ω,k)は複素振幅であるから振幅λq(ω,k)と時刻τq(ω,k)を用いてrsAq(ω,k)=λq(ω,k)exp[iωτq(ω,k)]と表すことができる。従って、Ｎ番目の反射音に対するｑ番目（ｑ≠Ｎ）の反射音の位相差は、Ｎ番目の反射音の到来振幅rsAN(ω,k)に対するｑ番目（ｑ≠Ｎ）の反射音の到来振幅rsAq(ω,k)の偏角arg(rsAq(ω,k)／rsAN(ω,k))として与えられる。偏角arg(rsAq(ω,k)／rsAN(ω,k))を周波数のインデックスωで除することによりＮ番目の反射音とｑ番目（ｑ≠Ｎ）の反射音との時刻差λq(ω,k)-λN(ω,k)を求めることができる。つまり、λq(ω,k)-λN(ω,k)=arg(rsAq(ω,k)／rsAN(ω,k))／ωである。時刻λN(ω,k)に対する時刻λq(ω,k)が「早い」か「遅い」かの区別を行うため負号を導入する。つまり、-arg(rsAq(ω,k)／rsAN(ω,k))／ω＜０であれば時刻λN(ω,k)に対する時刻λq(ω,k)が「早い」ことを表し、-arg(rsAq(ω,k)／rsAN(ω,k))／ω＞０であれば時刻λN(ω,k)に対する時刻λq(ω,k)が「遅い」ことを表す。-arg(rsAq(ω,k)／rsAN(ω,k))／ω=０の場合、時刻差が無いことを表す。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して-arg(rsAq(ω,k)／rsAN(ω,k))／ωの総和を計算することにより、周波数による影響をキャンセルすることができる。即ち、相対的な到来時刻差rsDqc(k)は、式（１０ｂ）によって与えられる。集合Ωの要素の個数（濃度）を|Ω|とする。
【数１０】
【００４０】
なお、（(処理１)）にて基準反射音のインデックスを決定するが、基準反射音を各ｑについて共通とする場合には、（(処理１)）はひとたび行われれば十分である。この場合、先ず((処理１ａ))を行い、次いで差集合{１，…，Ｑ}−{Ｎ}の各要素ｑｃについて上述の((処理２))を行えばよい。
【００４１】
((処理１ａ))
推定されたＱ個の反射音から基準反射音のインデックスＮを決める。基準反射音の決定方法は種々考えられるが、例えば、Ｑ個の反射音のうちパワーが強い反射音（直接音）を選択することが好適である。この際、解析する周波数帯域に含まれる周波数のインデックスωに関して到来振幅の大きさの総和を計算することにより、周波数による影響をキャンセルすることができる。この場合、Ｎは式（１０ｃ）により決定される。
【数１１】
【００４２】
なお、マイクロホンアレーから音源ないし仮想音源までの音源距離disq(ω,k)は、相対的な到来時刻差を定数倍した量であるから、式（１０ｄ）によって与えられる。ｃは音速である。
【数１２】
【００４３】
以下、相対的な到来時刻差を求めるに必要な反射音情報rs→(ω,k)＝[rs1→(ω,k),…,rsQ→(ω,k)]Tを求めるいくつかの実施形態を説明する。いずれの実施形態であっても、第１実施形態のステップＳ６の処理を行うことで、相対的な到来時刻差rsDq(k)（１≦ｑ≦Ｑ）を求めることができる。
【００４４】
《第２実施形態》
第１実施形態ではテンプレート情報S→(ω)を用いて反射音情報rs→(ω,k)を求めたが、Ｐ個のテンプレートSp→(ω)の集合であるテンプレート情報S→(ω)を事前に求めておくことは必ずしも必須ではない。テンプレート情報S→(ω)を事前に求めておかない実施形態を第２実施形態として説明する。
【００４５】
第２実施形態では、第１実施形態におけるステップＳ１−Ｓ４の各処理が実施されるが、第１実施形態におけるステップＳｐの処理が不要であり、さらに第１実施形態のステップＳ５の処理に替えてステップＳ５ａの処理が行われる（図１３参照）。そこで、第１実施形態と同じ事項については重複説明を省略し、第１実施形態と異なる事項について説明する。
【００４６】
記号の定義を行う。ｑ番目の残差信号をEq→(ω,k)=[Eq1(ω,k),…,EqM(ω,k)]T、ｑ番目の反射音（q=1の場合は直接音を表す）をAq(ω,k)Rq→(ω,θq→(ω,k))とする。反射音を構成するRq→(ω,θq→(ω,k))=[R1(ω,θq→(ω,k)),…,RM(ω,θq→(ω,k))]Tは、空間中の任意の位置[x,y,z]と各マイクロホンとの間の周波数ごとの伝達特性を模擬した関数（以下、伝達特性関数という）であり、各マイクロホンに対する伝達特性を模擬した関数であれば何でもよい。通常、伝達特性関数を構成する各伝達特性Rm(ω,θq→(ω,k))とテンプレートの各要素Spm(ω)の算出式とは同じである。この場合、方向情報θq→(ω,k)で表される方向に在る位置[x,y,z]とｍ番目の受音点[um,vm,wm]との間の周波数ごとの伝達特性Rm(ω,θq→(ω,k))は式（１１）で表される。なお、方向情報θq→(ω,k)で表される方向に在る位置[x,y,z]は、例えば、上記座標系原点から十分に離れた球面上の位置とすればよい。位置[x,y,z]を原点から十分に離れた位置とする理由は既述のとおりであり、詳しくは位置[x,y,z]はマイクロホンアレーが配置されている局所領域にて音源ないし仮想音源からの直接音ないし反射音を平面波として模擬できる距離にある空間中の任意の位置であることが好ましい。なお、３次元直交座標系と球座標系とは相互に転換可能であるから（座標変換）、式（１１）の右辺は位置[x,y,z]でなく方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]を用いて、例えば式（１１ａ）のように表すこともできる。ここで、ｄはマイクロホン間隔であり、マイクロホンアレーをΦ行Ξ列（Φ×Ξ＝Ｍ）の２次元マイクロホンアレーとし、ｍ番目のマイクロホンの位置をφ行ξ列（１≦φ≦Φ，１≦ξ≦Ξ）にあるとする。
【数１３】
【００４７】
反射音を構成するAq(ω,k)は、音源２００自身が持つ位相や壁での反射、距離による減衰といったテンプレートRq→(ω,θq→(ω,k))と反射音との相違を表し、到来振幅に相当する。ｑの昇順に残差信号から反射音を減算する上述の方法を式で表すと式（１２）のようになる。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数１４】
【００４８】
次に、反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する方法について説明する。
ｑ番目の最適化された反射音Aq(ω,k)Rq→(ω,θq→(ω,k))は、式（１２）で表されるｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って決定される。具体的には、伝達特性関数Rq→(ω,θq→(ω,k))が方向情報θq→(ω,k)で決定されることに注意すると、ｑ番目の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を表現するパラメータAq(ω,k)，θq→(ω,k)の最適値Aq,opt(ω,k)，θq,opt→(ω,k)は式（１３）によって得られる。なお、記号Ｈは共役転置を表す。
【数１５】
【００４９】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（１４）、式（１５）で与えられる。
【数１６】
【００５０】
式（１３）の具体的な計算方法は種々考えられるが、ここではその一例を示す。方向情報θq→(ω,k)を、例えばビームフォーマ法などの到来方位推定方法によって定める。ビームフォーマ法は、指向性ビームを空間走査し、得られた電力スペクトルから電力が大きくなる方向を探索する方法である。ここでは、ビームフォーマ法によりＰ個の到来方向が推定できたとする。
【００５１】
実際には、ビームフォーマ法によって得られる電力スペクトルは到来方向に対して急峻でないことがあり、このような場合、例えば、予め定めたスペクトル強度以上のスペクトル強度を示す電力スペクトルに対応する方向の範囲にて予め定めた間隔で到来方向を定めればよい。具体例として、極角５°、方位角１０°から２０°の範囲で予め定めたスペクトル強度以上のスペクトル強度を示す電力スペクトルが得られたとすると、予め定めた間隔２°ごとに到来方向を定めるとして、（極角５°，方位角１０°），（極角５°，方位角１２°），（極角５°，方位角１４°），（極角５°，方位角１６°），（極角５°，方位角１８°），（極角５°，方位角２０°）を到来方向とすればよい。
【００５２】
また、電力スペクトルが或る方向にて急峻なピークを示したとしても、単純に当該方向を到来方向の一つとして定めるのではなく、当該方向の所定の範囲で到来方向を定めてもよい。具体例として、極角３０°、方位角５０°で急峻なピークを示す電力スペクトルが得られたとすると、所定の範囲（極角±４°，方位角±４°，間隔２°）で到来方向を定めるとして、（極角２６°，方位角４６°），（極角２８°，方位角４６°），（極角３０°，方位角４６°），（極角３２°，方位角４６°），（極角３４°，方位角４６°），（極角２６°，方位角４８°），（極角２８°，方位角４８°），（極角３０°，方位角４８°），（極角３２°，方位角４８°），（極角３４°，方位角４８°），（極角２６°，方位角５０°），（極角２８°，方位角５０°），（極角３０°，方位角５０°），（極角３２°，方位角５０°），（極角３４°，方位角５０°），（極角２６°，方位角５２°），（極角２８°，方位角５２°），（極角３０°，方位角５２°），（極角３２°，方位角５２°），（極角３４°，方位角５２°），（極角２６°，方位角５４°），（極角２８°，方位角５４°），（極角３０°，方位角５４°），（極角３２°，方位角５４°），（極角３４°，方位角５４°）を到来方向とすればよい。なお、第１実施形態ではＰは固定値であったが、第２実施形態ではＰはビームフォーマ法などの到来方位推定方法による推定成果に依存する値であることに留意されたい。
【００５３】
ビームフォーマ法によって得られたＰ個の到来方向についてテンプレートを生成する。テンプレートの各要素の算出式は例えば式（６）である。このＰ個のテンプレート（テンプレート情報S→(ω)）を用いて第１実施形態で説明した処理を行えばよい。
【００５４】
《第３実施形態》
第３実施形態は、第２実施形態と異なり、テンプレート情報を利活用して反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する。式（１３）の具体的な計算方法の一例を示す。下記に説明する最適化方法はｑの昇順に各ｑに対して適用される。
【００５５】
§１ 方向情報の初期値設定
最初に、方向情報θq→(ω,k)の初期値θini,q→(ω,k)をテンプレート情報S→(ω)を用いて決定する。このために、推定されるべき到来方向に最も近いと考えられる方向に対応するテンプレートを決定し、この決定されたテンプレートに対応する方向情報を方向情報θq→(ω,k)の初期値θini,q→(ω,k)とすればよい。
【００５６】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、ｑ＋１番目の残差信号Eq+1→(ω,k)は式（１６）のように表される。ただしE1→(ω,k)=X→(ω,k)である。
【数１７】
【００５７】
反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（１６）に基づくｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の=処理１=と=処理２=はｑの昇順に各ｑについて行われる。
【００５８】
=処理１=
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（１８）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目のテンプレートSp→(ω)が残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（１７）により求められる。なお、この段階では、式（１７）左辺のｑは意味を持たないことに留意されたい。
【数１８】
【００５９】
=処理２=
集合Λの要素の個数（濃度）を|Λ|とすると、式（１７）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）を用いて、テンプレートSg(ω,q)→(ω)のインデックスを表すg(ω,q)は、残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするインデックスとして式（１８）により得られる。
【数１９】
【００６０】
従って、方向情報θq→(ω,k)の初期値θini,q→(ω,k)は、式（１８）により得られるg(ω,q)をインデックスに持つテンプレートSg(ω,q)→(ω)に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００６１】
§２ 反射音の最適化
次に、方向情報θq→(ω,k)の初期値θini,q→(ω,k)を起点として、式（１２）で表されるｑ＋１番目の残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とするように、反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を最適化する。反射音は、係数Aq(ω,k)とRq→(ω,θq→(ω,k))の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。この最適化方法は様々あるが、そのうちの一つの方法（勾配法）について述べる。例示する方法では、方向情報θq→(ω,k)の補正と係数Aq(ω,k)の補正が交互に所定回数（δ回）反復して行われることにより反射音Aq(ω,k)Rq→(ω,θq→(ω,k))が最適化される。δは例えば５０程度の値とされるが１でもよい。
【００６２】
§２．１ 方向情報の補正
方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（１９）による更新によって行われる。初めて§２．１の処理を行う場合、式（１９）右辺の方向情報θq→(ω,k)は§１の処理で得られた初期値θini,q→(ω,k)であり、§２．１の処理が初めてではない場合、式（１９）右辺の方向情報θq→(ω,k)は直前の§２．１の処理で得られた方向情報とする。また、初めて§２．１の処理を行う場合、パワー(Eq+1→(ω,k))HEq+1→(ω,k)の計算に用いられる係数Aq(ω,k)は式（１７）で得られたAq(ω,k,p)とし、§２．１の処理が初めてではない場合、パワー(Eq+1→(ω,k))HEq+1→(ω,k)の計算に用いられる係数Aq(ω,k)は直前の§２．２の処理（後述する）で得られた係数Aq(ω,k)とする。ステップ幅α1およびα2は小さい正の定数であり、収束速度などを考慮して決定されるが、例えばそれぞれ0.1程度の値とされる。
【数２０】
【００６３】
§２．２ 係数の補正
係数Aq(ω,k)の補正は、最小二乗法に基づき、式（２０）に従って新たな係数Aq(ω,k)を求めることにより行われる。式（２０）で用いるRq→(ω,θq→(ω,k))は§２．１の処理で得られた方向情報θq→(ω,k)と式（６）から得られる。
【数２１】
【００６４】
δ回の反復処理が終了した時点で得られている係数Aq(ω,k)と方向情報θq→(ω,k)がAq,opt(ω,k)とθq,opt→(ω,k)であり、ｑ番目の反射音情報成分rsq→(ω,k)となる。すなわち、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（２１）、式（２２）で与えられる。
【数２２】
【００６５】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。なお、δ＝１に設定されている場合、係数の補正を行わないことにより、反射音情報として到来方向のみを求めることができる。
【００６６】
《第４実施形態》
次に、第１、第２、第３実施形態と異なり、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))を一括して最適化する方法について説明する。
【００６７】
Ｑ個の最適化された反射音Aq(ω,k)Rq→(ω,θq→(ω,k))は、観測信号からＱ個の反射音を除去して得られる残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)（式（２３ａ）参照）を最小とする基準に従って決定される。具体的には、伝達特性関数Rq→(ω,θq→(ω,k))が方向情報θq→(ω,k)で決定されることに注意すると、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）を表現する各パラメータAq(ω,k)，θq→(ω,k)（１≦ｑ≦Ｑ）の最適値Aq,opt(ω,k)，θq,opt→(ω,k)（１≦ｑ≦Ｑ）は式（２３ｂ）によって得られる。なお、記号Ｈは共役転置を表す。式（２３ｂ）にて、｛(Aq,opt(ω,k)，θq,opt→(ω,k))｝q∈{1,…,Q}は、｛(A1,opt(ω,k)，θ1,opt→(ω,k)),…,(Aq,opt(ω,k)，θq,opt→(ω,k)),…,(AQ,opt(ω,k)，θQ,opt→(ω,k))｝を表し、｛(Aq(ω,k)，θq→(ω,k))｝q∈{1,…,Q}は、｛(A1(ω,k)，θ1→(ω,k)),…,(Aq(ω,k)，θq→(ω,k)),…,(AQ(ω,k)，θQ→(ω,k))｝を表す。
【数２３】
【００６８】
このとき、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（２４）、式（２５）で与えられる。
【数２４】
【００６９】
式（２３ｂ）の具体的な計算方法は種々考えられるが、ここではその一例を示す。
【００７０】
§１．１ 方向情報の初期値設定
最初に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を設定する。Ｑ個の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定方法として、この実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を、観測信号X→(ω,k)とテンプレート情報S→(ω)を用いて決定するビーム探索的決定方法を説明する。この方法によると、Ｑ個の推定されるべき到来方向それぞれに最も近いと考えられるＱ個の方向に対応するＱ個のテンプレートを決定し、この決定されたＱ個のテンプレートに対応する方向情報をＱ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）とすればよい。なお、この場合、ＰとＱとの間にＱ＜Ｐなる関係がある。
【００７１】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、ｑ番目の反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、観測信号X→(ω,k)からｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)を除去して得られる残差信号Eq→(ω,k)は式（２６）のように表される。
【数２５】
【００７２】
ｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（２６）に基づく残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。
【００７３】
反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。まず、ｐ番目のテンプレートSp→(ω)が残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（２７）により求められる。なお、この段階では、式（２７）左辺のｑは意味を持たないことに留意されたい。
【数２６】
【００７４】
次に、式（２７）に基づき得られたＰ個の係数Aq(ω,k,p)（１≦ｐ≦Ｐ）の中から、その絶対値の大きい方から順にＱ個の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（２８）参照）。記号Λは、インデックスｐの全体の集合{1,…,p,…,P}から式（２８）により決定されたインデックスの集合を除いた集合であり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}である。
【数２７】
【００７５】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（２８）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００７６】
§１．２ 係数Aq(ω,k)の初期値設定
次に、Ｑ個の係数Aq(ω,k)の初期値Aini,q(ω,k)を設定する。Ｑ個の初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）の決定方法として種々のものが考えられるが、ここでは一例として、パワー最小化基準でＱ個の初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）を決定する方法を説明する。まず、Aq(ω,k,p)=0（１≦ｐ≦Ｐ）とする。そして、初期値Aini,q(ω,k)（１≦ｑ≦Ｑ）は、残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)を最小化するように、最小二乗法に基づき、式（２９）により求められる。式（２９）にてFq→(ω,k)は式（３０）で与えられる。式（３０）にてΥ={1,…,q-1,q+1,…,Q}であり、Fq→(ω,k)は観測信号からｑ番目の反射音を除去した残差信号である。なお、Ｑ個の方向情報θq→(ω,k)として§１．１で決定されたＱ個の初期値θini,q→(ω,k)を用いる。式（２９）で用いるRq→(ω,θq→(ω,k))は方向情報の初期値θini,q→(ω,k)と式（６）から得られる。
【数２８】
【００７７】
§２ 反射音の最適化
次に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を起点として、式（２３ａ）で表される残差信号E→(ω,k)のパワー(E→(ω,k))HE→(ω,k)を最小とするように、Ｑ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）を一括して最適化する。各反射音は、係数Aq(ω,k)とRq→(ω,θq→(ω,k))の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。この最適化方法は様々あるが、そのうちの一つの方法（勾配法）について述べる。例示する方法では、方向情報θq→(ω,k)の補正と係数Aq(ω,k)の補正が交互に所定回数（δ回）反復して行われることによりＱ個の反射音Aq(ω,k)Rq→(ω,θq→(ω,k))（１≦ｑ≦Ｑ）が最適化される。δは例えば１００程度の値とされるが１でもよい。
【００７８】
§２．１ 方向情報の補正
Ｑ個の方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]（１≦ｑ≦Ｑ）の補正は、式（３１）による更新によって行われる。各ｑ（１≦ｑ≦Ｑ）について、初めて§２．１の処理を行う場合、式（３１）右辺の方向情報θq→(ω,k)は§１．１の処理で得られた初期値θini,q→(ω,k)であり、§２．１の処理が初めてではない場合、式（３１）右辺の方向情報θq→(ω,k)は直前の§２．１の処理で得られた方向情報とする。また、初めて§２．１の処理を行う場合、パワー(Fq→(ω,k))HFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として式（２９）で得られたAini,q(ω,k)（１≦ｑ≦Ｑ）を用い、§２．１の処理が初めてではない場合、パワー(Fq→(ω,k))HFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として直前の§２．２の処理（後述する）で得られた係数Aq(ω,k)（１≦ｑ≦Ｑ）を用いる。ステップ幅α1およびα2は小さい正の定数であり、収束速度などを考慮して決定されるが、例えばそれぞれ0.1程度の値とされる。
【数２９】
【００７９】
§２．２ 係数の補正
Ｑ個の係数Aq(ω,k)（１≦ｑ≦Ｑ）の補正は、最小二乗法に基づき、式（３２）に従って新たな係数Aq(ω,k)（１≦ｑ≦Ｑ）を求めることにより行われる。式（３２）で用いるRq→(ω,θq→(ω,k))は§２．１の処理で得られた方向情報θq→(ω,k)と式（６）から得られる。式（３２）にてFq→(ω,k)の計算に用いられる係数Aq(ω,k)（ｑ∈Υ）として、初めて§２．２の処理を行う場合、§１．１の処理で得られた初期値Aini,q(ω,k)を用い、§２．２の処理が初めてではない場合、直前の§２．２の処理で得られた係数Aq(ω,k)（１≦ｑ≦Ｑ）を用いる。
【数３０】
【００８０】
δ回の反復処理が終了した時点で得られている係数Aq(ω,k)と方向情報θq→(ω,k)のＱ個の組み合わせ(Aq(ω,k)，θq→(ω,k))（１≦ｑ≦Ｑ）が｛(Aq,opt(ω,k)，θq,opt→(ω,k))｝q∈{1,…,Q}であり、Ｑ個の反射音情報成分rsq→(ω,k)（１≦ｑ≦Ｑ）となる。すなわち、ｑ番目の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]は式（３３）、式（３４）で与えられる。
【数３１】
【００８１】
以上の過程により、Ｑ個の反射音情報成分rsq→(ω,k)＝[rsAq(ω,k),rsBq(ω,k)]（ｑ＝１，…，Ｑ）が求められる。なお、δ＝１に設定されている場合、係数の補正を行わないことにより、反射音情報として到来方向のみを求めることができる。
【００８２】
《第５実施形態》
第５実施形態は、「§１．１ 方向情報の初期値設定」が第４実施形態と異なる。そこで、第４実施形態と同じ事項については重複説明を省略し、第４実施形態と異なる事項について説明する。
【００８３】
§１．１ 方向情報の初期値設定
最初に、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を設定する。Ｑ個の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定方法として、この実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）を、観測信号X→(ω,k)とテンプレート情報S→(ω)を用いて決定する一般化調和解析的決定方法を説明する。この方法によると、Ｑ個の推定されるべき到来方向それぞれに最も近いと考えられるＱ個の方向に対応するＱ個のテンプレートを決定し、この決定されたＱ個のテンプレートに対応する方向情報をＱ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）とすればよい。なお、この場合、ＰとＱとの間にＱ＜Ｐなる関係がある。
【００８４】
そこで、テンプレート情報の中から上述のようなテンプレートを決定するために、便宜上、ｑ番目の反射音をAq(ω,k,g(ω,q))Sg(ω,q)→(ω)と表すことにする。ここで、g(ω,q)は、テンプレート情報の中でｑ番目の反射音を最も精度良く表現できるテンプレートのインデックスを表す。反射音を構成する係数Aq(ω,k,g(ω,q))は、音源２００自身が持つ位相や壁での反射、距離による減衰などによるテンプレートSg(ω,q)→(ω)と反射音との相違を表す。この場合、観測信号から１番目からｑ番目までのｑ個の反射音を除去して得られる残差信号Eq+1→(ω,k)は式（３５）のように表される。ただし、１≦ｑ≦Ｑであり、E1→(ω,k)=X→(ω,k)である。
【数３２】
【００８５】
ｑ番目の反射音Aq(ω,k,g(ω,q))Sg(ω,q)→(ω)は、式（３５）に基づく残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小とする基準に従って推定される。推定方法は様々あるが、そのうちの一つの方法について述べる。反射音は、Aq(ω,k,g(ω,q))とSg(ω,q)→(ω)の２つの要素で構成されるので、２つの要素に対して最適化することが必要となる。後述の＜処理１＞と＜処理２＞はｑの昇順に各ｑについて行われる。
【００８６】
＜処理１＞
記号Λはインデックスｐの全体の集合{1,…,p,…,P}から後述する式（３７）により決定されたインデックスの集合を除いた集合である。つまり、Λ＝{1,…,p,…,P}-{g(ω,1),…,g(ω,q-1)}とする。ただし、初めて＜処理１＞を行うときはΛ＝{1,…,p,…,P}である。
ｐ番目（ｐ∈Λ）のテンプレートSp→(ω)が式（３５）に基づく残差信号Eq+1→(ω,k)のパワー(Eq+1→(ω,k))HEq+1→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)は、最小二乗法に基づき、式（３６）により求められる。なお、この段階では、式（３６）左辺のｑは意味を持たないことに留意されたい。
【数３３】
【００８７】
＜処理２＞
集合Λの要素の個数（濃度）を|Λ|とすると、式（３６）に基づき得られた|Λ|個の係数Aq(ω,k,p)（ｐ∈Λ）の中から、その絶対値が最大の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（３７）参照）。
【数３４】
【００８８】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（３７）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００８９】
第５実施形態では、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）の決定に一般化調和解析的決定方法を用いた。一般化調和解析的決定方法によると、初期値計算量がビーム探索的決定方法に比べて増大するが、推定されるべき反射音（正解）に近い初期値を設定できる可能性が高く、この場合、推定精度の向上を望めるだけでなく、反復処理の回数を減らすことができる。
【００９０】
《第６実施形態》
第６実施形態では、「§１．１ 方向情報の初期値設定」が第４実施形態と異なる。そこで、第４実施形態と同じ事項については重複説明を省略し、第４実施形態と異なる事項について説明する。第６実施形態のコンセプトは、第４実施形態よりも簡便でありながら、上述の初期値の偏りを防止することにある。
【００９１】
§１．１ 方向情報の初期値設定
まず、方向情報を構成するθpol(ω)とθazi(ω)のいずれか一方について、偏りが無いように複数の方向を決定する。この例では、極角θpol(ω)について偏りが無いようにβ個（β≧2）の方向{θ1,pol(ω),…,θβ,pol(ω)}を決定する。通常、極角は0°≦θpol(ω)≦180°を満たすから、例えば等間隔10°ごとに方向を定めることにより{θ1,pol(ω),…,θβ,pol(ω)}={0,10,20,…,180}となる（β=19）。テンプレート情報に含まれるテンプレートのうち集合{θ1,pol(ω),…,θβ,pol(ω)}の要素のいずれかを極角θpol(ω)として持つ方向情報（位置）に対応するテンプレートのインデックスの集合をΨとする。この際、集合Ψの要素の個数（濃度）|Ψ|が、Ｑ≦|Ψ|＜Ｐを満たすようになることが好ましい。集合Ψは、テンプレート情報に含まれるテンプレートに対応する方向情報（位置）の集合の真部分集合である。
【００９２】
そして、ｐ番目のテンプレートSp→(ω)が残差信号Eq→(ω,k)のパワー(Eq→(ω,k))HEq→(ω,k)を最小化するための最適なテンプレートであると仮定した場合の係数Aq(ω,k,p)を、最小二乗法に基づき、式（３８ａ）により求める。ただし、式（３８ａ）の右辺で用いるテンプレートのインデックスｐはｐ∈Ψである。なお、この段階では、式（３８ａ）左辺のｑは意味を持たないことに留意されたい。
【数３５】
【００９３】
次に、ｐ∈Ψなる条件の下で式（３８ａ）に基づき得られた|Ψ|個の係数Aq(ω,k,p)（１≦ｐ≦|Ψ|）の中から、その絶対値の大きい方から順にＱ個の係数Aq(ω,k,p)のインデックスｑ（１≦ｑ≦Ｑ）を決定する（式（３８ｂ）参照）。記号Γは、集合Ψから式（３８ｂ）により決定されたインデックスの集合を除いた集合であり、Γ＝Ψ-{g(ω,1),…,g(ω,q-1)}である。
【数３６】
【００９４】
従って、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）は、式（３８）により得られるＱ個のg(ω,q)（１≦ｑ≦Ｑ）をインデックスに持つテンプレートSg(ω,q)→(ω)（１≦ｑ≦Ｑ）に対応する方向情報θg(ω,q)→(ω)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）として与えられる。すなわち、θini,q→(ω,k)＝[θg(ω,q),pol(ω),θg(ω,q),azi(ω)]（１≦ｑ≦Ｑ）である。初期値θini,q→(ω,k)はフレームインデックスｋに依存しないことに留意されたい。
【００９５】
これらの実施形態のほか、Ｑ個の方向情報θq→(ω,k)の初期値θini,q→(ω,k)（１≦ｑ≦Ｑ）をランダムに設定する実施形態も許容される。
【００９６】
＜変形例＞
上述の実施形態では周波数ごとに観測信号X→(ω,k)を用いて反射音情報rs→(ω,k)を推定したが、周波数ごとに反射音情報を推定すると、一意に推定されるべき仮想音源の方向（推定到来方向）以外の方向に関する情報も含んでしまうことがあり、この結果、反射音情報に誤差が生じることがありうる。例えば、図１４（ａ）に示すように推定到来方向に関する情報だけを抽出できることが望ましいが、実際には図１４（ｂ）に示すように推定到来方向以外の方向に関する情報が混在してしまうことがありえる。
【００９７】
そこで変形例では、全周波数に亘り一括してパワーを算出することで、反射音情報の推定誤差を小さくする。つまり、図１５に示すように、残差信号のパワーを全周波数に亘り統括することによって、推定到来方向以外の方向の影響を極力減らすことができる。一般的に推定到来方向以外の方向では各周波数でのパワーにバラつきが生じるので、残差信号のパワーを全周波数に亘り統括することにより、推定到来方向のパワーに比してそれ以外の方向のパワーの相対的な影響を低減することができる。なお、図１５では、縦軸のパワーは相対値を示しているので各グラフのスケールが同じであるわけではないことに留意されたい。
【００９８】
この変形例での処理は次のとおりである。解析する周波数帯域に含まれる周波数のインデックスωの集合をΩとする。例えば、音声信号を扱うのであれば、1.0〜3.0kHz帯域に対応するインデックスの集合をΩとすればよい。そして、テンプレートSg(q,k)→(ω)のインデックスg(q,k)を式（１０）の替わりに式（３９）によって求める。
【数３７】
【００９９】
テンプレートSg(ω,q)→(ω)のインデックスg(ω,q)を式（１８）の替わりに式（４０）によって求める。また、方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（１９）の替わりに式（４１）による更新によって行われる。
【数３８】
【０１００】
テンプレートSg(ω,q)→(ω)のインデックスg(ω,q)を式（２８）や式（３７）の替わりに式（４２）によって求める。同様に、式（３８）の替わりに式（４３）によって求める。また、方向情報θq→(ω,k)＝[θq,pol(ω,k),θq,azi(ω,k)]の補正は、式（３１）の替わりに式（４４）による更新によって行われる。
【数３９】
【実施例】
【０１０１】
本発明を用いて相対的な到来時刻差を推定した実験結果を示す。直方体の部屋で１００本のマイクロホンを１０行１０列に等間隔に２次元的に並べ、壁際に配置した環境をシミュレートした。マイクロホン間隔ｄは４ｃｍである。詳細な実験条件は図１１に示してある。テンプレートの座標として、図１０に示したように、マイクロホンアレーを取り囲む半球上に等間隔となるように配置した。1.0-3.0kHz帯域を解析した場合の実験結果を図１２に示す。基準反射音として直接音を選択した。図が煩雑になることを避けるため、ここでは反射音１−３についてのみ図示している。反射音１について、直接音からの到来時刻差は周波数平均で1.4msecであり、反射音２について、直接音からの到来時刻差は周波数平均で2.5msecであり、反射音３について、直接音からの到来時刻差は周波数平均で4.7msecであった。低周波数帯域では距離に対する位相の変化が小さいために到来時刻差の推定値の誤差が大きくなる傾向があるため、周波数帯域全体に亘って平均をとることによって到来時刻差を求めることがよい。
【０１０２】
＜応用例＞
反射音情報は人間が生活する上で、非常に重要な音声情報である。例えば、視覚障害者は、タッピングによって発した音源信号が壁や天井等で反射して耳で観測することにより、環境を把握している。また、日常会話でも、適度な反射が生じる部屋で会話することと、反射音が比較的少ない環境で会話することでは会話のしやすさに相違が生じる。以下、本発明により推定された反射音情報を用いたサービス例について述べる。
１つ目は、会議システムに本発明を組み込んだ例である。指向性音源の向きに応じて反射音の振幅は変化するので、反射音情報が分かると、どの方向に音源が向いているのかを推定することができる。会議システムに音源向きの推定装置を組み込めば、誰に向かって発言したのかを提示することに応用できる。
２つ目は、自由な位置で映像や音声を鑑賞できるシステムである。遠方にある音は直接到来する音源のパワーが小さいので収音することが困難である。反射音情報が分かると、直接音だけでなく、反射音も強調収音できるので、遠方の音を強調することが可能となる。また、音声処理の分野では、方向別に音源の強調収音は可能であるが、距離別に音声を強調収音することは非常に難しいとされている。反射音情報が分かると、距離に対応する物理的な特徴量が得られるので、距離別に収音することが可能となる。遠方の音を収音したり、方向別、距離別に収音することができれば、視聴者の選択した位置に対応した音場を擬似的に生成することが可能となる。
【０１０３】
音声通信システムにおいて、反射音情報を推定することは、直接音だけでは得られなかった音場の情報を得ることにつながる。反射音情報が分かれば、これまでの音声強調技術ではできなかったような遠方音の収音や距離別の収音に結びついたり、従来の収音技術では推定できなかった音場の情報（例えば音源の向き）を推定できる。こういった音場の情報の推定は、これまでの技術では実現できなかった音声処理装置の開発に繋がる。反射音情報の推定に関する従来技術は、インパルス応答を求めるために特殊な信号を観測する必要があったが、本発明は音声信号のような一般的な観測信号で反射音情報を得られるという利点を持つ。
【０１０４】
＜反射音情報推定装置のハードウェア構成例＞
上述の実施形態に関わる反射音情報推定装置は、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ＣＰＵ（Central Processing Unit）〔キャッシュメモリなどを備えていてもよい。〕、メモリであるＲＡＭ（Random Access Memory）やＲＯＭ（Read Only Memory）と、ハードディスクである外部記憶装置、並びにこれらの入力部、出力部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置間のデータのやり取りが可能なように接続するバスなどを備えている。また必要に応じて、反射音情報推定装置に、ＣＤ−ＲＯＭなどの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【０１０５】
反射音情報推定装置の外部記憶装置には、反射音情報を推定するためのプログラム並びにこのプログラムの処理において必要となるデータなどが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。以下、データやその格納領域のアドレスなどを記憶する記憶装置を単に「記憶部」と呼ぶことにする。
【０１０６】
反射音情報推定装置の記憶部には、アナログ信号に対してＡＤ変換を行うためのプログラム、フレーム分割処理を行うためのプログラム、フレームごとのデジタル信号を周波数領域の観測信号に変換するためのプログラム、テンプレート情報を生成するためのプログラム、周波数領域の観測信号とテンプレート情報を用いて反射音情報を推定するためのプログラムが記憶されている。
【０１０７】
反射音情報推定装置では、記憶部に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭに読み込まれて、ＣＰＵで解釈実行・処理される。この結果、ＣＰＵが所定の機能（ＡＤ変換部、フレーム分割部、周波数領域変換部、テンプレート生成部、反射音情報推定部、到来時刻差推定部）を実現することで反射音情報の推定が実現される。
【０１０８】
＜補記＞
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１０９】
また、上記実施形態において説明したハードウェアエンティティ（反射音情報推定装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１１０】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１１３】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【特許請求の範囲】
【請求項１】
所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（以下、伝達特性関数という）に複素振幅を乗じて表される複数の反射音について、基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する反射音情報推定装置であって、
Ｑを２以上の予め定められた整数として、Ｑ個の上記反射音それぞれの複素振幅（以下、到来振幅という）を入力とし、上記基準反射音以外の反射音（以下、対象反射音という）それぞれについて、上記基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、上記基準反射音に対する対象反射音の相対的な到来時刻差とする到来時刻差推定部
を含む反射音情報推定装置。
【請求項２】
請求項１に記載の反射音情報推定装置において、
上記到来時刻差推定部は、
Ｑ個の上記反射音の中から各上記対象反射音に共通の上記基準反射音を決定する
ことを特徴とする反射音情報推定装置。
【請求項３】
請求項１に記載の反射音情報推定装置において、
上記到来時刻差推定部は、
上記対象反射音ごとにＱ個の上記反射音の中から上記基準反射音を決定する
ことを特徴とする反射音情報推定装置。
【請求項４】
請求項１から請求項３のいずれかに記載の反射音情報推定装置において、
Ｑ個の上記反射音の到来振幅はそれぞれ、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は上記観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものである
ことを特徴とする反射音情報推定装置。
【請求項５】
請求項１から請求項３のいずれかに記載の反射音情報推定装置において、
Ｑ個の上記反射音の到来振幅は、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）から相異なるＱ個の上記位置に対応するＱ個の反射音を減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅である
ことを特徴とする反射音情報推定装置。
【請求項６】
所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（以下、伝達特性関数という）に複素振幅を乗じて表される複数の反射音について、基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する反射音情報推定方法であって、
Ｑを２以上の予め定められた整数として、Ｑ個の上記反射音それぞれの複素振幅（以下、到来振幅という）を用いて、上記基準反射音以外の反射音（以下、対象反射音という）それぞれについて、上記基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、上記基準反射音に対する対象反射音の相対的な到来時刻差とする到来時刻差推定過程
を有する反射音情報推定方法。
【請求項７】
請求項６に記載の反射音情報推定方法において、
上記到来時刻差推定過程では、
Ｑ個の上記反射音の中から各上記対象反射音に共通の上記基準反射音を決定する
ことを特徴とする反射音情報推定方法。
【請求項８】
請求項６に記載の反射音情報推定方法において、
上記到来時刻差推定過程では、
上記対象反射音ごとにＱ個の上記反射音の中から上記基準反射音を決定する
ことを特徴とする反射音情報推定方法。
【請求項９】
請求項６から請求項８のいずれかに記載の反射音情報推定方法において、
Ｑ個の上記反射音の到来振幅はそれぞれ、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は上記観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものである
ことを特徴とする反射音情報推定方法。
【請求項１０】
請求項６から請求項８のいずれかに記載の反射音情報推定方法において、
Ｑ個の上記反射音の到来振幅は、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）から相異なるＱ個の上記位置に対応するＱ個の反射音を減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅である
ことを特徴とする反射音情報推定方法。
【請求項１１】
コンピュータに、請求項６から請求項１０のいずれかに記載の反射音情報推定方法の処理を実行させるためのプログラム。
【請求項１】
所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（以下、伝達特性関数という）に複素振幅を乗じて表される複数の反射音について、基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する反射音情報推定装置であって、
Ｑを２以上の予め定められた整数として、Ｑ個の上記反射音それぞれの複素振幅（以下、到来振幅という）を入力とし、上記基準反射音以外の反射音（以下、対象反射音という）それぞれについて、上記基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、上記基準反射音に対する対象反射音の相対的な到来時刻差とする到来時刻差推定部
を含む反射音情報推定装置。
【請求項２】
請求項１に記載の反射音情報推定装置において、
上記到来時刻差推定部は、
Ｑ個の上記反射音の中から各上記対象反射音に共通の上記基準反射音を決定する
ことを特徴とする反射音情報推定装置。
【請求項３】
請求項１に記載の反射音情報推定装置において、
上記到来時刻差推定部は、
上記対象反射音ごとにＱ個の上記反射音の中から上記基準反射音を決定する
ことを特徴とする反射音情報推定装置。
【請求項４】
請求項１から請求項３のいずれかに記載の反射音情報推定装置において、
Ｑ個の上記反射音の到来振幅はそれぞれ、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は上記観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものである
ことを特徴とする反射音情報推定装置。
【請求項５】
請求項１から請求項３のいずれかに記載の反射音情報推定装置において、
Ｑ個の上記反射音の到来振幅は、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）から相異なるＱ個の上記位置に対応するＱ個の反射音を減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅である
ことを特徴とする反射音情報推定装置。
【請求項６】
所定の周波数帯域に含まれる周波数ごとに空間中の任意の位置と複数のマイクロホンとの間の伝達特性を模擬した関数（以下、伝達特性関数という）に複素振幅を乗じて表される複数の反射音について、基準となる反射音に対する他の反射音の相対的な到来時刻差を推定する反射音情報推定方法であって、
Ｑを２以上の予め定められた整数として、Ｑ個の上記反射音それぞれの複素振幅（以下、到来振幅という）を用いて、上記基準反射音以外の反射音（以下、対象反射音という）それぞれについて、上記基準反射音の到来振幅に対する対象反射音の到来振幅の偏角を周波数で除したものに対する周波数についての相加平均を、上記基準反射音に対する対象反射音の相対的な到来時刻差とする到来時刻差推定過程
を有する反射音情報推定方法。
【請求項７】
請求項６に記載の反射音情報推定方法において、
上記到来時刻差推定過程では、
Ｑ個の上記反射音の中から各上記対象反射音に共通の上記基準反射音を決定する
ことを特徴とする反射音情報推定方法。
【請求項８】
請求項６に記載の反射音情報推定方法において、
上記到来時刻差推定過程では、
上記対象反射音ごとにＱ個の上記反射音の中から上記基準反射音を決定する
ことを特徴とする反射音情報推定方法。
【請求項９】
請求項６から請求項８のいずれかに記載の反射音情報推定方法において、
Ｑ個の上記反射音の到来振幅はそれぞれ、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）を用いて、ｑを１以上Ｑ以下の各整数として、各ｑについて、伝達特性関数に複素振幅を乗じたものをｑ番目の最小の残差信号（ただし、１番目の最小の残差信号は上記観測信号とする）から減じて得られるｑ＋１番目の残差信号のパワーが最小になるように当該伝達特性関数および当該複素振幅を推定し、推定された当該複素振幅をｑ番目の反射音の到来振幅とすることによって得られたものである
ことを特徴とする反射音情報推定方法。
【請求項１０】
請求項６から請求項８のいずれかに記載の反射音情報推定方法において、
Ｑ個の上記反射音の到来振幅は、
音声信号をＭ個のマイクロホンで収音して得られるＭ個の収音信号がそれぞれ周波数領域に変換された信号（以下、観測信号という）から相異なるＱ個の上記位置に対応するＱ個の反射音を減じて得られる残差信号のパワーが最小になるようにＱ個の伝達特性関数およびＱ個の複素振幅を一括して推定することによって得られたＱ個の複素振幅である
ことを特徴とする反射音情報推定方法。
【請求項１１】
コンピュータに、請求項６から請求項１０のいずれかに記載の反射音情報推定方法の処理を実行させるためのプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図８】
【図９】
【図１６】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図８】
【図９】
【図１６】
【公開番号】特開２０１２−３９２７６（Ｐ２０１２−３９２７６Ａ）
【公開日】平成２４年２月２３日（２０１２．２．２３）
【国際特許分類】
【出願番号】特願２０１０−１７６０１９（Ｐ２０１０−１７６０１９）
【出願日】平成２２年８月５日（２０１０．８．５）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ・発行者名 社団法人日本音響学会 ・刊行物名 日本音響学会２０１０年春季研究発表会 〔講演論文集〕 講演要旨・講演論文ＣＤ−ＲＯＭ ・発行年月日 ２０１０年３月１日
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
【公開日】平成２４年２月２３日（２０１２．２．２３）
【国際特許分類】
【出願日】平成２２年８月５日（２０１０．８．５）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ・発行者名 社団法人日本音響学会 ・刊行物名 日本音響学会２０１０年春季研究発表会 〔講演論文集〕 講演要旨・講演論文ＣＤ−ＲＯＭ ・発行年月日 ２０１０年３月１日
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【Ｆターム（参考）】
[ Back to top ]