(57)【要約】【課題】原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから、原音声信号を精度よく分離して、音声信号の明瞭度を向上させることができるようにする。【解決手段】パラメータ更新部２４により、パワースペクトル密度系列モデルΦxk[t]と、観測時間周波数成分ｙk[t]との距離が小さくなるように、原音声信号のパワースペクトル密度Φ(s)k[t]の係数行列の各要素Ｇ(s)m[t]と、雑音信号のパワースペクトル密度Φ(n)k[t]の基底行列の各要素Ｂ(n)q,k及び係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]を更新することを繰り返す。【選択図】図２
2015049406.tif 000002
【特許請求の範囲】
【請求項１】
音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分ｙk[t]を出力する時間周波数解析部と、
前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、Ｍ個の基底スペクトルｍを表す非負値の各要素Ｂ(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(s)m,kのゲインＧ(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素Ｇ(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Ｑ個の基底スペクトルqを表す非負値の各要素Ｂ(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素Ｂ(n)q,k、前記係数行列の各要素Ｇ(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素Ｂ(s)m,kに、予め求められた値を設定するパラメータ初期値設定部と、
（τ、m）の全ての組み合わせにおける、前記基底行列の要素Ｂ(s)m,k、前記係数行列の要素Ｇ(s)m[t]、及び前記パワーＨk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]とを更新するパラメータ更新部と、
予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、
を含む音響信号解析装置。
【請求項２】
前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離の尺度としてβ−ｄｉｖｅｒｇｅｎｃｅを用いた
請求項１記載の音響信号解析装置。
【請求項３】
補助変数更新部を更に含み、
前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離を表す関数を、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]と、（k、m、t、τ）の全ての組み合わせについての補助変数ζk,m,t,τと、（k、q、t）の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]を用いて表され、かつ、前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]とのβ−ｄｉｖｅｒｇｅｎｃｅの上限関数である補助関数とし、
前記補助変数更新部は、前記補助関数を小さくするように、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、（k、m、t、τ）の全ての組み合わせについての補助変数ζk,m,t,τと、（k、q、t）の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]を更新し、
前記パラメータ更新部は、前記補助関数を小さくするように、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、（k、m、t、τ）の全ての組み合わせについての補助変数ζk,m,t,τと、（k、q、t）の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]とに基づいて、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]とを更新する
請求項２記載の音響信号解析装置。
【請求項４】
時間周波数解析部によって、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分ｙk[t]を出力し、
パラメータ初期値設定部によって、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、Ｍ個の基底スペクトルｍを表す非負値の各要素Ｂ(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(s)m,kのゲインＧ(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素Ｇ(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Ｑ個の基底スペクトルqを表す非負値の各要素Ｂ(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素Ｂ(n)q,k、前記係数行列の各要素Ｇ(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素Ｂ(s)m,kに、予め求められた値を設定し、
パラメータ更新部によって、（τ、m）の全ての組み合わせにおける、前記基底行列の要素Ｂ(s)m,k、前記係数行列の要素Ｇ(s)m[t]、及び前記パワーＨk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]とを更新し、
終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う
音響信号解析方法。
【請求項５】
コンピュータを、請求項１〜請求項３の何れか１項に記載の音響信号解析装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、音響信号の時系列データから、原音声信号を強調するための音響信号解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音声強調技術は、例えば、ハンズフリー通話システムにおいて通信音声の了解性を向上させたり、自動会議録作成システムにおいて音声認識率を向上させたりするための重要な要素技術である。
【０００３】
雑音や残響環境の情報は事前に知り得ないことが多いため、観測信号のみから音声信号を獲得する問題を扱う必要があるが、音声信号に残響と雑音が重畳されて観測信号が得られるプロセスを順問題と捉えると、観測信号から音声信号を得る問題は逆問題と見なすことができる。
【０００４】
雑音や室内伝達系の情報が未知の場合、この逆問題には解が無数に存在しうるため、解を絞り込むための何らかの仮定が必要である。音声強調に関する従来アプローチでは、雑音に定常性、室内伝達系に時不変性などの制約を仮定し、その仮定の下で当該逆問題が定式化されるものが多かった。
【０００５】
例えば、古典的な雑音除去手法であるスペクトラルサブトラクション法では、雑音の定常性が仮定され、既知の非音声区間から得られる雑音のパワー（または振幅）スペクトルを観測信号のパワー（または振幅）スペクトルから減算することで音声信号のパワー（または振幅）スペクトルを推定する。また、逆フィルタリングに基づく残響除去法では、時不変かつ最小位相な室内伝達系が仮定され、音声に関する仮定やモデル（調波性、自己回帰モデル、自己相関関数コードブック等）を手がかりに最適な残響除去フィルタを推定する。
【０００６】
こうした従来手法は、仮定が成立する状況においては高い性能を発揮する一方で、実環境においては、雑音が例えば携帯電話の着信音及び背景音楽などのように非定常なものである場合や、僅かでも音源移動がある場合等では、雑音や室内伝達系に対する上述の仮定が成立しない場合があり、このような環境の下では高い性能を発揮できなくなるという問題があった。
【０００７】
こうした室内伝達系の変化に対する頑健な残響除去の手法として、音声スペクトログラムと室内インパルス応答のスペクトログラムの時間方向の畳み込みモデルに基づき、音声スペクトログラムのスパース性を手がかりに音声スペクトログラムを推定する方法が提案されている（非特許文献１）。
【０００８】
しかし、この手法では残響のような乗法的な雑音の混入プロセスしか想定されておらず、加法的な雑音の混入プロセスは陽に想定されていない（雑音が存在しない状況を仮定している）ため、観測信号に雑音が含まれる場合には高い性能を発揮できないという問題があった。
【０００９】
一方、非定常な雑音を抑圧するための手法として、スペクトログラムを非負値行列と見立て、これを二つの非負値行列に分解する非負値行列因子分解アプローチが近年注目されている（非特許文献２）。非負値行列因子分解を用いれば、観測スペクトログラムを構成する基底スペクトルと各時刻におけるそれらの結合係数を同時に得ることができる。これにより、例えばクリーンな音声サンプルからあらかじめ基底スペクトルを獲得（学習）しておき、雑音つき音声の観測スペクトログラムを基底行列と係数行列に分解する際、基底スペクトルセットの一部をあらかじめ学習した基底スペクトルのセットに固定しておくことで、音声とそれ以外（すなわち雑音）のスペクトログラムに分解する。
【先行技術文献】
【非特許文献】
【００１０】
【非特許文献１】H. Kameoka, T. Nakatani, T. Yoshioka, “Robust Speech Dereverberation Based on Non-negativity and Sparse Nature of Speech Spectrograms,” in Proc. 2008 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP2009), pp. 45−48, 2009.
【非特許文献２】P. Smaragdis, B. Raj and M. V. Shashanka, “Supervised and semi-supervised separation of sounds from single-channel mixtures,” in Proc. ICA2007, pp. 414−421, 2007.
【発明の概要】
【発明が解決しようとする課題】
【００１１】
しかし、上記の非特許文献２の手法では加法的な雑音の混入プロセスしか想定されておらず、残響のような乗法的な雑音の混入プロセスは陽に想定されていない（残響がない状況を仮定している）ため、観測信号に残響が含まれる場合には高い性能を発揮できないという問題があった。
【００１２】
本発明は、上記の事情を鑑みてなされたもので、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１３】
上記の目的を達成するために本発明に係る音響信号解析装置は、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分ｙk[t]を出力する時間周波数解析部と、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、Ｍ個の基底スペクトルｍを表す非負値の各要素Ｂ(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(s)m,kのゲインＧ(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素Ｇ(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Ｑ個の基底スペクトルqを表す非負値の各要素Ｂ(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素Ｂ(n)q,k、前記係数行列の各要素Ｇ(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素Ｂ(s)m,kに、予め求められた値を設定するパラメータ初期値設定部と、（τ、m）の全ての組み合わせにおける、前記基底行列の要素Ｂ(s)m,k、前記係数行列の要素Ｇ(s)m[t]、及び前記パワーＨk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]とを更新するパラメータ更新部と、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、を含んで構成されている。
【００１４】
本発明に係る音響信号解析方法は、時間周波数解析部によって、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分ｙk[t]を出力し、パラメータ初期値設定部によって、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、Ｍ個の基底スペクトルｍを表す非負値の各要素Ｂ(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(s)m,kのゲインＧ(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素Ｇ(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Ｑ個の基底スペクトルqを表す非負値の各要素Ｂ(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素Ｂ(n)q,k、前記係数行列の各要素Ｇ(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素Ｂ(s)m,kに、予め求められた値を設定し、パラメータ更新部によって、（τ、m）の全ての組み合わせにおける、前記基底行列の要素Ｂ(s)m,k、前記係数行列の要素Ｇ(s)m[t]、及び前記パワーＨk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分ｙk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分ｙk[t]と、前記基底行列の各要素Ｂ(s)m,kと、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素Ｇ(s)m[t]と、前記基底行列の各要素Ｂ(n)q,kと、前記係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]とを更新し、終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う。
【００１５】
本発明に係るプログラムは、上記の音響信号解析装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１６】
以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、原音声信号sのパワースペクトル密度Φ(s)k[t]の基底行列及び係数行列と、室内インパルス応答のパワーＨk[τ]と、雑音信号nのパワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、観測時間周波数成分ｙk[t]との距離が小さくなるように、パワースペクトル密度Φ(s)k[t]の係数行列と、パワースペクトル密度Φ(n)k[t]の基底行列及び係数行列と、室内インパルス応答のパワーＨk[τ]を更新することを繰り返すことにより、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる、という効果が得られる。
【図面の簡単な説明】
【００１７】
【図１】マイクロホンの位置を変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分及び位相成分を示す図である。
【図２】本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。
【図３】本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。
【図４】（Ａ）メル周波数ケプストラム係数歪みの評価結果を示すグラフ、（Ｂ）バークスペクトル歪みスコアの評価結果を示すグラフ、及び（Ｃ）信号対干渉比の評価結果を示すグラフである。
【発明を実施するための形態】
【００１８】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１９】
＜発明の原理＞
【００２０】
まず、本発明の原理について説明する。
【００２１】
＜音声強調問題の定式化＞
【００２２】
室内伝達系が時不変の場合、原音声信号をｓ[t]、室内インパルス応答をｈ[t]、雑音信号をｎ[t] とすると、音響信号としての観測信号ｙ[t]は時間領域において、以下の（１）式のように表される。
【００２３】
【数１】
2015049406.tif 000003
【００２４】
ただし、tは離散時刻のインデックスである。右辺第1項Στｓ[t-τ]ｈ[τ]は残響音声信号を表す。ここで、観測信号、原信号、室内インパルス応答、雑音信号の短時間フーリエ変換をそれぞれｙk[t]、ｓk[t]、ｈk[t]、ｎk[t]とする。ただし、k=1,...,Kは周波数のインデックス、t=1,...,Tはフレーム（時刻）のインデックスを表す。時間領域における信号同士の畳み込みは、特定の条件の下でそれぞれの信号の各周波数における狭帯域信号同士の畳み込みで近似できることが知られている。もしこの近似が成立するような短時間フーリエ変換の分析条件を選べば、ｙk[t]は（２）式のように表される。
【００２５】
【数２】
2015049406.tif 000004
【００２６】
以上では時不変な系を考えたが、実際には音源が移動することがあるため、室内伝達系の時不変性の仮定は必ずしも成り立たない。ここで、時変な室内伝達系を仮定することは、（２）式において室内インパルス応答ｈk[τ]が時刻tごとに変化しうることを許容することに相当する。すなわち、以下の（３）式のように表される。
【００２７】
【数３】
2015049406.tif 000005
【００２８】
非定常雑音、時変残響環境下における音声強調問題は、（３）式をもとに観測時間周波数成分Ｙ={ｙk[t]}1≦k≦K,1≦t≦TからＳ={ｓk[t]}1≦k≦K,1≦t≦Tを推定する問題として定式化される。（３）式では室内インパルス応答が時刻ごとに変化することを許容しているため、（２）式に比べて未知パラメータが圧倒的に多くなっており、同じ右辺を与えるｓとｈとｎの組み合わせは無数に存在する。そこで以下では、その無数の組み合わせの中から正しい解を得るためにどのような仮定を置くべきか、について説明する。
【００２９】
＜室内伝達系の「半時変性」＞
【００３０】
室内インパルス応答が時間的に変化する主要因の１つは話者の移動である。なぜなら、話者の移動に伴い、話者からマイクロホンまでの（壁からの反射音の到来経路を含む）あらゆる音響経路の距離が一斉に変化するからである。
【００３１】
ここで、単一の音響経路のみに着目すると、伝達系のインパルス応答はデルタ関数となるが、経路の距離に変化があった場合、当該インパルス応答のフーリエ変換の位相成分には音響経路の距離変化に伴う到来時間の変化、振幅成分には音響経路の距離変化に伴う強度の変化がそれぞれ反映されることになる。音の強度は伝播距離に対して反比例的であるため、伝播距離がある程度大きい場合には伝播距離に変化があっても振幅成分はほとんど変化しないことになる。
【００３２】
一方、位相成分については、到来時間と周波数に比例した量だけ変化するため、到来時間の変化がたとえわずかであっても、特に高周波域においては著しく変化しうる。
【００３３】
よって、単一の音響経路の伝達系のインパルス応答の短時間フーリエ変換には、話者の移動に伴って振幅は変化しにくく位相は変化しやすい、という性質があることが予想される。
【００３４】
実際の室内インパルス応答はあらゆる音響経路の伝達系のインパルス応答の重ね合わせとなるので、一つ一つのインパルス応答が上述のような傾向があるならば、その重ね合わせもおおよそ同様な傾向があるのではないかと考えられる。
【００３５】
図１は、話者の代わりにマイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分と位相成分を比較したものである。
【００３６】
上記図１の上段は、室内インパルス応答の短時間フーリエ変換の振幅成分（Ａ）及び位相成分（Ｂ）を示したものであり、上記図１の下段は、マイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分(Ｃ)と位相成分(Ｄ)を示したものである。なお、室内インパルス応答のデータはRWCP実環境音声・音響データベースのマイクロホンアレーデータベースに収録されている、残響時間0.3[s]の残響可変室（パネル）で測定されたものである。
【００３７】
上記図１から、振幅成分はさほど変化していないのに対し位相成分は著しく変化していることが分かる。
【００３８】
以上より、発明者らは、話者の移動等、環境の軽微な変化については、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変と仮定した特殊な系により扱えるのではないかと考えた。以後、このような系を「半時変系」と呼ぶこととする。
【００３９】
＜観測信号の生成モデル＞
【００４０】
以上の２つの仮定をベースに、観測信号の生成モデルを定式化する。（３）式に半時変性の仮定を導入すると、観測信号の短時間フーリエ変換は（４）式のように表される。
【００４１】
【数４】
2015049406.tif 000006
【００４２】
ただし、φk,t[τ]とＨk[t]は室内インパルス応答の短時間フーリエ変換の位相成分と振幅成分であり、Ｈk[t]=|ｈk,t[τ]|である。ここで、φk,t[τ]をk,n,mごとに独立に区間[0;2π)上の一様分布に従う確率変数とする。さらに、原音声信号、雑音については、それぞれ平均が０の複素正規分布（５）式及び（６）式に従う確率変数とする。
【００４３】
【数５】
2015049406.tif 000007
【００４４】
ただし、Φ(s)k[t]、Φ(n)k[t]はそれぞれ原音声信号および雑音信号の時刻t、k番目の周波数ビンにおけるパワースペクトル密度である。ここで、ＮC(z;0,λ)は複素正規分布の確率密度関数を表し、（７）式によって与えられる。
【００４５】
【数６】
2015049406.tif 000008
【００４６】
証明は省略するが、原音声信号と雑音のガウス性の仮定と、正規分布の再生性より、観測信号の短時間フーリエ変換ｙk[t]はk、tごとに独立に、（８）式に示す複素正規分布に従うことが示される。
【００４７】
【数７】
2015049406.tif 000009
【００４８】
次に、原音声信号のパワースペクトル密度系列Φ(s)k[t]に対し、非負値行列積表現を導入する。基底スペクトルをＭ個とし、m番目の基底スペクトルをＢ(s)m,kとする。また、時刻tにおける基底スペクトルＢ(s)m,kのゲインをＧ(s)m[t]とすると、原音声信号のパワースペクトル密度系列は（９）式及び（１０）式で示される。
【００４９】
【数８】
2015049406.tif 000010
【００５０】
なお、基底スペクトルとそのゲインは非負であることに注意する．また、スケールの任意性を除くため、（１１）式に示す条件を満たすものとする。
【００５１】
【数９】
2015049406.tif 000011
【００５２】
音声の基底スペクトルＢ(s)m,kは未知変数として観測信号から他の未知パラメータとしてブラインドで推定することもありえるが、ある程度の量のクリーン音声サンプルからあらかじめ学習しておくことを想定する。その意味で、本実施の形態に係る手法（以下、提案法と称する）はセミブラインドな音声強調手法である。基底の事前学習は、クリーン音声データのスペクトログラムに対してＮＭＦを適用することで行うことができる。
【００５３】
＜パラメータ推定アルゴリズム＞
【００５４】
以上で立てた観測時間周波数成分Ｙ={ｙk[t]}1≦k≦K,1≦t≦Tの確率密度関数は、未知パラメータの尤度関数に対応する。（８）式より、観測時間周波数成分Ｙが与えられたもとでの未知パラメータの対数尤度は、定数項・符号を無視すれば、観測スペクトログラム|ｙk[t]|2と（１２）式で示されるパワースペクトル密度系列モデルとの距離を表す板倉斎藤距離（１３）式と等しくなる。
【００５５】
【数１０】
2015049406.tif 000012
【００５６】
従って、提案法は観測スペクトログラム|ｙk[t]|2とパワースペクトル密度系列モデルΦxk[t]との最適フィッティング問題に帰着する。ここで、雑音信号のパワースペクトル密度系列Φ(n)k[t]に関しても、（１４）式のようにモデル化される。
【００５７】
【数１１】
2015049406.tif 000013
【００５８】
なお、板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))は、距離尺度を表すβ-divergenceにおいてβを０にした場合に得られる距離に相当する。
【００５９】
（１３）式で示される板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))は直接最小化することが困難であるため、板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))を下回らず１点で接する上限関数を設定する。導出は省略するが、（１５）式に示す上限関数が板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))の上限関数となる。
【００６０】
【数１２】
2015049406.tif 000014
【００６１】
すなわち、板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))の代わりに（１５）式に示された上限関数を小さくするようにＨ,Ｇ(s),Ｂ(n),Ｇ(n)および補助変数ζ,ξ,ηを交互に更新することで目的関数Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))を小さくしていくことができる。この際、Ｈ,Ｇ(s),Ｂ(n),Ｇ(n)の各パラメータは非負値を保ったまま更新するようにする。なお、（１５）式におけるΦ(x)k[t]は（１６）式によって示される。
【００６２】
【数１３】
2015049406.tif 000015
【００６３】
そして、この上限関数Ｄ0の更新則は（１７）式〜（２３）式で示される。
【００６４】
【数１４】
2015049406.tif 000016
＜システム構成＞
【００６５】
次に、原音声信号に雑音及び残響成分が重畳された観測信号を解析して、観測信号に含まれる原音声信号を推定する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００６６】
図２に示すように、本発明の実施の形態に係る音響信号解析装置は、ＣＰＵと、ＲＡＭと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【００６７】
音響信号解析装置１００は、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。
【００６８】
入力部１０により、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが入力される。記憶部３０は、入力部１０により入力された観測信号の時系列データを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、本処理ルーチンで用いる各パラメータの初期値を記憶している。
【００６９】
演算部２０は、時間周波数解析部２１と、初期設定部２２と、補助変数更新部２３と、パラメータ更新部２４と、終了判定部２５と、信号変換部２６とを備えている。
【００７０】
時間周波数解析部２１は、例えばマイクロホンの時系列信号としての観測された観測信号ｙ[t]を入力として、観測信号ｙ[t]の観測時間周波数成分Ｙ={ｙk[t]}（周波数k＝１，・・・，K、時刻t＝１，・・・，T）を計算する。また、計算した観測時間周波数成分Ｙを、記憶部３０に記憶しておく。より詳細には、時間周波数解析部２１は、例えばマイクロホンで観測された観測信号の時系列データを入力として、短時間フーリエ変換（Short-Time Fourier Transform；ＳＴＦＴ）を用いて時間周波数解析を行うことにより観測時間周波数成分Ｙを計算する。
【００７１】
初期設定部２２は、後述する処理で用いる各パラメータＧ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]の各初期値及びＢ(s)m,k の値を設定する。なお、Ｂ(s)m,k以外のパラメータの初期値は、例えば乱数を用いて適当な値に設定すればよい。この場合、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]の各パラメータの初期値は非負値となるように設定する。
【００７２】
そして、Ｂ(s)m,kの値には、一定量のクリーン音声サンプルから予め学習して得た値を設定する。なお、スケールの任意性を排除するため、Ｂ(s)m,kの値は上記（１１）式を満たすように予め学習されているものとする。
【００７３】
補助変数更新部２３は、（k、m、t、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているＢ(s)m,k、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]に基づいて、上記（２１）式に従って上限関数Ｄ0を小さくするように補助変数ζk,m,t,τを更新し、記憶部３０に格納する。
【００７４】
また、補助変数更新部２３は、（k、q、t）の全ての組み合わせの各々について、記憶部３０に記憶されているＢ(s)m,k、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]に基づいて、上記（２２）式に従って上限関数Ｄ0を小さくするように補助変数ξk,q,tを更新し、記憶部３０に格納する。
【００７５】
更に、補助変数更新部２３は、各時刻t及び各周波数kについて、記憶部３０に記憶されているＢ(s)m,k、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]に基づいて、上記（１６）式、上記（２３）式に従って上限関数Ｄ0を小さくするように補助変数ηk[t]を更新し、記憶部３０に格納する。
【００７６】
パラメータ更新部２４は、（m、t）の全ての組み合わせの各々について、記憶部３０に記憶されているｙk[t]、Ｂ(s)m,k、Ｈk[τ]、ζk,m,t,τ、ηk[t]に基づいて、上記（１７）式に従って、上限関数Ｄ0を小さくするように、基底スペクトルＢ(s)m,kのゲイン係数Ｇ(s)m[t]を更新し、記憶部３０に格納する。
【００７７】
また、パラメータ更新部２４は、（q、k）の全ての組み合わせの各々について、記憶部３０に記憶されているｙk[t]、Ｇ(n)q[t]、ξk,q,t、ηk[t]に基づいて、上記（１８）式に従って、上限関数Ｄ0を小さくするように、雑音信号nのパワースペクトル密度Φ(n)k[t]の基底行列の要素Ｂ(n)q,kを更新し、記憶部３０に格納する。
【００７８】
また、パラメータ更新部２４は、（q、t）の全ての組み合わせの各々について、記憶部３０に記憶されているｙk[t]、Ｂ(n)q,k、ξk,q,t、ηk[t]に基づいて、上記（１９）式に従って、上限関数Ｄ0を小さくするように、基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を更新し、記憶部３０に格納する。
【００７９】
更に、パラメータ更新部２４は、（k、τ）の全ての組み合わせの各々について、記憶部３０に記憶されているｙk[t]、Ｂ(s)m,k、Ｇ(s)m[t]、ζk,m,t,τ、ηk[t]に基づいて、上記（２０）式に従って、上限関数Ｄ0を小さくするように、室内インパルス応答のパワーＨk[τ]を更新し、記憶部３０に格納する。
【００８０】
終了判定部２５は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、補助変数更新部２３及びパラメータ更新部２４の各処理を繰り返す。終了判定部２５は、終了条件を満足したと判定した場合には、信号変換部２６による処理に移行する。信号変換部２６は、記憶部３０に記憶されているＢ(s)m,k 及びＧ(s)m[t]に基づいて、雑音信号を取り除くことにより推定される原音声信号（以下、推定原音声信号と称する）に変換し、出力部４０により、推定原音声信号を出力する。
【００８１】
なお、終了条件としては、繰り返し回数がＬ-1回目の上限関数（１５）式の値と、繰り返し回数がＬ回目の上限関数（１５）式の値との差が、予め定めた閾値よりも小さくなったことを用いればよい。あるいは、終了条件として、繰り返し回数が、予め定められた上限回数に到達したことを用いてもよい。
【００８２】
＜音響信号解析装置の作用＞
【００８３】
次に、本実施の形態に係る音響信号解析装置１００の作用について説明する。まず、解析対象の信号として、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが音響信号解析装置１００に入力され、記憶部３０に格納される。そして、音響信号解析装置１００において、図３に示す音響信号解析処理ルーチンが実行される。
【００８４】
まず、ステップＳ１０１において、記憶部３０から、観測信号ｙ[t]を読み込み、当該観測信号ｙ[t]に対して、短時間フーリエ変換を用いた時間周波数分析を行い、観測信号ｙ[t]の観測時間周波数成分Ｙ={ｙk[t]}（k＝１，・・・，K、t＝１，・・・，T）を算出すると共に、得られた観測時間周波数成分Ｙを記憶部３０に記憶する。
【００８５】
そして、ステップＳ１０２において、乱数を用いて、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]の各初期値を設定して、記憶部３０に記憶すると共に、一定量のクリーン音声サンプルから予め学習して得た各周波数kのＭ個の基底スペクトルｍからなる基底行列の各要素をＢ(s)m,kとして設定して、記憶部３０に記憶する。
【００８６】
次にステップＳ１０３では、上記ステップＳ１０２で設定されたＢ(s)m,k、Ｇ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]、又は上記ステップＳ１０２で設定されたＢ(s)m,k及び後述するステップＳ１０４、Ｓ１０５、Ｓ１０６、Ｓ１０７で更新されたＧ(s)m[t]、Ｈk[τ]、Ｂ(n)q,k 、Ｇ(n)q[t]に基づいて、上記（２１）式に従って、補助変数ζk,m,t,τを各（k、m、t、τ）の組み合わせについて算出すると共に、上記（２２）式に従って、補助変数ξk,q,tを各（k、q、t）の組み合わせについて算出し、更に、上記（１６）式、（２３）式に従って、補助変数ηk[t]を各時刻t及び各周波数kについて算出して、記憶部３０に格納する。
【００８７】
ステップＳ１０４では、上記ステップＳ１０１で生成された観測信号ｙ[t]の観測時間周波数成分Ｙと、上記ステップＳ１０２で設定されたＢ(s)m,kと、上記ステップＳ１０２で設定されたＨk[τ]、又は後述するステップＳ１０７で更新されたＨk[τ]と、上記ステップＳ１０３で更新された補助変数ζk,m,t,τ、ηk[t]に基づいて、上記（１７）式に従って、基底スペクトルＢ(s)m,kの非負値のゲイン係数Ｇ(s)m[t]を各（m、t）の組み合わせについて算出して、記憶部３０に格納する。
【００８８】
ステップＳ１０５では、上記ステップＳ１０１で生成された観測信号ｙ[t]の観測時間周波数成分Ｙと、上記ステップＳ１０２で設定されたＧ(n)q[t]、又は後述するステップＳ１０６で更新されたＧ(n)q[t]と、上記ステップＳ１０３で更新された補助変数ξk,q,t、ηk[t]に基づいて、上記（１８）式に従って、雑音信号nのパワースペクトル密度Φ(n)k[t]の基底行列の要素Ｂ(n)q,kを各（q、k）の組み合わせについて算出して、記憶部３０に格納する。
【００８９】
ステップＳ１０６では、上記ステップＳ１０１で生成された観測信号ｙ[t]の観測時間周波数成分Ｙと、上記ステップＳ１０２で設定されたＢ(n)q,k、又は上記ステップＳ１０５で更新されたＢ(n)q,kと、上記ステップＳ１０３で更新された補助変数ξk,q,t、ηk[t]に基づいて、上記（１９）式に従って、基底スペクトルＢ(n)q,kのゲインＧ(n)q[t]を各（q、t）の組み合わせについて算出して、記憶部３０に格納する。
【００９０】
そして、ステップＳ１０７では、上記ステップＳ１０１で生成された観測信号ｙ[t]の観測時間周波数成分Ｙと、上記ステップＳ１０２で設定されたＢ(s)m,kと、上記ステップＳ１０２で設定されたＧ(s)m[t]、又は上記ステップＳ１０４で更新されたＧ(s)m[t]と、上記ステップＳ１０３で更新された補助変数ζk,m,t,τ、ηk[t]に基づいて、上記（２０）式に従って、室内インパルス応答のパワーＨk[τ]を各（k、τ）の組み合わせについて算出して、記憶部３０に格納する。
【００９１】
次のステップＳ１０８では、上記ステップＳ１０１で生成された観測信号ｙ[t]の観測時間周波数成分Ｙと、上記ステップＳ１０２で設定されたＢ(s)m,kと、上記ステップＳ１０３で更新された補助変数ζk,m,t,τ、ξk,q,t、ηk[t]と、上記ステップＳ１０４で更新されたＧ(s)m[t]と、上記ステップＳ１０５で更新されたＢ(n)q,kと、上記ステップＳ１０６で更新されたＧ(n)q[t]と、上記ステップＳ１０７で更新されたＨk[τ]に基づいて、上記（１５）式に従って上限関数Ｄ0の値を算出して、記憶部３０に記憶する。そして、前回のステップＳ１０８で算出した上限関数Ｄ0の値を記憶部３０から読み込み、今回のステップＳ１０８で算出した上限関数Ｄ0の値と、前回のステップＳ１０８で算出した上限関数Ｄ0の値との差分が、予め記憶部３０に記憶されている予め定められた閾値よりも小さいか否かを判定し、差分が予め定められた閾値以上の場合には、終了条件を満足していないと判断して、上記ステップＳ１０３へ戻り、上記ステップＳ１０３〜ステップＳ１０７の処理を繰り返す。一方、差分が予め定められた閾値未満の場合には、終了条件を満足したと判断して、ステップＳ１０９で、上記ステップＳ１０２で設定されたＢ(s)m,kと、上記ステップＳ１０４で最終的に更新されたＧ(s)m[t]に基づいて推定原音声信号を算出すると共に、出力部４０より推定原音声信号を出力して音響信号解析処理ルーチンを終了する。
【００９２】
＜実施結果＞
【００９３】
次に、本実施の形態に係る手法の有用性を示す目的で、残響室内（残響時間１．３秒）での移動音声信号に、ＰＨＳの着信音、背景雑音、当該ＰＨＳの着信音及び背景雑音を組み合わせた音、粒子落下音、ＢＧＭ音、他者の音声といった非定常雑音を重畳したものを観測信号として用いた場合の、提案法での評価実験の結果について説明する。
【００９４】
図４は、提案法を用いた場合の評価実験結果を示した図である。上記図４（Ａ）はメル周波数ケプストラム係数歪み（Mel-Frequency Cepstral Coefficients Distortion：ＭＦＣＣＤ）を評価基準とした実験結果、上記図４（Ｂ）はバークスペクトル歪みスコア（Bark Spectral Distortion Score：ＢＳＤＳ）を評価基準とした実験結果、上記図４（Ｃ）は信号対干渉比（Source to Interference Ratio：ＳＩＲ）を評価基準とした実験結果を示している。なお、ＭＦＣＣＤ及びＢＳＤＳは値が小さくなる程、音声の明瞭度が良くなることを示す値であり、一方、ＳＩＲは値が大きくなる程、音声の明瞭度が良くなることを示す値である。
【００９５】
上記図４から分かるように、移動音声信号にＢＧＭ又は他者の音声を重畳した場合のＭＦＣＣＤの評価結果を除いて、観測信号に含まれる音声の明瞭度に比べ提案法による推定原音声信号の明瞭度の方が向上するという結果を得られた。特に、移動音声信号にＰＨＳの着信音を重畳した場合のＳＩＲに関する評価結果では、４倍以上の改善が見られた。
【００９６】
以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、上記（１５）式で示される上限関数の値が小さくなるように、係数行列の各要素Ｇ(s)m[t]と、基底行列の各要素Ｂ(n)q,kと、係数行列の各要素Ｇ(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーＨk[τ]と、補助変数ζk,m,t,τと、補助変数ξk,q,tと、補助変数ηk[t]とを繰り返し更新することにより、音声信号に雑音及び残響成分が重畳した観測信号から精度よく音声信号を推定して、音声信号の明瞭度を向上させることができる。
【００９７】
このように、本発明で提案する手法では、音源の移動等に関して、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変とする観測信号の生成モデルを立てて、観測信号に含まれる原音声信号を高い精度で推定する
【００９８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００９９】
例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１００】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【０１０１】
また、本実施の形態に係る手法では、パワースペクトル密度系列モデルΦxk[t]と観測時間周波数成分ｙk[t]との距離の尺度として板倉斎藤距離Ｉ(Ｈ,Ｇ(s),Ｂ(n),Ｇ(n))を用いたが、これに限らず、β-divergenceの一種であるGeneralized Kullback-Leibler距離やFrobeniusノルムといった尺度を用いるようにしてもよい。
【符号の説明】
【０１０２】
１０   入力部
２０   演算部
２１   時間周波数解析部
２２   初期設定部
２３   補助変数更新部
２４   パラメータ更新部
２５   終了判定部
２６   信号変換部
３０   記憶部
４０   出力部
１００ 音響信号解析装置
【図２】
2015049406.tif 000017
【図３】
2015049406.tif 000018
【図１】
2015049406.tif 000019
【図４】
2015049406.tif 000020