(57)【要約】【課題】少ない事前情報のもとで雑音を効果的に低減することができる。【解決手段】周波数領域変換部は音響信号をチャネル毎に周波数領域信号に変換し、フィルタ係数算出部は周波数領域変換部が変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタを予め定めた個数のフレームからなる区間毎に少なくとも２組算出し、出力信号算出部は周波数領域変換部が変換した周波数領域信号とフィルタ係数算出部が算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する。【選択図】図２
2014026115.tif 000002
【特許請求の範囲】
【請求項１】
音響信号をチャネル毎に周波数領域信号に変換する周波数領域変換部と、
前記周波数領域変換部が変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出するフィルタ係数算出部と、
前記周波数領域変換部が変換した周波数領域信号と前記フィルタ係数算出部が算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する出力信号算出部と、
を備えることを特徴とする音響信号処理装置。
【請求項２】
前記フィルタ係数算出部において、前記チャネル間の伝達特性の差は位相差であり、前記フィルタは前記位相差に基づく遅延和要素であって、
前記位相差の初期値として、チャネル及びフレーム毎に乱数を設定する初期値設定部を、
備えることを特徴とする請求項１に記載の音響信号処理装置。
【請求項３】
前記初期値設定部において、前記位相差の初期値として設定する乱数は位相領域における乱数であり、
前記フィルタ係数算出部は、前記初期値設定部が設定した初期値を用いて、前記残差の大きさを極小化する遅延和要素を与える位相差を再帰的に算出することを特徴とする請求項２に記載の音響信号処理装置。
【請求項４】
前記少なくとも２組のフィルタを要素とするフィルタ行列を特異値分解して特異ベクトルを算出する特異ベクトル算出部を備え、
前記出力信号算出部は、前記特異ベクトル算出部が算出した特異ベクトルと前記周波数領域信号を要素とする入力信号ベクトルに基づいて前記出力信号を算出することを特徴とする請求項１ないし３のいずれかに記載の音響信号処理装置。
【請求項５】
前記出力信号算出部は、前記特異ベクトル算出部が算出した特異ベクトルのうち最も大きい特異値から降順に予め定めた個数の特異値に各々対応する特異ベクトルに基づいて前記出力信号を算出することを特徴とする請求項４に記載の音響信号処理装置。
【請求項６】
音響信号処理装置における音響信号処理方法であって、
前記音響信号処理装置は、音響信号をチャネル毎に周波数領域信号に変換する第１の過程と、
前記第１の過程で変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出する第２の過程と、
前記第１の過程で変換した周波数領域信号と前記第２の過程で算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する第３の過程と、
を有することを特徴とする音響信号処理方法。
【請求項７】
音響信号処理装置のコンピュータに、
音響信号をチャネル毎に周波数領域信号に変換する第１の手順、
前記第１の手順で変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出する第２の手順、
前記第１の手順で変換した周波数領域信号と前記第２の手順で算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する第３の手順、
を実行させるための音響信号処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号処理装置、音響信号処理方法、及び音響信号処理プログラムに関する。
【背景技術】
【０００２】
収録した音響信号から、ある音源による成分、その他の音源による成分、及び雑音による成分を分離する音源分離技術が提案されている。例えば、特許文献１に記載の音源方向推定装置では、消去又は集束させる音を選択するために、当該音源方向推定装置は、音響信号を入力する音響信号入力手段を備え、入力された音響信号の相関行列を算出する。このような音源分離技術では、音源からマイクロホンまでの伝達特性を事前に高い精度で同定しておかなければ、一定の分離精度を得ることができなかった。
しかし、実環境において伝達関数を高い精度で同定することは現実的に困難であった。また、音源分離技術は、人型ロボットが周囲の音声を収録する際、動作中に発生する雑音（例えば、モータの動作音等）の除去することに応用されることが期待されている。しかし、動作中に雑音のみを同定することは困難であった。
そこで、予め設定すべき事前情報が少ない能動雑音制御技術（Ａｃｔｉｖｅ  Ｎｏｉｓｅ  Ｃｏｎｔｒｏｌ、ＡＮＣ）が提案されている。ＡＮＣは、適応フィルタを用いて雑音に対して位相が反転した逆位相波を用いて雑音を低減する技術である。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１０−２８１８１６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ＡＮＣには、適応フィルタを動作して得られるフィルタ係数は、必ずしも大局的な最適解にはならず雑音だけではなく目的音まで抑圧するという問題があった。
本発明は上記の点に鑑みてなされたものであり、少ない事前情報のもとで雑音を効果的に低減する音響信号処理装置、音響信号処理方法、及び音響信号処理プログラムを提供する。
【課題を解決するための手段】
【０００５】
（１）本発明は上記の課題を解決するためになされたものであり、本発明の一態様は、音響信号をチャネル毎に周波数領域信号に変換する周波数領域変換部と、前記周波数領域変換部が変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出するフィルタ係数算出部と、前記周波数領域変換部が変換した周波数領域信号と前記フィルタ係数算出部が算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する出力信号算出部と、を備えることを特徴とする音響信号処理装置である。
【０００６】
（２）本発明のその他の態様は、前記フィルタ係数算出部において、前記チャネル間の伝達特性の差は位相差であり、前記フィルタは前記位相差に基づく遅延和要素であって、前記位相差の初期値として、チャネル及びフレーム毎に乱数を設定する初期値設定部を、備えることを特徴とする。
【０００７】
（３）本発明のその他の態様は、前記初期値設定部において、前記位相差の初期値として設定する乱数は位相領域における乱数であり、前記フィルタ係数算出部は、前記初期値設定部が設定した初期値を用いて、前記残差の大きさを極小化する遅延和要素を与える位相差を再帰的に算出することを特徴とする。
【０００８】
（４）本発明のその他の態様は、前記少なくとも２組のフィルタを要素とするフィルタ行列を特異値分解して特異ベクトルを算出する特異ベクトル算出部を備え、前記出力信号算出部は、前記特異ベクトル算出部が算出した特異ベクトルと前記周波数領域信号を要素とする入力信号ベクトルに基づいて前記出力信号を算出することを特徴とする。
【０００９】
（５）本発明のその他の態様は、前記出力信号算出部は、前記特異ベクトル算出部が算出した特異ベクトルのうち最も大きい特異値から降順に予め定めた個数の特異値に各々対応する特異ベクトルに基づいて前記出力信号を算出することを特徴とする。
【００１０】
（６）本発明のその他の態様は、音響信号処理装置における音響信号処理方法であって、前記音響信号処理装置は、音響信号をチャネル毎に周波数領域信号に変換する第１の過程と、前記第１の過程で変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出する第２の過程と、前記第１の過程で変換した周波数領域信号と前記第２の過程で算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する第３の過程と、を有することを特徴とする音響信号処理方法である。
【００１１】
（７）本発明のその他の態様は、音響信号処理装置のコンピュータに、音響信号をチャネル毎に周波数領域信号に変換する第１の手順、前記第１の手順で変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、前記音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて算出した残差の大きさが極小化されるように前記フィルタのフィルタ係数を予め定めた個数のフレームからなる区間毎に少なくとも２組算出する第２の手順、前記第１の手順で変換した周波数領域信号と前記第２の手順で算出した少なくとも２組のフィルタのフィルタ係数に基づいて周波数領域の出力信号を算出する第３の手順、を実行させるための音響信号処理プログラムである。
【発明の効果】
【００１２】
本発明の一態様（１）、（６）又は（７）によれば、少ない事前情報のもとで雑音を効果的に低減することができる。
本発明のその他の態様（２）によれば、事前情報を容易に生成でき、フィルタ係数を算出する処理量を低減することができる。
本発明のその他の態様（３）によれば、雑音を低減するための遅延和要素についてチャネル間の縮退を回避することができるため、雑音を効果的に低減することができる。
本発明のその他の態様（４）によれば、特定の方向から到来する音波に基づく音響信号に対して雑音を低減することができる。
本発明のその他の態様（５）によれば、特定の方向から到来する雑音をより少ない演算量で顕著に低減することができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の第１の実施形態に係る音響信号処理を示す概念図である。
【図２】本実施形態に係る音響信号処理システムの構成を示す概略図である。
【図３】本実施形態に係る音響信号処理を示すフローチャートである。
【図４】本発明の第２の実施形態に係る音響信号処理システムの構成を示す概略図である。
【図５】本実施形態に係る音響信号処理を示すフローチャートである。
【図６】信号入力部と雑音源及び音源の配置例を示す平面図である。
【図７】信号入力部の構成例を示す概略図である。
【図８】実験に用いた雑音のスペクトルの一例を示す図である。
【図９】実験に用いた目的音のスペクトルの一例を示す図である。
【図１０】繰り返しによる位相の変化の一例を示す図である。
【図１１】特異値の区間数による依存性の一例を示す図である。
【図１２】特異値の区間数による依存性の他の例を示す図である。
【図１３】出力音響信号のスペクトログラムの一例を示す図である。
【図１４】出力音響信号のスペクトログラムの他の例を示す図である。
【図１５】出力音響信号のスペクトログラムのさらに他の例を示す図である。
【図１６】平均ＭＵＳＩＣスペクトルの一例を示す図である。
【図１７】本実施形態に係る方向算出部が定めた音源の方向の一例を示す図である。
【図１８】従来のＭＵＳＩＣ法を用いて推定した音源の方向の一例を示す図である。
【発明を実施するための形態】
【００１４】
（第１の実施形態）
本実施形態に係る音響信号処理は、多チャネルの音響信号をチャネル毎に周波数領域に変換した周波数領域信号に、複数チャネルの信号の遅延和を残差とし、残差の大きさを極小化する遅延和要素から成る遅延和要素行列を算出する。そして遅延和要素行列を特異値分解して得られたユニタリ行列もしくは特異ベクトルを、入力された音響信号に基づく入力信号ベクトルに乗じて出力信号ベクトルを算出する処理である。本音響信号処理では、遅延和要素を算出する際、初期値に乱数を与え残差の大きさを極小化するように再帰的に演算を行う。
【００１５】
そこで、本実施形態に係る音響信号処理の概要について図１を用いて説明する。
図１は、本実施形態に係る音響信号処理を示す概念図である。
図１において、水平方向は時刻を示す。図１の最上行は、あるチャネルの入力音響信号ｙの波形である。チャネル数Ｍは、１よりも大きい予め定めた整数（例えば、８）である。この行において上下方向は振幅を示す。この波形の中央部は、他の区間よりも入力音響信号の振幅が大きく、目的音が主である区間である。この区間の前後の区間は、雑音が主である区間である。
【００１６】
図１の最上行から２行目は、標本化フレーム（ｓａｍｐｌｅｄ  ｆｒａｍｅｓ）の概略を示す図である。標本化フレームとは、フレームｋ毎の周波数領域で表された周波数領域係数ｙｋから抽出（標本化、ｓａｍｐｌｉｎｇ）するフレームである。標本化フレームは、フレーム数Ｌ（Ｌは、０よりも大きい整数）毎に予め定められている。図中に左右方向に、一まとまりに並列している縦棒のそれぞれは、周波数領域係数ｙｋから標本化フレーム毎に抽出された周波数領域係数ｙｋ，ｙｋ＋Ｌ，…を示す。つまり、Ｌフレーム毎に順にｐ個の周波数領域係数ｙｋ，ｙｋ＋Ｌ，…を各チャネルについて抽出される。ｐは、予め定めた整数（図１に示す例では５）である。そして、各チャネルについて周波数領域係数ｙｋを要素としてｐ個含む入力信号行列Ｙｋ１，Ｙｋ２，…をｐ・Ｌフレームの区間毎にＱ（図１に示す例では５）個ずつ生成する。
図１の最上行から３行目の下向きの矢印ｄ１〜ｄ５は、各矢印の起点に示す入力信号行列Ｙｋ１，Ｙｋ２，…に基づいて、残差の大きさを極小化するようなフィルタ係数による遅延要素ベクトルｃｋ１，ｃｋ２，…をそれぞれ算出する遅延和算出処理を示す。遅延要素ベクトルｃｋ１等は、入力信号行列Ｙｋ１等に対してそれぞれチャネル間の位相差を補償する遅延和要素を表すフィルタを与える非零のベクトルである。
【００１７】
図１の最下行の下向きの矢印は、算出した遅延要素ベクトルｃｋ１，ｃｋ２，…を標本化フレーム間で統合して得られる遅延和要素行列Ｃを特異値分解し、ユニタリ行列Ｖｃを算出することを示す。特異値分解において、０又は０よりも予め定めた閾値よりも大きい特異値にそれぞれ対応するＭ’（Ｍ’は、１又は１より大きくＭより小さい整数、例えば５）個の右特異ベクトルｖ１，ｖ２，…，ｖＭ’が算出される。ユニタリ行列Ｖｃは、算出されたＭ’個の右特異ベクトルを対応する特異値の降順に統合した行列［ｖ１，ｖ２，…，ｖＭ’］である。本実施形態では、ユニタリ行列Ｖｃの共役転置行列ＶｃＨを、各チャネルにおける周波数領域係数ｙｋを要素とする入力信号ベクトルｙに乗じて、周波数領域における出力信号ｚｋを要素とする出力信号ベクトルｚが得られる。これにより、Ｍ−Ｍ’個の雑音成分が低減し、それぞれ位置の異なるＭ’個の音源それぞれから到来した周波数領域の信号が抽出される。なお、本実施形態では、図１に示す処理を周波数毎に行う。
【００１８】
（音響信号処理システムの構成）
次に、本実施形態に係る音響信号処理システム１の構成について説明する。
図２は、本実施形態に係る音響信号処理システム１の構成を示す概略図である。
音響信号処理システム１は、信号入力部１１、音響信号処理装置１２及び信号出力部１３を含んで構成される。なお、以下の説明では、特に断らない限りベクトル、行列を［…］と示す。ここで、ベクトルを、例えば［ｙ］と小文字で示し、行列を、例えば［Ｙ］と大文字で示す。
【００１９】
信号入力部１１は、Ｍチャネルの音響信号を取得し、取得したＭチャネルの音響信号を音響信号処理装置１２に出力する。信号入力部１１は、マイクロホンアレイと変換部を備える。マイクロホンアレイは、例えば、各々異なる位置に設置されたＭ個のマイクロホン１１１−１〜１１１−Ｍを含んで構成される。各マイクロホン１１１−１等は、到達した音波を電気信号であるアナログ音響信号に変換して変換部に出力する。変換部は、入力されたアナログ音響信号をＡＤ（Ａｎａｌｏｇ−ｔｏ−Ｄｉｇｉｔａｌ、アナログディジタル）変換してチャネル毎にディジタル音響信号を生成する。変換部は、生成したディジタル信号を音響信号処理装置１２にチャネル毎に出力する。信号入力部１１に係るマイクロホンアレイの構成例については後述する。なお、信号入力部１１は、通信回線を通じて遠隔地の通信機器、又はデータ記憶装置からＭチャネルの音響信号を入力する入力インタフェース部であってもよい。
【００２０】
信号出力部１３は、音響信号処理装置１２が出力したＭ’チャネルの出力音響信号を音響信号処理システム１の外部に出力する。信号出力部１３は、例えば、Ｍ’チャネルのうち任意のチャネルの出力音響信号に基づく音を再生する音響再生部である。また、信号出力部１３は、Ｍ’チャネルの出力音響信号をデータ記憶装置又は通信回線を通じて遠隔地の通信機器に出力する出力インタフェース部であってもよい。
【００２１】
音響信号処理装置１２は、周波数領域変換部１２１、入力信号行列生成部１２２、初期値設定部１２３、遅延和要素行列算出部（フィルタ係数算出部）１２４、特異ベクトル算出部１２５、出力信号ベクトル算出部（出力信号算出部）１２６、及び時間領域変換部１２７を含んで構成される。
【００２２】
周波数領域変換部１２１は、信号入力部１１から入力されたＭチャネルの音響信号を、各チャネルについてフレーム毎に時間領域から周波数領域に変換して周波数領域係数を算出する。周波数領域変換部１２１は、周波数領域への変換において、例えば、高速フーリエ変換（Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ、ＦＦＴ）を用いる。周波数領域変換部１２１は、フレーム毎に算出した周波数領域係数を入力信号行列生成部１２２及び出力信号ベクトル算出部１２６に出力する。なお、入力信号行列生成部１２２、初期値設定部１２３、遅延和要素行列算出部１２４、特異ベクトル算出部１２５及び出力信号ベクトル算出部１２６は、以下に述べる処理を各周波数について行う。
【００２３】
入力信号行列生成部１２２は、周波数領域変換部１２１からフレーム毎に入力されたＭチャネルの周波数領域係数に基づいて入力信号行列［Ｙｋ］を生成する。ここで、入力信号行列生成部１２２は、標本数ｐ、フレーム間隔Ｌを予め設定しておく。入力信号行列生成部１２２は、入力されたチャネルｍ（ｍは、０より大きくＭと等しいかＭより小さい整数）の周波数領域係数ｙｍｋをＬフレーム毎に１回ずつ、ｐ回抽出する。入力信号行列生成部１２２は、抽出した周波数領域係数ｙｍｋを、チャネルｍを行方向に、標本数ｐを列方向に配列してｐ・Ｌフレームからなる各区間についてＭ行Ｌ列の入力信号行列［Ｙｋ］を生成する。従って、入力信号行列［Ｙｋ］は式（１）で表される。
【００２４】
【数１】
2014026115.tif 000003
【００２５】
入力信号行列生成部１２２は、生成した区間毎の入力信号行列［Ｙｋ］を区間毎に遅延和要素行列算出部１２４に出力する。
なお、入力信号行列生成部１２２はＬフレーム毎に周波数領域係数ｙｍｋを抽出せずに、１フレーム毎に周波数領域係数ｙｍｋを抽出してもよい。上述のようにＬフレーム毎に周波数領域係数ｙｍｋを抽出した場合、極力異なる時刻に取得された周波数領域係数ｙｍｋを用いて、後述する遅延要素ベクトルの解として、より安定な解を求めることができる。
【００２６】
初期値設定部１２３には、予め定めた区間数Ｑが設定されており、Ｑ個の遅延要素ベクトル［ｃｋ］の初期値を設定する。遅延要素ベクトル［ｃｋ］は、フレームｋにおいて予め定めたチャネル（例えば、チャネル１）と他のチャネルｍとの間の位相差θｍ，ｋを要素として有するベクトルである。遅延要素ベクトル［ｃｋ］は、一般に式（２）で表される。
【００２７】
【数２】
2014026115.tif 000004
【００２８】
式（２）においてωは角周波数である。従って、位相差θｍ，ｋの初期値は、（Ｍ−１）・Q個存在する。
初期値設定部１２３は、（Ｍ−１）・Q個の初期値θｍ，ｋを［−π，π）の範囲の乱数として設定する。望ましい位相角について事前に情報がない場合、この乱数として一様乱数を用いることができるが、この場合、遅延要素ベクトル［ｃｋ］の各要素値（但し、チャネル１を除く）は、単位円上において位相角の方向に一様に分布する乱数、つまり位相角領域での一様乱数となる。
初期値設定部１２３は、設定したＱ個の遅延要素ベクトル［ｃｋ］の初期値を遅延和要素行列算出部１２４に出力する。
【００２９】
遅延和要素行列算出部１２４は、入力信号行列生成部１２２から入力された区間毎の入力信号行列［Ｙｋ］と初期値設定部１２３から入力された区間毎の遅延要素ベクトル［ｃｋ］の初期値に基づいて、遅延要素ベクトル［ｃｋ］を算出する。ここで、遅延和要素行列算出部１２４は、残差ベクトル［εｋ］の大きさであるノルム｜［εｋ］｜が極小化されるように遅延要素ベクトル［ｃｋ］を算出する。残差ベクトル［εｋ］は、入力信号行列［Ｙｋ］に遅延要素ベクトル［ｃｋ］から成る遅延和フィルタを作用して得られるベクトルである。つまり、遅延和要素行列算出部１２４は、遅延和の大きさが零になる方向である死角に対応した遅延要素ベクトル［ｃｋ］を求める。言い換えれば、遅延要素ベクトル［ｃｋ］は死角制御型ビームフォーマを要素として有するベクトルである。また、遅延要素ベクトル［ｃｋ］は、各チャネルの周波数領域係数ｙｍｋにそれぞれ乗じられる係数を有するフィルタ係数群とみることができる。
【００３０】
遅延和要素行列算出部１２４は、ノルム｜［εｋ］｜が極小化される遅延要素ベクトル［ｃｋ］を算出するために、例えば、最小平均二乗（Ｌｅａｓｔ  Ｍｅａｎ  Ｓｑｕａｒｅ）法等の既知の方法を用いる。例えば、遅延和要素行列算出部１２４は、最小平均二乗法を用いて、式（３）に示すように、現在の繰り返し（ｉｔｅｒａｔｉｏｎ）ｔにおける位相θｍ，ｋ（ｔ）に基づいて、次の繰り返しｔ＋１における位相θｍ，ｋ（ｔ＋１）を再帰的に算出する。
【００３１】
【数３】
2014026115.tif 000005
【００３２】
式（３）において、［θｋ（ｔ＋１）］は、繰り返しｔ＋１におけるフレームｋに係る各チャネルｍの位相θｍ，ｋを要素とするベクトルである。αは、予め定めた正の実数（例えば、０．０００１２）である。式（３）を用いて位相θｍ，ｋ（ｔ＋１）を算出する方法は勾配法と呼ばれる方法である。
遅延和要素行列算出部１２４は、Ｑ個の区間毎に算出した遅延要素ベクトル［ｃｋ］を区間の順序に行方向に配置してＱ行Ｍ列の遅延和要素行列［Ｃ］を生成する。
遅延和要素行列算出部１２４は、Ｑ個の区間毎に生成した遅延和要素行列［Ｃ］を特異ベクトル算出部１２５に出力する。
【００３３】
上述したように、初期値設定部１２３では、位相差θｍ，ｋの初期値に乱数を与え、与えられた位相差θｍ，ｋの初期値に基づいて複数の遅延要素ベクトル［ｃｋ］の初期値が得られる。遅延和要素行列算出部１２４は、複数の遅延要素ベクトル［ｃｋ］のそれぞれについて残差を極小化する解の候補を算出する。これらの遅延要素ベクトル［ｃｋ］を算出するために用いる入力信号行列［Ｙｋ］は、それぞれ異なる時間の区間毎に入力された音響信号に基づく。本実施形態では、上述のように初期値に乱数を与え再帰的に位相差を算出する処理法をモンテカルロ的パラメータ探索法と呼ぶ。
このように、初期値に乱数を与えることで複数の遅延要素ベクトル［ｃｋ］を縮退することなく算出することができるため特定の方向から到来する雑音を抑圧するベクトル空間を表すための十分な解が得られる。また、雑音は定常的に発生するのに対し、人間の発話等の目的音は一時的に発生する傾向がある。上述のように、複数の区間にわたって算出した遅延要素ベクトル［ｃｋ］には、雑音のみが到来する区間において算出されたものが主であり、目的音と雑音の両者が到来する区間において算出されたものは比較的少ない。言い換えれば、目的音を抑圧する遅延要素ベクトル［ｃｋ］は、ごく一部に限られる。
【００３４】
特異ベクトル算出部１２５は、遅延和要素行列算出部１２４からＱ個の区間毎に入力された遅延和要素行列［Ｃ］を特異値分解してＱ行Ｍ列の特異値行列［Σ］を算出する。特異値分解とは、特異値行列［Σ］の他、式（４）の関係を満足するＱ行Ｑ列のユニタリ行列ＵとＭ行Ｍ列のユニタリ行列Ｖを算出する演算である。
【００３５】
【数４】
2014026115.tif 000006
【００３６】
式（４）において、［Ｖ］Ｈは、行列［Ｖ］の共役転置行列（ｃｏｎｊｕｇａｔｅ  ｔｒａｎｓｐｏｓｅ  ｍａｔｒｉｘ）である。行列［Ｖ］は、各列に特異値σ１，…，σＭそれぞれに対応するＭ個の右特異ベクトル［v１］，…，［ｖＭ］を有する。順序を示すインデックス１，…，Ｍは、特異値σ１，…，σＭの降順である。特異ベクトル算出部１２５は、このＭ個の右特異ベクトルから、Ｍ’個（Ｍ’は、Ｍと等しいかＭより小さく、０よりも大きい予め定めた整数）の右特異ベクトル［v１］，…，［ｖＭ’］を選択する。これにより、ゼロ又はゼロに近似する特異値に対応する特異ベクトルを排除する。なお、特異ベクトル算出部１２５は、このＭ個の右特異ベクトルから、特異値が予め定めた閾値σｔｈよりも大きい特異値にそれぞれ対応したＭ’個の右特異ベクトル［v１］，…，［ｖＭ’］を選択してもよい。
特異ベクトル算出部１２５は、選択したＭ’個の右特異ベクトル［v１］，…，［ｖＭ’］を特異値の降順に列方向に配列してＭ行Ｍ’列の行列［Ｖｃ］を生成し、生成した行列［Ｖｃ］の共役転置行列［Ｖｃ］Ｈを生成する。特異ベクトル算出部１２５は、生成した共役転置行列［Ｖｃ］ＨをＱ個の区間毎に出力信号ベクトル算出部１２６に出力する。
【００３７】
出力信号ベクトル算出部１２６は、周波数領域変換部１２１からフレーム毎に入力されたＭチャネルの周波数領域係数に基づいて入力信号ベクトル［ｙｋ］を生成する。出力信号ベクトル算出部１２６は、入力された各チャネルｍのフレームｋ毎の周波数領域係数ｙｍｋを列方向に配列してＭ列の入力信号ベクトル［ｙｋ］を生成する。出力信号ベクトル算出部１２６は、生成したＭ列の入力信号ベクトル［ｙｋ］に特異ベクトル算出部１２５から入力されたＭ’行Ｍ列の共役転置行列［Ｖｃ］Ｈを乗算してＭ’列の出力信号ベクトル［ｚｋ］を算出する。各列の成分は、チャネル毎の出力周波数領域係数を示す。即ち、右特異ベクトル［v１］，…，［ｖＭ’］の各々は、入力信号ベクトル［ｙｋ］に対するフィルタ係数とみなすことができる。出力信号ベクトル算出部１２６は、算出した出力信号ベクトル［ｚｋ］を時間領域変換部１２７に出力する。
なお、出力信号ベクトル算出部１２６は、入力信号ベクトル［ｙｋ］に右特異ベクトル［v１］，…，［ｖＭ’］を転置したベクトル［v１］Ｈ，…，［ｖＭ’］Ｈのうちいずれか１個を乗算して出力周波数領域係数ｚｋ（スカラー量）を算出してもよい。出力信号ベクトル算出部１２６は、算出した出力周波数領域係数を時間領域変換部１２７に出力する。ここで、入力信号ベクトル［ｙｋ］に乗算するベクトルとして、最大の特異値σ１に対応するベクトル［v１］Ｈを用いる。共役転置行列［Ｖｃ］Ｈは、雑音成分を最小化する成分を要素として含むベクトル［v１］Ｈ，…，［ｖＭ’］Ｈからなる行列である。特異値σ１，…，σＭ’は、各ベクトル［v１］Ｈ，…，［ｖＭ’］Ｈが遅延和要素行列に寄与する度合いを示すから、最も雑音成分を最小化する成分の比率が高いベクトル［v１］Ｈを用いることにより、雑音を効果的に抑圧することができる。
【００３８】
時間領域変換部１２７は、出力信号ベクトル算出部１２６から入力された出力信号ベクトル［ｚｋ］が有する出力周波数領域係数を、各フレームについてチャネル毎に周波数領域から時間領域に変換して時間領域の出力音響信号を算出する。時間領域変換部１２７は、時間領域への変換において、例えば、高速逆フーリエ変換（Ｉｎｖｅｒｓｅ  Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ、ＩＦＦＴ）を用いる。時間領域変換部１２７は、算出したチャネル毎の出力音響信号を信号出力部１３に出力する。
【００３９】
（音響信号処理）
次に、本実施形態に係る音響信号処理について説明する。
図３は、本実施形態に係る音響信号処理を示すフローチャートである。
（ステップＳ１０１）信号入力部１１は、Ｍチャネルの音響信号を取得し、取得したＭチャネルの音響信号を音響信号処理装置１２に出力する。その後、ステップＳ１０２に進む。
（ステップＳ１０２）周波数領域変換部１２１は、信号入力部１１から入力されたＭチャネルの音響信号を、各チャネルについてフレーム毎に時間領域から周波数領域に変換して周波数領域係数を算出する。周波数領域変換部１２１は、算出した周波数領域係数を入力信号行列生成部１２２及び出力信号ベクトル算出部１２６に出力する。その後、ステップＳ１０３に進む。
【００４０】
（ステップＳ１０３）入力信号行列生成部１２２は、周波数領域変換部１２１からフレーム毎に入力されたＭチャネルの周波数領域係数に基づいてｐ・Ｌフレームからなる各区間について入力信号行列［Ｙｋ］を生成する。入力信号行列生成部１２２は、生成した区間毎の入力信号行列［Ｙｋ］を遅延和要素行列算出部１２４に出力する。その後、ステップＳ１０４に進む。
（ステップＳ１０４）初期値設定部１２３は、（Ｍ−１）・Q個の初期値θｍ，ｋを［−π，π）の範囲で乱数として設定し、各々（Ｍ−１）個の初期値θｍ，ｋに基づいてQ個の遅延要素ベクトル［ｃｋ］の初期値を設定する。初期値設定部１２３は、設定したＱ個の遅延要素ベクトル［ｃｋ］の初期値を遅延和要素行列算出部１２４に出力する。その後、ステップＳ１０５に進む。
【００４１】
（ステップＳ１０５）遅延和要素行列算出部１２４は、入力信号行列生成部１２２から入力された入力信号行列［Ｙｋ］と初期値設定部１２３から入力された区間毎の遅延要素ベクトル［ｃｋ］の初期値に基づいて、遅延要素ベクトル［ｃｋ］を算出する。ここで、遅延和要素行列算出部１２４は、残差ベクトル［εｋ］のノルム｜［εｋ］｜が極小化されるように遅延要素ベクトル［ｃｋ］を算出する。遅延和要素行列算出部１２４は、Ｑ個の遅延要素ベクトル［ｃｋ］を順に行方向に配置して遅延和要素行列［Ｃ］を生成する。遅延和要素行列算出部１２４は、生成した遅延和要素行列［Ｃ］を特異ベクトル算出部１２５に出力する。その後、ステップＳ１０６に進む。
（ステップＳ１０６）特異ベクトル算出部１２５は、遅延和要素行列算出部１２４から入力された遅延和要素行列［Ｃ］を特異値分解して、特異値行列［Σ］、ユニタリ行列Ｕ及びユニタリ行列Ｖを算出する。特異ベクトル算出部１２５は、特異値σ１，…，σＭの降順にユニタリ行列Ｖから選んだＭ’個の右特異ベクトル［v１］，…，［ｖＭ’］を列方向に配列して行列［Ｖｃ］を生成する。特異ベクトル算出部１２５は、生成した行列［Ｖｃ］の共役転置行列［Ｖｃ］Ｈを出力信号ベクトル算出部１２６に出力する。その後、ステップＳ１０７に進む。
【００４２】
（ステップＳ１０７）出力信号ベクトル算出部１２６は、周波数領域変換部１２１からフレーム毎に入力されたＭチャネルの周波数領域係数に基づいて入力信号ベクトル［ｙｋ］を生成する。出力信号ベクトル算出部１２６は、生成した入力信号ベクトル［ｙｋ］に特異ベクトル算出部１２５から入力されたＭ’行Ｍ列の共役転置行列［Ｖｃ］Ｈを乗算してＭ’列の出力信号ベクトル［ｚｋ］を算出する。出力信号ベクトル算出部１２６は、算出した出力信号ベクトル［ｚｋ］を時間領域変換部１２７に出力する。その後、ステップＳ１０８に進む。
（ステップＳ１０８）時間領域変換部１２７は、出力信号ベクトル算出部１２６から入力された出力信号ベクトル［ｚｋ］が有する出力周波数領域係数を、各フレームについてチャネル毎に周波数領域から時間領域に変換して時間領域の出力音響信号を算出する。時間領域変換部１２７は、算出したチャネル毎の出力音響信号を信号出力部１３に出力する。その後、ステップＳ１０９に進む。
（ステップＳ１０９）信号出力部１３は、音響信号処理装置１２が出力したＭ’チャネルの出力音響信号を音響信号処理システム１の外部に出力する。その後、処理を終了する。
【００４３】
以上に、説明したように本実施形態では、音響信号をチャネル毎に周波数領域信号に変換する。本実施形態では、変換した周波数領域信号をフレーム毎に標本化した標本化信号に対して、遅延要素を並べたベクトル（遅延要素ベクトル）によって表現される音響信号のチャネル間の伝達特性の差を補償するフィルタに基づいて、算出した残差が極小化されるように前記遅延要素ベクトルを予め定めた個数のフレームの区間毎に少なくとも２組算出する。また、本実施形態では、変換した周波数領域信号と算出した少なくとも２組のフィルタ係数に基づいて周波数領域の出力信号を算出する。これにより、本実施形態では特定の方向から到来する雑音を極小化するフィルタが算出されるため、算出されたフィルタに基づいて、その方向から到来する雑音が抑制される。従って、少ない事前情報のもとで雑音を効果的に低減することができる。
【００４４】
また、本実施形態は、さらにチャネル間の伝達特性の差は位相差であり、フィルタは位相差に基づく遅延和であって、位相差の初期値として、チャネル及び予め定めた時間毎に位相領域での乱数を設定する。これにより、事前情報である位相差の初期値が容易に生成され、フィルタ係数を算出する処理量を低減することができる。
また、本実施形態は、少なくとも２組の遅延要素ベクトルを要素とする遅延和要素行列を特異値分解して特異ベクトルを算出し、算出した特異ベクトルと周波数領域信号を要素とする入力信号ベクトルに基づいて出力信号を算出する。本実施形態において、特異値分解の対象である遅延和要素行列は、入力信号ベクトルの雑音成分が極小化される遅延和要素に相当する要素ベクトルから成るので、算出された特異ベクトルと入力信号ベクトルの雑音成分はほぼ直交する。そのため、本実施形態によれば、特定の方向から到来する音波に基づく音響信号に対して雑音を低減することができる。
また、本実施形態は、算出した特異ベクトルのうち最も大きい特異値から降順に予め定めた個数の特異値に各々対応する特異ベクトルに基づいて前記出力信号を算出する。特異値は雑音成分を最小化する成分の比率を示すので、本実施形態によれば、特定の方向から到来する雑音をより少ない演算量で低減することができる。
【００４５】
（第２の実施形態）
次に本発明の第２の実施形態について説明する。
本実施形態に係る音響信号処理システム２の構成について同一の構成及び処理について同一の符号を付して説明する。
図４は、本実施形態に係る音響信号処理システム２の構成を示す概略図である。
音響信号処理システム２は、信号入力部１１、音響信号処理装置２２、信号出力部１３、及び方向出力部２３を含んで構成される。
音響信号処理装置２２は、周波数領域変換部１２１、入力信号行列生成部１２２、初期値設定部１２３、遅延和要素行列算出部１２４、特異ベクトル算出部１２５、出力信号ベクトル算出部１２６、時間領域変換部１２７に加え、方向推定部２２１を含んで構成される。
【００４６】
方向推定部２２１は、出力信号ベクトル算出部１２６が出力した出力信号ベクトル［ｚｋ］に基づいて音源の方向を推定し、推定した音源の方向を示す音源方向信号を方向出力部２３に出力する。方向推定部２２１は、音源の方向を推定する際、例えばＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ  Ｓｉｇｎａｌ  Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法を用いる。ＭＵＳＩＣ法は、雑音部分空間と信号部分空間が直交することを利用して音波の到来方向を音源の方向として推定する方法である。
ＭＵＳＩＣ法を用いる場合には、方向推定部２２１は、相関行列算出部２２１１、固有ベクトル算出部２２１２及び方向算出部２２１３を含んで構成される。相関行列算出部２２１１、固有ベクトル算出部２２１２及び方向算出部２２１３は、特に断らない限り周波数毎に処理を行う。
【００４７】
出力信号ベクトル算出部１２６は、出力信号ベクトル［ｚｋ］を相関行列算出部２２１１にも出力する。相関行列算出部２２１１は、出力信号ベクトル［ｚｋ］に基づいて式（５）を用いてＭ’行Ｍ’列の相関行列［Ｒｚｚ］を算出する。
【００４８】
【数５】
2014026115.tif 000007
【００４９】
即ち、この相関行列［Ｒｚｚ］は、チャネル間における出力信号値の積についての、予め定めたフレーム数にわたる時間平均値を要素とする行列である。相関行列算出部２２１１は、算出した相関行列［Ｒｚｚ］を固有ベクトル算出部２２１２に出力する。
【００５０】
固有ベクトル算出部２２１２は、相関行列算出部２２１１から入力された相関行列［Ｒｚｚ］を対角化してＭ’個の固有ベクトル［ｆ１］，…，［ｆＭ’ ］を算出する。固有ベクトル［ｆ１，…，［ｆＭ’ ］の順序は、それぞれ対応する固有値λ１，…，λＭ’の降順である。固有ベクトル算出部２２１２は、算出した固有ベクトル［ｆ１］，…，［ｆＭ’ ］を方向算出部２２１３に出力する。
【００５１】
方向算出部２２１３には、固有ベクトル算出部２２１２から固有ベクトル［ｆ１］，…，［ｆＭ’ ］が入力され、特異ベクトル算出部１２５から共役転置行列［Ｖｃ］Ｈが入力される。方向算出部２２１３は、ステアリングベクトル［ａ（φ）］を生成する。ステアリングベクトル［ａ（φ）］は、信号入力部１１が備えるマイクロホン１１１−１〜１１１−Ｍの代表点（例えば、重心点）から方向φにある音源からマイクロホン１１１−１〜１１１−Ｍの各々までの音波の伝達特性を表す係数を要素として有するベクトルである。ステアリングベクトル［ａ（φ）］は、例えば、［ａ１（φ），…，ａＭ（φ）］Ｈである。本実施形態では、係数ａ１（φ）〜ａＭ（φ）は、例えば、方向φにある音源からマイクロホン１１１−１〜１１１−Ｍの各々までの伝達関数を示す。そのために、方向算出部２２１３は、予め方向φと伝達関数ａ１（φ），…，ａＭ（φ）を対応付けて記憶させておいた記憶部を備える。
係数ａ１（φ）〜ａＭ（φ）は、方向φから到来する音波に対するチャネル間の位相差を示す大きさ１の係数であってもよい。例えば、マイクロホン１１１−１〜１１１−Ｍが一直線上に配列されており、方向φがその配列方向を基準とする角度である場合、係数ａｍ（φ）は、ｅｘｐ（−ｊωｄｍ，１ｓｉｎφ）である。ｄｍ，１は、マイクロホン１１１−ｍとマイクロホン１１１−１との間の距離である。従って、方向算出部２２１３は、マイクロホン間距離ｄｍ，１を予め設定しておけば、任意のステアリングベクトル［ａ（φ）］を算出することができる。
【００５２】
方向算出部２２１３は、算出したステアリングベクトル［ａ（φ）］、入力された共役転置行列［Ｖｃ］Ｈ及び固有ベクトル［ｆ１］，…，［ｆＭ’ ］に基づき、式（６）を用いてＭＵＳＩＣスペクトルＰ（φ）を周波数毎に算出する。
【００５３】
【数６】
2014026115.tif 000008
【００５４】
式（６）において、Ｍ’’は、推定対象となる音源の最大数を示す整数であって、０より大きくＭ’よりも小さい整数である。これにより、方向算出部２２１３は、算出したＭＵＳＩＣスペクトルＰ（φ）を予め設定した周波数帯域内で平均して平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)を算出する。予め設定した周波数帯域として、発話者が発する音声の音圧が大きい周波数帯域であり、かつ雑音の音圧が小さい周波数帯域を用いてもよい。例えば、周波数帯域は０．５〜２．８ｋＨｚである。
【００５５】
方向算出部２２１３は、算出したＭＵＳＩＣスペクトルＰ(φ)を広帯域信号に拡張して平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)を算出してもよい。そのために、方向算出部２２１３は、出力信号ベクトル算出部１２６から入力された出力信号ベクトルに基づいて、予め設定した閾値よりもＳ／Ｎ比が高い（即ち、ノイズが少ない）周波数ωを選択する。方向算出部２２１３は、選択した周波数ωにおいて固有ベクトル算出部２２１２が算出した最大固有値λ１の平方根にＭＵＳＩＣスペクトルＰ(φ)に対して、式（７）を用いて重み付け加算して広帯域のＭＵＳＩＣスペクトルＰａｖｇ（φ)を算出する。
【００５６】
【数７】
2014026115.tif 000009
【００５７】
式（４）において、Ωは周波数ωの集合を示し、｜Ω｜は集合Ωの要素数、ｋは周波数帯域を示すインデックスを示す。重み付け加算によって平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)には、周波数帯域ωにおけるＭＵＳＩＣスペクトルＰａｖｇ（φ)による成分が強く反映される。
方向算出部２２１３は、平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)のピーク値（極大値）を検知し、検知したピーク値に対応する方向φを最大Ｍ’’個選択する。この選択されたφが音源方向として推定される。
方向算出部２２１３は、選択した方向φを示す方向情報を方向出力部２３に出力する。
【００５８】
方向出力部２３は、方向算出部２２１３から入力された方向情報を音響信号処理システム２の外部に出力する。方向出力部２３は、方向情報をデータ記憶装置又は通信回線を通じて遠隔地の通信機器に出力する出力インタフェース部であってもよい。
【００５９】
（音響信号処理）
次に、本実施形態に係る音響信号処理について説明する。
図５は、本実施形態に係る音響信号処理を示すフローチャートである。
図５に示す音響信号処理は、図３に示すステップＳ１０１〜Ｓ１０９に対してステップＳ２０１〜Ｓ２０４が加わっている。本実施形態では、ステップＳ１０８〜Ｓ１０９を実行した後でステップＳ２０１〜Ｓ２０４を実行してもよいが、これには限られない。本実施形態では、ステップＳ１０８〜Ｓ１０９とステップＳ２０１〜Ｓ２０４とを並行して実行してもよいし、ステップＳ１０８〜Ｓ１０９をステップＳ２０１〜Ｓ２０４の後で実行してもよい。以下、ステップＳ１０８〜Ｓ１０９の後で、ステップＳ２０１〜Ｓ２０４を実行する場合を例にとって説明する。
【００６０】
（ステップＳ２０１）相関行列算出部２２１１は、出力信号ベクトル算出部が算出した出力信号ベクトル［ｚｋ］に基づいて、式（５）を用いてＭ’行Ｍ’列の相関行列［Ｒｚｚ］を算出する。相関行列算出部２２１１は、算出した相関行列［Ｒｚｚ］を固有ベクトル算出部２２１２に出力する。その後、ステップＳ２０２に進む。
（ステップＳ２０２）固有ベクトル算出部２２１２は、相関行列算出部２２１１から入力された相関行列［Ｒｚｚ］を対角化してＭ’個の固有ベクトル［ｆ１］，…，［ｆＭ’ ］を算出する。固有ベクトル算出部２２１２は、算出した固有ベクトル［ｆ１］，…，［ｆＭ’ ］を方向算出部２２１３に出力する。その後、ステップＳ２０３に進む。
【００６１】
（ステップＳ２０３）方向算出部２２１３は、ステアリングベクトル［ａ（φ）］を生成する。方向算出部２２１３は、生成したステアリングベクトル［ａ（φ）］、固有ベクトル算出部２２１２から入力された固有ベクトル［ｆ１］，…，［ｆＭ’ ］、及び特異ベクトル算出部１２５から入力された共役転置行列［Ｖｃ］Ｈに基づき、式（６）を用いてＭＵＳＩＣスペクトルＰ（φ）を周波数毎に算出する。方向算出部２２１３は、算出したＭＵＳＩＣスペクトルＰ（φ）を予め設定した周波数帯域内で平均して平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)を算出する。
方向算出部２２１３は、平均ＭＵＳＩＣスペクトルＰａｖｇ（φ)のピーク値を検知し、検知したピーク値に対応する方向φを定め、定めた方向φを示す方向情報を方向出力部２３に出力する。その後、ステップＳ２０４に進む。
（ステップＳ２０４）方向出力部２３は、方向算出部２２１３から入力された方向情報を音響信号処理システム２の外部に出力する。その後、処理を終了する。
【００６２】
（実験例）
次に、本実施形態に係る音響信号処理システム２を動作させて行った実験例について説明する。実験では、実験室内に配置された１個の雑音源３１に雑音を放射させ、１個の音源３２に目的音を放射させた。収録した雑音と目的音が混合した音響信号を信号入力部１１から入力して音響信号処理システム２を動作させた。
【００６３】
信号入力部１１と雑音源３１及び音源３２の配置例について説明する。
図６は、信号入力部１１と雑音源３１及び音源３２の配置例を示す平面図である。
図６に示す横長の矩形は、実験室の内壁面を示す。実験室の大きさは、縦３．５ｍ、横６．５ｍ、高さ２．７ｍの直方体である。雑音源３１は、実験室のほぼ中央部に配置されている。信号入力部１１の重心点は、雑音源３１から実験室の左端に１．０ｍに離れた位置に配置されている。信号入力部１１は、８個のマイクロホンを備えるマイクロホンアレイである。図６において、方向φは、信号入力部１１の重心点から雑音源への方向とは逆方向を基準とする方位角で表される。ここで、雑音源の方向は１８０°である。音源３２は、信号入力部１１の重心点から雑音源とは異なる方向φに１．０ｍ離れた位置に配置されている。
【００６４】
次に実験で用いた信号入力部１１の構成について説明する。
図７は、信号入力部１１の構成例を示す概略図である。
信号入力部１１は、水平面上に重心点を中心とする直径０．３ｍの円周上に等間隔（４５°）で８個の無指向性マイクロホン１１１−１〜１１１−Ｍが配置されている。
【００６５】
次に実験で用いた雑音の例について説明する。
図８は、実験で用いた雑音のスペクトルの一例を示す図である。
図８の横軸は周波数、縦軸はパワーを表す。実験で用いた雑音は、約２５０Ｈｚにパワーのピークを有し、このピークに係る周波数よりも高い周波数では、周波数が高くなるに従いパワーが単調に低くなる。この雑音は、約６００Ｈｚよりも低い周波数の低域成分が主である。
【００６６】
次に実験で用いた目的音の例について説明する。
図９は、実験で用いた目的音のスペクトルの一例を示す図である。
図９の横軸は周波数、縦軸はパワーを表す。実験で用いた目的音は、約３５０Ｈｚにパワーのピークを有する。このピークに係る周波数より高い周波数では、概ね周波数が高くなるに従いパワーが低くなる傾向があるが、単調にパワーが低くなるとは限らない。実験で用いた目的音は、約１３００Ｈｚ、３０００Ｈｚにおいて、パワーの概形において、それぞれ滑らかなボトム（極小）、ピーク（極大）を有する。なお、実験で用いた目的音として音楽が用いられたため、スペクトルは時刻の経過に伴って変動する。
【００６７】
その他、実験における条件は次の通りである。周波数領域変換部１２１、時間領域変換部１２７におけるＦＦＴ点数は１０２４である。ＦＦＴ点数とは、１フレームに含まれる信号のサンプル数である。シフト長、即ち、各フレームの先頭サンプルに係る隣接フレーム間のサンプル位置のずれは５１２である。周波数領域変換部１２１では、フレーム毎に抽出した音響信号に窓関数としてブラックマン窓（Ｂｌａｃｋｍａｎ  ｗｉｎｄｏｗ）をかけて生成した時間領域の信号を周波数領域係数に変換した。
【００６８】
（位相差の変化例）  
次に、遅延和要素行列算出部１２４があるフレームｋにおいて算出した位相差θｍ，ｋ（ｔ）の一例について説明する。以下の説明では、位相差θｍ，ｋ（ｔ）においてフレーム、繰り返しを示すインデックスｋ、ｔをそれぞれ省略して、チャネルｍについて算出した位相差をθｍ（ｍは、１より大きく８と等しいか８より小さい整数である）と表す。また基準とするチャネル１からの位相差をθ1として示すが、θ1はその定義より常に０であり、チャネル１の位相を任意に取ることができることから、これを０と定めれば、θｍを単に位相と呼んでも差し支えない。
【００６９】
図１０は、繰り返しｔによる位相差θｍの変化の一例を示す図である。
図１０において、縦軸は位相差（ラジアン）を示し、横軸は繰り返し（回数）を示す。
本実施形態では、各チャネルに係る位相差θ２，…，θ８の初期値（即ち、ｔ＝０）は、上述したようにランダムな値であるが繰り返しが増加すると単調に一定値に収束する。繰り返しｔが９０回を超えると、位相差θ２，…，θ８はそれぞれ一定値に達する。
【００７０】
（特異値の例）
次に、特異ベクトル算出部１２５が算出した特異値ｍの一例について説明する。
図１１は、特異値σｍの区間数Ｑによる依存性の一例を示す図である。
図１１において、縦軸は特異値を示し、横軸は区間数Ｑを示す。但し、図１１に示す特異値σ１，…，σ８は、上述のように位相差θｍの初期値としてランダムな値が設定され、位相差θｍが十分に収束した後に得られた遅延和要素行列Ｃに基づいて算出されたものである。
図１１に示すように、特異値σ１，…，σ８は、各次数ともに区間数Ｑが増加するほど増加する。区間数Ｑが８よりも小さい場合、ゼロ又はゼロに近似する特異値が、少なくとも１個ある。即ち、その少なくとも１個の特異値にそれぞれ対応する右特異ベクトルは、雑音を抑圧する効果をもたらさない。他方、区間数Ｑが２０よりも大きい場合、全ての特異値σ１，…，σ８は、１よりも大きくなる。即ち、各特異値にそれぞれ対応する右特異ベクトルは、雑音を抑圧する効果をもたらす。本実験では、雑音が１方向のみから到来するため、７個の特異値がゼロとは有意に異なる（非零）の特異値であり、１個の特異値はゼロ又はゼロに近似する特異値であるべきとも考えられる。しかし、８個の非零の特異値が得られたのは、実験室の内壁や設置物による反射のためであると考えられる。
【００７１】
次に、特異ベクトル算出部１２５が算出した特異値ｍの他の例について説明する。
図１２は、特異値σｍの区間数Ｑによる依存性の他の例を示す図である。
図１２における縦軸と横軸が示す関係は、図１１と同様である。この例では、位相差θｍの初期値としていずれもゼロが設定され、位相差θｍが十分に収束した後に得られた遅延和要素行列Ｃに基づいて算出されたものである。
図１２に示す特異値σ１，…，σ８は、各次数ともに区間数が増加するほど増加する。但し、特異値σ１が、他の特異値σ２，…，σ８よりも顕著に大きい値をとる。区間数が８０の場合でも、１を超える特異値は、特異値σ１の他にσ２，σ３の２個に過ぎない。区間数がより多いほど、１を超える特異値がより多く算出される可能性はあるが、処理量が過大になる。即ち、図１１、１２は、本実施形態のように位相差θｍの初期値としてランダムな値を設定することによって、効率的に特異ベクトルを算出でき、その算出した特異ベクトルを用いて十分な雑音抑圧性能を得ることができることを裏付けている。
【００７２】
（出力音響信号の例）  
次に、チャネルｍについて時間領域変換部１２７が算出した出力音響信号の例について説明する。
図１３は、出力音響信号のスペクトログラムの一例を示す図である。
図１３において、（ａ）は位相差θｍの初期値としていずれもゼロが設定された場合、（ｂ）は位相差θｍの初期値としてランダムな値が設定された場合、（ｃ）は位相差θｍの初期値として（ｂ）とは異なるランダムな値が設定された場合、（ｄ）は位相差θｍの初期値として（ｂ）、（ｃ）ともに異なるランダムな値が設定された場合を示す。（ａ）〜（ｄ）ともに、縦軸は周波数（Ｈｚ）、横軸は時刻（ｓ）を示し、出力音響信号のレベルを濃淡で示す。暗い領域ほどレベルが低く、明るい領域ほどレベルが高いことを示す。
また、（ａ）〜（ｄ）ともに、間欠的に広い周波数帯域にわたってレベルが高くなる時間帯があることを示す。この時間帯は、目的音が到来している時間帯であり、それ以外の時間帯は雑音のみが到来している時間帯であることを示す。図１３は、（ａ）において雑音のレベルが高い領域が最も広いことを示す。即ち、（ｂ）〜（ｄ）は、位相差θｍの初期値としてランダムな値が設定されることで雑音が効果的に抑圧されることを示す。
【００７３】
次に、ある区間においてチャネルｍについて時間領域変換部１２７が算出した出力音響信号の他の例について説明する。
図１４は、出力音響信号のスペクトログラムの他の例を示す図である。
但し、図１４に示す出力音響信号は、入力信号ベクトル［ｙｋ］と右特異ベクトル［v１］，…，［ｖ８］のうち各１個のみに基づいて得られた出力周波数領域係数ｚｋを時間領域に変換した信号である。これらを出力音響信号１〜８と呼ぶ。右特異ベクトル［v１］，…，［ｖ８］は、位相差θｍの初期値としてランダムな値を設定して算出した遅延和要素行列Ｃに基づく。
出力音響信号１〜８のスペクトログラムを、図１４（ａ）〜（ｈ）にそれぞれ示す。
図１４（ａ）〜（ｈ）のそれぞれについて、縦軸、横軸及び濃淡の関係は、図１３（ａ）〜（ｄ）と同様である。周囲よりも雑音のレベルが高い領域に注目すると、（ａ）〜（ｈ）間で雑音のレベルが高い領域の広さは概ね同等である反面、図１４（ｈ）に示される雑音のレベルが最も高い。つまり、図１４は、出力音響信号８に雑音の成分が集中し、出力音響信号１〜７では雑音が抑圧されていることを示す。
【００７４】
図１５は、出力音響信号のスペクトログラムのさらに他の例を示す図である。
図１５に示す出力音響信号も、出力音響信号１〜８と同様に、入力信号ベクトル［ｙｋ］と右特異ベクトル［v１］，…，［ｖ８］のうち各１個のみに基づいて得られた出力周波数領域係数ｚｋを時間領域に変換した信号である。これらを出力音響信号１’〜８’と呼ぶ。但し、右特異ベクトル［v１］，…，［ｖ８］は、位相差θｍの初期値としていずれもゼロを設定して算出した遅延和要素行列Ｃに基づく。
出力音響信号１’〜８’のスペクトログラムを、図１５（ａ）〜（ｈ）にそれぞれ示す。図１５（ａ）〜（ｈ）のそれぞれについて、縦軸、横軸及び濃淡の関係は、図１５（ａ）〜（ｈ）と同様である。これによれば、周囲よりも雑音のレベルが高い領域の広さ、その領域における雑音のレベルは、図１５（ａ）〜（ｈ）間でまちまちである。従って、位相差θｍの初期値としてゼロが設定されると、正しく遅延和要素行列Ｃを算出できないために、必ずしも効果的に雑音が抑圧されないことを示す。
【００７５】
（平均ＭＵＳＩＣスペクトルの例）  
次に、方向算出部２２１３が算出する平均ＭＵＳＩＣスペクトルＰａｖｇ（φ）の例について説明する。
図１６は、平均ＭＵＳＩＣスペクトルＰａｖｇ（φ）の一例を示す図である。
図１６の横軸は方位角（°）を示し、縦軸は平均ＭＵＳＩＣスペクトルＰａｖｇ（φ）のパワー（ｄＢ）を示す。
図１６は、方位角１８０°において平均ＭＵＳＩＣスペクトルＰａｖｇ（φ）のパワーが最大であるピークを示す。方向算出部２２１３は、このパワーが最大となるピークを与える方位角１８０°を音源の方向と定める。
【００７６】
（音源方向の例）  
次に、方向算出部２２１３が定めた音源の方向φの例について説明する。
図１７は、本実施形態に係る方向算出部２２１３が定めた音源の方向φの一例を示す図である。
上述したように、ＭＵＳＩＣスペクトルＰ（φ）を算出する際に用いる共役転置行列［Ｖｃ］Ｈを、Ｍ’個の右特異ベクトル［v１］，…，［ｖＭ’］を統合して生成する。
図１７（ａ）〜（ｆ）では、共役転置行列［Ｖｃ］Ｈに含まれている右特異ベクトルの数Ｍ’が８〜３個、それぞれの場合について方向φを示す。実験では、信号入力部１１からの方向０°、４５°、９０°、１３５°、１８０°、２２５°、２７０°、３１５°にそれぞれ異なる時刻に音源を設置し、音を発生させた。
【００７７】
図１７（ａ）〜（ｆ）において、横軸は時刻（ｓ）、縦軸は方位角（°）を示す。×印は、目的音を放射する音源の方向を示す。
図１７（ａ）は、右特異ベクトルの数Ｍ’が８個の場合、最も高精度に音源の方向φを推定できることを示す。
図１７（ｂ）〜（ｅ）は、右特異ベクトルの数Ｍ’が７〜４個の場合、ほぼ音源の方向φを現実の音源の方向に推定できることを示す。但し、現実に音源が存在しないにも関わらず、音源の方向φを約３３０°と推定することがある。
図１７（ｆ）は、右特異ベクトルの数Ｍ’が３個に減少すると、ほとんど音源の方向φを推定することができないことを示す。これは、出力音響信号のチャネル数が減少することで、特定の方向から到来する雑音を抑圧するベクトル空間を十分に利用できないことによる。
【００７８】
次に、上述した実験と同様な条件で従来のＭＵＳＩＣ法を用いて推定した音源の方向φの例について説明する。
図１８は、従来のＭＵＳＩＣ法を用いて推定した音源の方向φの一例を示す図である。
図１８における縦軸、横軸の関係は、図１７と同様である。
図１８は、方位角１８０°に設置された雑音源の方向が音源の方向φとして常に推定されることを示す。つまり、本実施形態のように雑音が抑圧されないことを示す。また、図１８は、音源の方向が１３５°、１８０°、２２５°の場合には、雑音源と区別できないことを示す。これは、雑音源が放射する雑音のスペクトルの周波数帯域と音源が放射する目的音のスペクトルの周波数帯域が互いに重なるため、両者を区別できないことによる。
裏返せば、本実施形態では従来のＭＵＳＩＣ法とは異なり、雑音源と同一又は近似する方向にある音源から到来する目的音の成分を抽出し、その方向を推定することができるという従来のＭＵＳＩＣ法では得られなかった効果を奏する。
【００７９】
以上に、説明したように本実施形態では、第１の実施形態の構成を備え、第１の実施形態で算出した出力信号に基づいて算出した相関行列を対角化して固有ベクトルを算出する。本実施形態では、算出した固有ベクトル、第１の実施形態で算出した特異ベクトル及び方向毎の伝達特性を示す伝達関数に基づいて方向毎のスペクトルを算出し、算出したスペクトルが極大となる方向を定める。
そのため、本実施形態では、第１の実施形態と同様な効果を奏することで、雑音を抑圧して目的音が残るため、残った目的音の方向を高精度に推定することができる。
【００８０】
なお、上述した実施形態における音響信号処理装置１２、２２の一部、例えば、周波数領域変換部１２１、入力信号行列生成部１２２、初期値設定部１２３、遅延和要素行列算出部１２４、特異ベクトル算出部１２５、出力信号ベクトル算出部１２６、時間領域変換部１２７、及び方向推定部２２１をコンピュータで実現するようにしても良い。その場合、この制御機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、音響信号処理装置１２、２２に内蔵されたコンピュータシステムであって、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
また、上述した実施形態における音響信号処理装置１２、２２の一部、または全部を、ＬＳＩ（Ｌａｒｇｅ  Ｓｃａｌｅ  Ｉｎｔｅｇｒａｔｉｏｎ）等の集積回路として実現しても良い。音響信号処理装置１２、２２の各機能ブロックは個別にプロセッサ化してもよいし、一部、または全部を集積してプロセッサ化しても良い。また、集積回路化の手法はＬＳＩに限らず専用回路、または汎用プロセッサで実現しても良い。また、半導体技術の進歩によりＬＳＩに代替する集積回路化の技術が出現した場合、当該技術による集積回路を用いても良い。
【００８１】
以上、図面を参照してこの発明の一実施形態について詳しく説明してきたが、具体的な構成は上述のものに限られることはなく、この発明の要旨を逸脱しない範囲内において様々な設計変更等をすることが可能である。
【符号の説明】
【００８２】
１、２…音響信号処理システム
１１…信号入力部、
１１１−１〜１１１−Ｍ…マイクロホン、
１２、２２…音響信号処理装置、
１２１…周波数領域変換部、１２２…入力信号行列生成部、１２３…初期値設定部、
１２４…遅延和要素行列算出部、１２５…特異ベクトル算出部、
１２６…出力信号ベクトル算出部、１２７…時間領域変換部、
１３…信号出力部、
２２１…方向推定部、
２２１１…相関行列算出部、２２１２…固有ベクトル算出部、２２１３…方向算出部、
２３…方向出力部、
【図２】
2014026115.tif 000010
【図３】
2014026115.tif 000011
【図４】
2014026115.tif 000012
【図５】
2014026115.tif 000013
【図６】
2014026115.tif 000014
【図１０】
2014026115.tif 000015
【図１１】
2014026115.tif 000016
【図１２】
2014026115.tif 000017
【図１】
2014026115.tif 000018
【図７】
2014026115.tif 000019
【図８】
2014026115.tif 000020
【図９】
2014026115.tif 000021
【図１３】
2014026115.tif 000022
【図１４】
2014026115.tif 000023
【図１５】
2014026115.tif 000024
【図１６】
2014026115.tif 000025
【図１７】
2014026115.tif 000026
【図１８】
2014026115.tif 000027