特許第6142402号(P6142402)IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 日本電信電話株式会社の特許一覧 ▶ 国立大学法人 東京大学の特許一覧

特許6142402音響信号解析装置、方法、及びプログラム
<>
  • 特許6142402-音響信号解析装置、方法、及びプログラム 図000016
  • 特許6142402-音響信号解析装置、方法、及びプログラム 図000017
  • 特許6142402-音響信号解析装置、方法、及びプログラム 図000018
  • 特許6142402-音響信号解析装置、方法、及びプログラム 図000019
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6142402
(24)【登録日】2017年5月19日
(45)【発行日】2017年6月7日
(54)【発明の名称】音響信号解析装置、方法、及びプログラム
(51)【国際特許分類】
   G10L 21/0264 20130101AFI20170529BHJP
   G10L 25/27 20130101ALI20170529BHJP
   G10L 21/0232 20130101ALI20170529BHJP
   G10L 21/0308 20130101ALI20170529BHJP
【FI】
   G10L21/0264 Z
   G10L25/27
   G10L21/0232
   G10L21/0308 Z
【請求項の数】5
【全頁数】19
(21)【出願番号】特願2013-181701(P2013-181701)
(22)【出願日】2013年9月2日
(65)【公開番号】特開2015-49406(P2015-49406A)
(43)【公開日】2015年3月16日
【審査請求日】2015年7月9日
【新規性喪失の例外の表示】特許法第30条第2項適用 2013年3月5日 一般社団法人 電子情報通信学会 発行の「2013年総合大会講演論文集」に発表
(73)【特許権者】
【識別番号】000004226
【氏名又は名称】日本電信電話株式会社
(73)【特許権者】
【識別番号】504137912
【氏名又は名称】国立大学法人 東京大学
(74)【代理人】
【識別番号】110001519
【氏名又は名称】特許業務法人太陽国際特許事務所
(72)【発明者】
【氏名】亀岡 弘和
(72)【発明者】
【氏名】池澤 浩気
(72)【発明者】
【氏名】北条 伸克
(72)【発明者】
【氏名】エスピ ミケル
(72)【発明者】
【氏名】嵯峨山 茂樹
【審査官】 上田 雄
(56)【参考文献】
【文献】 特開2009−128906(JP,A)
【文献】 特開2012−027196(JP,A)
【文献】 池澤浩気,外4名,非定常雑音・時変残響環境下でのパワースペクトログラム領域セミブラインド音声強調,電子情報通信学会総合大会講演論文集,日本,一般社団法人電子情報通信学会,2013年 3月 5日,基礎・境界講演論文集,p.72
【文献】 安良岡直希,調波パラメトリックNMF及びIダイバージェンス規準残響推定に基づく音響信号モデリングとフレーズ置換への応用,修士論文,日本,京都大学大学院情報学研究科,2011年,pp.1-47
(58)【調査した分野】(Int.Cl.,DB名)
G10L 21/00−25/93
(57)【特許請求の範囲】
【請求項1】
音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分yk[t]を出力する時間周波数解析部と、
前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、M個の基底スペクトルmを表す非負値の各要素B(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルB(s)m,kのゲインG(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素G(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Q個の基底スペクトルqを表す非負値の各要素B(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルB(n)q,kのゲインG(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素B(n)q,k、前記係数行列の各要素G(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素B(s)m,kに、予め求められた値を設定するパラメータ初期値設定部と、
(τ、m)の全ての組み合わせにおける、前記基底行列の要素B(s)m,k、前記係数行列の要素G(s)m[t]、及び前記パワーHk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新するパラメータ更新部と、
予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、
を含む音響信号解析装置。
【請求項2】
前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離の尺度としてβ−divergenceを用いた
請求項1記載の音響信号解析装置。
【請求項3】
補助変数更新部を更に含み、
前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離を表す関数を、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]と、(k、m、t、τ)の全ての組み合わせについての補助変数ζk,m,t,τと、(k、q、t)の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]を用いて表され、かつ、前記パワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]とのβ−divergenceの上限関数である補助関数とし、
前記補助変数更新部は、前記補助関数を小さくするように、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、(k、m、t、τ)の全ての組み合わせについての補助変数ζk,m,t,τと、(k、q、t)の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]を更新し、
前記パラメータ更新部は、前記補助関数を小さくするように、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、(k、m、t、τ)の全ての組み合わせについての補助変数ζk,m,t,τと、(k、q、t)の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]とに基づいて、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新する
請求項2記載の音響信号解析装置。
【請求項4】
時間周波数解析部によって、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分yk[t]を出力し、
パラメータ初期値設定部によって、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、M個の基底スペクトルmを表す非負値の各要素B(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルB(s)m,kのゲインG(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素G(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Q個の基底スペクトルqを表す非負値の各要素B(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルB(n)q,kのゲインG(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素B(n)q,k、前記係数行列の各要素G(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素B(s)m,kに、予め求められた値を設定し、
パラメータ更新部によって、(τ、m)の全ての組み合わせにおける、前記基底行列の要素B(s)m,k、前記係数行列の要素G(s)m[t]、及び前記パワーHk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新し、
終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う
音響信号解析方法。
【請求項5】
コンピュータを、請求項1〜請求項3の何れか1項に記載の音響信号解析装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、音響信号の時系列データから、原音声信号を強調するための音響信号解析装置、方法、及びプログラムに関する。
【背景技術】
【0002】
音声強調技術は、例えば、ハンズフリー通話システムにおいて通信音声の了解性を向上させたり、自動会議録作成システムにおいて音声認識率を向上させたりするための重要な要素技術である。
【0003】
雑音や残響環境の情報は事前に知り得ないことが多いため、観測信号のみから音声信号を獲得する問題を扱う必要があるが、音声信号に残響と雑音が重畳されて観測信号が得られるプロセスを順問題と捉えると、観測信号から音声信号を得る問題は逆問題と見なすことができる。
【0004】
雑音や室内伝達系の情報が未知の場合、この逆問題には解が無数に存在しうるため、解を絞り込むための何らかの仮定が必要である。音声強調に関する従来アプローチでは、雑音に定常性、室内伝達系に時不変性などの制約を仮定し、その仮定の下で当該逆問題が定式化されるものが多かった。
【0005】
例えば、古典的な雑音除去手法であるスペクトラルサブトラクション法では、雑音の定常性が仮定され、既知の非音声区間から得られる雑音のパワー(または振幅)スペクトルを観測信号のパワー(または振幅)スペクトルから減算することで音声信号のパワー(または振幅)スペクトルを推定する。また、逆フィルタリングに基づく残響除去法では、時不変かつ最小位相な室内伝達系が仮定され、音声に関する仮定やモデル(調波性、自己回帰モデル、自己相関関数コードブック等)を手がかりに最適な残響除去フィルタを推定する。
【0006】
こうした従来手法は、仮定が成立する状況においては高い性能を発揮する一方で、実環境においては、雑音が例えば携帯電話の着信音及び背景音楽などのように非定常なものである場合や、僅かでも音源移動がある場合等では、雑音や室内伝達系に対する上述の仮定が成立しない場合があり、このような環境の下では高い性能を発揮できなくなるという問題があった。
【0007】
こうした室内伝達系の変化に対する頑健な残響除去の手法として、音声スペクトログラムと室内インパルス応答のスペクトログラムの時間方向の畳み込みモデルに基づき、音声スペクトログラムのスパース性を手がかりに音声スペクトログラムを推定する方法が提案されている(非特許文献1)。
【0008】
しかし、この手法では残響のような乗法的な雑音の混入プロセスしか想定されておらず、加法的な雑音の混入プロセスは陽に想定されていない(雑音が存在しない状況を仮定している)ため、観測信号に雑音が含まれる場合には高い性能を発揮できないという問題があった。
【0009】
一方、非定常な雑音を抑圧するための手法として、スペクトログラムを非負値行列と見立て、これを二つの非負値行列に分解する非負値行列因子分解アプローチが近年注目されている(非特許文献2)。非負値行列因子分解を用いれば、観測スペクトログラムを構成する基底スペクトルと各時刻におけるそれらの結合係数を同時に得ることができる。これにより、例えばクリーンな音声サンプルからあらかじめ基底スペクトルを獲得(学習)しておき、雑音つき音声の観測スペクトログラムを基底行列と係数行列に分解する際、基底スペクトルセットの一部をあらかじめ学習した基底スペクトルのセットに固定しておくことで、音声とそれ以外(すなわち雑音)のスペクトログラムに分解する。
【先行技術文献】
【非特許文献】
【0010】
【非特許文献1】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.
【非特許文献2】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.
【発明の概要】
【発明が解決しようとする課題】
【0011】
しかし、上記の非特許文献2の手法では加法的な雑音の混入プロセスしか想定されておらず、残響のような乗法的な雑音の混入プロセスは陽に想定されていない(残響がない状況を仮定している)ため、観測信号に残響が含まれる場合には高い性能を発揮できないという問題があった。
【0012】
本発明は、上記の事情を鑑みてなされたもので、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【0013】
上記の目的を達成するために本発明に係る音響信号解析装置は、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分yk[t]を出力する時間周波数解析部と、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、M個の基底スペクトルmを表す非負値の各要素B(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルB(s)m,kのゲインG(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素G(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Q個の基底スペクトルqを表す非負値の各要素B(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルB(n)q,kのゲインG(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素B(n)q,k、前記係数行列の各要素G(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素B(s)m,kに、予め求められた値を設定するパラメータ初期値設定部と、(τ、m)の全ての組み合わせにおける、前記基底行列の要素B(s)m,k、前記係数行列の要素G(s)m[t]、及び前記パワーHk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新するパラメータ更新部と、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、を含んで構成されている。
【0014】
本発明に係る音響信号解析方法は、時間周波数解析部によって、音響信号の時系列データを入力として、各時刻tにおける各周波数kの観測時間周波数成分yk[t]を出力し、パラメータ初期値設定部によって、前記音響信号に含まれる原音声信号sの各時刻t及び各周波数kのパワースペクトル密度Φ(s)k[t]を、M個の基底スペクトルmを表す非負値の各要素B(s)m,kからなる基底行列と、各時刻tにおける基底スペクトルB(s)m,kのゲインG(s)m[t]を非負値である各要素とする係数行列との積として表した場合の前記係数行列の各要素G(s)m[t]、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]、前記音響信号に含まれる雑音信号nの各時刻t及び各周波数kのパワースペクトル密度Φ(n)k[t]を、Q個の基底スペクトルqを表す非負値の各要素B(n)q,kからなる基底行列と、各時刻tにおける基底スペクトルB(n)q,kのゲインG(n)q[t]を非負値である各要素とする係数行列との積として表した場合の前記基底行列の各要素B(n)q,k、前記係数行列の各要素G(n)q[t]の各々の初期値を設定すると共に、前記基底行列の各要素B(s)m,kに、予め求められた値を設定し、パラメータ更新部によって、(τ、m)の全ての組み合わせにおける、前記基底行列の要素B(s)m,k、前記係数行列の要素G(s)m[t]、及び前記パワーHk[τ]と、前記パワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、前記観測時間周波数成分yk[t]との距離が小さくなるように、各時刻t及び各周波数kの前記観測時間周波数成分yk[t]と、前記基底行列の各要素B(s)m,kと、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦxk[t]とに基づいて、前記係数行列の各要素G(s)m[t]と、前記基底行列の各要素B(n)q,kと、前記係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新し、終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う。
【0015】
本発明に係るプログラムは、上記の音響信号解析装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【0016】
以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、原音声信号sのパワースペクトル密度Φ(s)k[t]の基底行列及び係数行列と、室内インパルス応答のパワーHk[τ]と、雑音信号nのパワースペクトル密度Φ(n)k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦxk[t]と、観測時間周波数成分yk[t]との距離が小さくなるように、パワースペクトル密度Φ(s)k[t]の係数行列と、パワースペクトル密度Φ(n)k[t]の基底行列及び係数行列と、室内インパルス応答のパワーHk[τ]を更新することを繰り返すことにより、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる、という効果が得られる。
【図面の簡単な説明】
【0017】
図1】マイクロホンの位置を変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分及び位相成分を示す図である。
図2】本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。
図3】本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。
図4】(A)メル周波数ケプストラム係数歪みの評価結果を示すグラフ、(B)バークスペクトル歪みスコアの評価結果を示すグラフ、及び(C)信号対干渉比の評価結果を示すグラフである。
【発明を実施するための形態】
【0018】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【0019】
<発明の原理>
【0020】
まず、本発明の原理について説明する。
【0021】
<音声強調問題の定式化>
【0022】
室内伝達系が時不変の場合、原音声信号をs[t]、室内インパルス応答をh[t]、雑音信号をn[t] とすると、音響信号としての観測信号y[t]は時間領域において、以下の(1)式のように表される。
【0023】
【数1】
【0024】
ただし、tは離散時刻のインデックスである。右辺第1項Στs[t-τ]h[τ]は残響音声信号を表す。ここで、観測信号、原信号、室内インパルス応答、雑音信号の短時間フーリエ変換をそれぞれyk[t]、sk[t]、hk[t]、nk[t]とする。ただし、k=1,...,Kは周波数のインデックス、t=1,...,Tはフレーム(時刻)のインデックスを表す。時間領域における信号同士の畳み込みは、特定の条件の下でそれぞれの信号の各周波数における狭帯域信号同士の畳み込みで近似できることが知られている。もしこの近似が成立するような短時間フーリエ変換の分析条件を選べば、yk[t]は(2)式のように表される。
【0025】
【数2】
【0026】
以上では時不変な系を考えたが、実際には音源が移動することがあるため、室内伝達系の時不変性の仮定は必ずしも成り立たない。ここで、時変な室内伝達系を仮定することは、(2)式において室内インパルス応答hk[τ]が時刻tごとに変化しうることを許容することに相当する。すなわち、以下の(3)式のように表される。
【0027】
【数3】
【0028】
非定常雑音、時変残響環境下における音声強調問題は、(3)式をもとに観測時間周波数成分Y={yk[t]}1≦k≦K,1≦t≦TからS={sk[t]}1≦k≦K,1≦t≦Tを推定する問題として定式化される。(3)式では室内インパルス応答が時刻ごとに変化することを許容しているため、(2)式に比べて未知パラメータが圧倒的に多くなっており、同じ右辺を与えるsとhとnの組み合わせは無数に存在する。そこで以下では、その無数の組み合わせの中から正しい解を得るためにどのような仮定を置くべきか、について説明する。
【0029】
<室内伝達系の「半時変性」>
【0030】
室内インパルス応答が時間的に変化する主要因の1つは話者の移動である。なぜなら、話者の移動に伴い、話者からマイクロホンまでの(壁からの反射音の到来経路を含む)あらゆる音響経路の距離が一斉に変化するからである。
【0031】
ここで、単一の音響経路のみに着目すると、伝達系のインパルス応答はデルタ関数となるが、経路の距離に変化があった場合、当該インパルス応答のフーリエ変換の位相成分には音響経路の距離変化に伴う到来時間の変化、振幅成分には音響経路の距離変化に伴う強度の変化がそれぞれ反映されることになる。音の強度は伝播距離に対して反比例的であるため、伝播距離がある程度大きい場合には伝播距離に変化があっても振幅成分はほとんど変化しないことになる。
【0032】
一方、位相成分については、到来時間と周波数に比例した量だけ変化するため、到来時間の変化がたとえわずかであっても、特に高周波域においては著しく変化しうる。
【0033】
よって、単一の音響経路の伝達系のインパルス応答の短時間フーリエ変換には、話者の移動に伴って振幅は変化しにくく位相は変化しやすい、という性質があることが予想される。
【0034】
実際の室内インパルス応答はあらゆる音響経路の伝達系のインパルス応答の重ね合わせとなるので、一つ一つのインパルス応答が上述のような傾向があるならば、その重ね合わせもおおよそ同様な傾向があるのではないかと考えられる。
【0035】
図1は、話者の代わりにマイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分と位相成分を比較したものである。
【0036】
上記図1の上段は、室内インパルス応答の短時間フーリエ変換の振幅成分(A)及び位相成分(B)を示したものであり、上記図1の下段は、マイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分(C)と位相成分(D)を示したものである。なお、室内インパルス応答のデータはRWCP実環境音声・音響データベースのマイクロホンアレーデータベースに収録されている、残響時間0.3[s]の残響可変室(パネル)で測定されたものである。
【0037】
上記図1から、振幅成分はさほど変化していないのに対し位相成分は著しく変化していることが分かる。
【0038】
以上より、発明者らは、話者の移動等、環境の軽微な変化については、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変と仮定した特殊な系により扱えるのではないかと考えた。以後、このような系を「半時変系」と呼ぶこととする。
【0039】
<観測信号の生成モデル>
【0040】
以上の2つの仮定をベースに、観測信号の生成モデルを定式化する。(3)式に半時変性の仮定を導入すると、観測信号の短時間フーリエ変換は(4)式のように表される。
【0041】
【数4】
【0042】
ただし、φk,t[τ]とHk[t]は室内インパルス応答の短時間フーリエ変換の位相成分と振幅成分であり、Hk[t]=|hk,t[τ]|である。ここで、φk,t[τ]をk,n,mごとに独立に区間[0;2π)上の一様分布に従う確率変数とする。さらに、原音声信号、雑音については、それぞれ平均が0の複素正規分布(5)式及び(6)式に従う確率変数とする。
【0043】
【数5】
【0044】
ただし、Φ(s)k[t]、Φ(n)k[t]はそれぞれ原音声信号および雑音信号の時刻t、k番目の周波数ビンにおけるパワースペクトル密度である。ここで、NC(z;0,λ)は複素正規分布の確率密度関数を表し、(7)式によって与えられる。
【0045】
【数6】
【0046】
証明は省略するが、原音声信号と雑音のガウス性の仮定と、正規分布の再生性より、観測信号の短時間フーリエ変換yk[t]はk、tごとに独立に、(8)式に示す複素正規分布に従うことが示される。
【0047】
【数7】
【0048】
次に、原音声信号のパワースペクトル密度系列Φ(s)k[t]に対し、非負値行列積表現を導入する。基底スペクトルをM個とし、m番目の基底スペクトルをB(s)m,kとする。また、時刻tにおける基底スペクトルB(s)m,kのゲインをG(s)m[t]とすると、原音声信号のパワースペクトル密度系列は(9)式及び(10)式で示される。
【0049】
【数8】
【0050】
なお、基底スペクトルとそのゲインは非負であることに注意する.また、スケールの任意性を除くため、(11)式に示す条件を満たすものとする。
【0051】
【数9】
【0052】
音声の基底スペクトルB(s)m,kは未知変数として観測信号から他の未知パラメータとしてブラインドで推定することもありえるが、ある程度の量のクリーン音声サンプルからあらかじめ学習しておくことを想定する。その意味で、本実施の形態に係る手法(以下、提案法と称する)はセミブラインドな音声強調手法である。基底の事前学習は、クリーン音声データのスペクトログラムに対してNMFを適用することで行うことができる。
【0053】
<パラメータ推定アルゴリズム>
【0054】
以上で立てた観測時間周波数成分Y={yk[t]}1≦k≦K,1≦t≦Tの確率密度関数は、未知パラメータの尤度関数に対応する。(8)式より、観測時間周波数成分Yが与えられたもとでの未知パラメータの対数尤度は、定数項・符号を無視すれば、観測スペクトログラム|yk[t]|2と(12)式で示されるパワースペクトル密度系列モデルとの距離を表す板倉斎藤距離(13)式と等しくなる。
【0055】
【数10】
【0056】
従って、提案法は観測スペクトログラム|yk[t]|2とパワースペクトル密度系列モデルΦxk[t]との最適フィッティング問題に帰着する。ここで、雑音信号のパワースペクトル密度系列Φ(n)k[t]に関しても、(14)式のようにモデル化される。
【0057】
【数11】
【0058】
なお、板倉斎藤距離I(H,G(s),B(n),G(n))は、距離尺度を表すβ-divergenceにおいてβを0にした場合に得られる距離に相当する。
【0059】
(13)式で示される板倉斎藤距離I(H,G(s),B(n),G(n))は直接最小化することが困難であるため、板倉斎藤距離I(H,G(s),B(n),G(n))を下回らず1点で接する上限関数を設定する。導出は省略するが、(15)式に示す上限関数が板倉斎藤距離I(H,G(s),B(n),G(n))の上限関数となる。
【0060】
【数12】
【0061】
すなわち、板倉斎藤距離I(H,G(s),B(n),G(n))の代わりに(15)式に示された上限関数を小さくするようにH,G(s),B(n),G(n)および補助変数ζ,ξ,ηを交互に更新することで目的関数I(H,G(s),B(n),G(n))を小さくしていくことができる。この際、H,G(s),B(n),G(n)の各パラメータは非負値を保ったまま更新するようにする。なお、(15)式におけるΦ(x)k[t]は(16)式によって示される。
【0062】
【数13】
【0063】
そして、この上限関数D0の更新則は(17)式〜(23)式で示される。
【0064】
【数14】

<システム構成>
【0065】
次に、原音声信号に雑音及び残響成分が重畳された観測信号を解析して、観測信号に含まれる原音声信号を推定する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【0066】
図2に示すように、本発明の実施の形態に係る音響信号解析装置は、CPUと、RAMと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したROMとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【0067】
音響信号解析装置100は、入力部10と、演算部20と、記憶部30と、出力部40とを備えている。
【0068】
入力部10により、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが入力される。記憶部30は、入力部10により入力された観測信号の時系列データを記憶する。また、記憶部30は、後述する各処理での結果を記憶すると共に、本処理ルーチンで用いる各パラメータの初期値を記憶している。
【0069】
演算部20は、時間周波数解析部21と、初期設定部22と、補助変数更新部23と、パラメータ更新部24と、終了判定部25と、信号変換部26とを備えている。
【0070】
時間周波数解析部21は、例えばマイクロホンの時系列信号としての観測された観測信号y[t]を入力として、観測信号y[t]の観測時間周波数成分Y={yk[t]}(周波数k=1,・・・,K、時刻t=1,・・・,T)を計算する。また、計算した観測時間周波数成分Yを、記憶部30に記憶しておく。より詳細には、時間周波数解析部21は、例えばマイクロホンで観測された観測信号の時系列データを入力として、短時間フーリエ変換(Short-Time Fourier Transform;STFT)を用いて時間周波数解析を行うことにより観測時間周波数成分Yを計算する。
【0071】
初期設定部22は、後述する処理で用いる各パラメータG(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]の各初期値及びB(s)m,k の値を設定する。なお、B(s)m,k以外のパラメータの初期値は、例えば乱数を用いて適当な値に設定すればよい。この場合、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]の各パラメータの初期値は非負値となるように設定する。
【0072】
そして、B(s)m,kの値には、一定量のクリーン音声サンプルから予め学習して得た値を設定する。なお、スケールの任意性を排除するため、B(s)m,kの値は上記(11)式を満たすように予め学習されているものとする。
【0073】
補助変数更新部23は、(k、m、t、τ)の全ての組み合わせの各々について、記憶部30に記憶されているB(s)m,k、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]に基づいて、上記(21)式に従って上限関数D0を小さくするように補助変数ζk,m,t,τを更新し、記憶部30に格納する。
【0074】
また、補助変数更新部23は、(k、q、t)の全ての組み合わせの各々について、記憶部30に記憶されているB(s)m,k、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]に基づいて、上記(22)式に従って上限関数D0を小さくするように補助変数ξk,q,tを更新し、記憶部30に格納する。
【0075】
更に、補助変数更新部23は、各時刻t及び各周波数kについて、記憶部30に記憶されているB(s)m,k、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]に基づいて、上記(16)式、上記(23)式に従って上限関数D0を小さくするように補助変数ηk[t]を更新し、記憶部30に格納する。
【0076】
パラメータ更新部24は、(m、t)の全ての組み合わせの各々について、記憶部30に記憶されているyk[t]、B(s)m,k、Hk[τ]、ζk,m,t,τ、ηk[t]に基づいて、上記(17)式に従って、上限関数D0を小さくするように、基底スペクトルB(s)m,kのゲイン係数G(s)m[t]を更新し、記憶部30に格納する。
【0077】
また、パラメータ更新部24は、(q、k)の全ての組み合わせの各々について、記憶部30に記憶されているyk[t]、G(n)q[t]、ξk,q,t、ηk[t]に基づいて、上記(18)式に従って、上限関数D0を小さくするように、雑音信号nのパワースペクトル密度Φ(n)k[t]の基底行列の要素B(n)q,kを更新し、記憶部30に格納する。
【0078】
また、パラメータ更新部24は、(q、t)の全ての組み合わせの各々について、記憶部30に記憶されているyk[t]、B(n)q,k、ξk,q,t、ηk[t]に基づいて、上記(19)式に従って、上限関数D0を小さくするように、基底スペクトルB(n)q,kのゲインG(n)q[t]を更新し、記憶部30に格納する。
【0079】
更に、パラメータ更新部24は、(k、τ)の全ての組み合わせの各々について、記憶部30に記憶されているyk[t]、B(s)m,k、G(s)m[t]、ζk,m,t,τ、ηk[t]に基づいて、上記(20)式に従って、上限関数D0を小さくするように、室内インパルス応答のパワーHk[τ]を更新し、記憶部30に格納する。
【0080】
終了判定部25は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、補助変数更新部23及びパラメータ更新部24の各処理を繰り返す。終了判定部25は、終了条件を満足したと判定した場合には、信号変換部26による処理に移行する。信号変換部26は、記憶部30に記憶されているB(s)m,k 及びG(s)m[t]に基づいて、雑音信号を取り除くことにより推定される原音声信号(以下、推定原音声信号と称する)に変換し、出力部40により、推定原音声信号を出力する。
【0081】
なお、終了条件としては、繰り返し回数がL-1回目の上限関数(15)式の値と、繰り返し回数がL回目の上限関数(15)式の値との差が、予め定めた閾値よりも小さくなったことを用いればよい。あるいは、終了条件として、繰り返し回数が、予め定められた上限回数に到達したことを用いてもよい。
【0082】
<音響信号解析装置の作用>
【0083】
次に、本実施の形態に係る音響信号解析装置100の作用について説明する。まず、解析対象の信号として、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが音響信号解析装置100に入力され、記憶部30に格納される。そして、音響信号解析装置100において、図3に示す音響信号解析処理ルーチンが実行される。
【0084】
まず、ステップS101において、記憶部30から、観測信号y[t]を読み込み、当該観測信号y[t]に対して、短時間フーリエ変換を用いた時間周波数分析を行い、観測信号y[t]の観測時間周波数成分Y={yk[t]}(k=1,・・・,K、t=1,・・・,T)を算出すると共に、得られた観測時間周波数成分Yを記憶部30に記憶する。
【0085】
そして、ステップS102において、乱数を用いて、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]の各初期値を設定して、記憶部30に記憶すると共に、一定量のクリーン音声サンプルから予め学習して得た各周波数kのM個の基底スペクトルmからなる基底行列の各要素をB(s)m,kとして設定して、記憶部30に記憶する。
【0086】
次にステップS103では、上記ステップS102で設定されたB(s)m,k、G(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]、又は上記ステップS102で設定されたB(s)m,k及び後述するステップS104、S105、S106、S107で更新されたG(s)m[t]、Hk[τ]、B(n)q,k 、G(n)q[t]に基づいて、上記(21)式に従って、補助変数ζk,m,t,τを各(k、m、t、τ)の組み合わせについて算出すると共に、上記(22)式に従って、補助変数ξk,q,tを各(k、q、t)の組み合わせについて算出し、更に、上記(16)式、(23)式に従って、補助変数ηk[t]を各時刻t及び各周波数kについて算出して、記憶部30に格納する。
【0087】
ステップS104では、上記ステップS101で生成された観測信号y[t]の観測時間周波数成分Yと、上記ステップS102で設定されたB(s)m,kと、上記ステップS102で設定されたHk[τ]、又は後述するステップS107で更新されたHk[τ]と、上記ステップS103で更新された補助変数ζk,m,t,τ、ηk[t]に基づいて、上記(17)式に従って、基底スペクトルB(s)m,kの非負値のゲイン係数G(s)m[t]を各(m、t)の組み合わせについて算出して、記憶部30に格納する。
【0088】
ステップS105では、上記ステップS101で生成された観測信号y[t]の観測時間周波数成分Yと、上記ステップS102で設定されたG(n)q[t]、又は後述するステップS106で更新されたG(n)q[t]と、上記ステップS103で更新された補助変数ξk,q,t、ηk[t]に基づいて、上記(18)式に従って、雑音信号nのパワースペクトル密度Φ(n)k[t]の基底行列の要素B(n)q,kを各(q、k)の組み合わせについて算出して、記憶部30に格納する。
【0089】
ステップS106では、上記ステップS101で生成された観測信号y[t]の観測時間周波数成分Yと、上記ステップS102で設定されたB(n)q,k、又は上記ステップS105で更新されたB(n)q,kと、上記ステップS103で更新された補助変数ξk,q,t、ηk[t]に基づいて、上記(19)式に従って、基底スペクトルB(n)q,kのゲインG(n)q[t]を各(q、t)の組み合わせについて算出して、記憶部30に格納する。
【0090】
そして、ステップS107では、上記ステップS101で生成された観測信号y[t]の観測時間周波数成分Yと、上記ステップS102で設定されたB(s)m,kと、上記ステップS102で設定されたG(s)m[t]、又は上記ステップS104で更新されたG(s)m[t]と、上記ステップS103で更新された補助変数ζk,m,t,τ、ηk[t]に基づいて、上記(20)式に従って、室内インパルス応答のパワーHk[τ]を各(k、τ)の組み合わせについて算出して、記憶部30に格納する。
【0091】
次のステップS108では、上記ステップS101で生成された観測信号y[t]の観測時間周波数成分Yと、上記ステップS102で設定されたB(s)m,kと、上記ステップS103で更新された補助変数ζk,m,t,τ、ξk,q,t、ηk[t]と、上記ステップS104で更新されたG(s)m[t]と、上記ステップS105で更新されたB(n)q,kと、上記ステップS106で更新されたG(n)q[t]と、上記ステップS107で更新されたHk[τ]に基づいて、上記(15)式に従って上限関数D0の値を算出して、記憶部30に記憶する。そして、前回のステップS108で算出した上限関数D0の値を記憶部30から読み込み、今回のステップS108で算出した上限関数D0の値と、前回のステップS108で算出した上限関数D0の値との差分が、予め記憶部30に記憶されている予め定められた閾値よりも小さいか否かを判定し、差分が予め定められた閾値以上の場合には、終了条件を満足していないと判断して、上記ステップS103へ戻り、上記ステップS103〜ステップS107の処理を繰り返す。一方、差分が予め定められた閾値未満の場合には、終了条件を満足したと判断して、ステップS109で、上記ステップS102で設定されたB(s)m,kと、上記ステップS104で最終的に更新されたG(s)m[t]に基づいて推定原音声信号を算出すると共に、出力部40より推定原音声信号を出力して音響信号解析処理ルーチンを終了する。
【0092】
<実施結果>
【0093】
次に、本実施の形態に係る手法の有用性を示す目的で、残響室内(残響時間1.3秒)での移動音声信号に、PHSの着信音、背景雑音、当該PHSの着信音及び背景雑音を組み合わせた音、粒子落下音、BGM音、他者の音声といった非定常雑音を重畳したものを観測信号として用いた場合の、提案法での評価実験の結果について説明する。
【0094】
図4は、提案法を用いた場合の評価実験結果を示した図である。上記図4(A)はメル周波数ケプストラム係数歪み(Mel-Frequency Cepstral Coefficients Distortion:MFCCD)を評価基準とした実験結果、上記図4(B)はバークスペクトル歪みスコア(Bark Spectral Distortion Score:BSDS)を評価基準とした実験結果、上記図4(C)は信号対干渉比(Source to Interference Ratio:SIR)を評価基準とした実験結果を示している。なお、MFCCD及びBSDSは値が小さくなる程、音声の明瞭度が良くなることを示す値であり、一方、SIRは値が大きくなる程、音声の明瞭度が良くなることを示す値である。
【0095】
上記図4から分かるように、移動音声信号にBGM又は他者の音声を重畳した場合のMFCCDの評価結果を除いて、観測信号に含まれる音声の明瞭度に比べ提案法による推定原音声信号の明瞭度の方が向上するという結果を得られた。特に、移動音声信号にPHSの着信音を重畳した場合のSIRに関する評価結果では、4倍以上の改善が見られた。
【0096】
以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、上記(15)式で示される上限関数の値が小さくなるように、係数行列の各要素G(s)m[t]と、基底行列の各要素B(n)q,kと、係数行列の各要素G(n)q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、補助変数ζk,m,t,τと、補助変数ξk,q,tと、補助変数ηk[t]とを繰り返し更新することにより、音声信号に雑音及び残響成分が重畳した観測信号から精度よく音声信号を推定して、音声信号の明瞭度を向上させることができる。
【0097】
このように、本発明で提案する手法では、音源の移動等に関して、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変とする観測信号の生成モデルを立てて、観測信号に含まれる原音声信号を高い精度で推定する
【0098】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【0099】
例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、WWWシステムを利用している場合であれば、ホームページ提供環境(あるいは表示環境)も含むものとする。
【0100】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【0101】
また、本実施の形態に係る手法では、パワースペクトル密度系列モデルΦxk[t]と観測時間周波数成分yk[t]との距離の尺度として板倉斎藤距離I(H,G(s),B(n),G(n))を用いたが、これに限らず、β-divergenceの一種であるGeneralized Kullback-Leibler距離やFrobeniusノルムといった尺度を用いるようにしてもよい。
【符号の説明】
【0102】
10 入力部
20 演算部
21 時間周波数解析部
22 初期設定部
23 補助変数更新部
24 パラメータ更新部
25 終了判定部
26 信号変換部
30 記憶部
40 出力部
100 音響信号解析装置
図1
図2
図3
図4