JP6142402B2 - 音響信号解析装置、方法、及びプログラム - Google Patents

音響信号解析装置、方法、及びプログラム Download PDF

Info

Publication number
JP6142402B2
JP6142402B2 JP2013181701A JP2013181701A JP6142402B2 JP 6142402 B2 JP6142402 B2 JP 6142402B2 JP 2013181701 A JP2013181701 A JP 2013181701A JP 2013181701 A JP2013181701 A JP 2013181701A JP 6142402 B2 JP6142402 B2 JP 6142402B2
Authority
JP
Japan
Prior art keywords
time
signal
unit
frequency
acoustic signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2013181701A
Other languages
English (en)
Other versions
JP2015049406A (ja
Inventor
弘和 亀岡
弘和 亀岡
浩気 池澤
浩気 池澤
伸克 北条
伸克 北条
ミケル エスピ
ミケル エスピ
茂樹 嵯峨山
茂樹 嵯峨山
Original Assignee
日本電信電話株式会社
国立大学法人 東京大学
国立大学法人 東京大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 日本電信電話株式会社, 国立大学法人 東京大学, 国立大学法人 東京大学 filed Critical 日本電信電話株式会社
Priority to JP2013181701A priority Critical patent/JP6142402B2/ja
Publication of JP2015049406A publication Critical patent/JP2015049406A/ja
Application granted granted Critical
Publication of JP6142402B2 publication Critical patent/JP6142402B2/ja
Application status is Active legal-status Critical
Anticipated expiration legal-status Critical

Links

Images

Description

本発明は、音響信号解析装置、方法、及びプログラムに係り、特に、音響信号の時系列データから、原音声信号を強調するための音響信号解析装置、方法、及びプログラムに関する。

音声強調技術は、例えば、ハンズフリー通話システムにおいて通信音声の了解性を向上させたり、自動会議録作成システムにおいて音声認識率を向上させたりするための重要な要素技術である。

雑音や残響環境の情報は事前に知り得ないことが多いため、観測信号のみから音声信号を獲得する問題を扱う必要があるが、音声信号に残響と雑音が重畳されて観測信号が得られるプロセスを順問題と捉えると、観測信号から音声信号を得る問題は逆問題と見なすことができる。

雑音や室内伝達系の情報が未知の場合、この逆問題には解が無数に存在しうるため、解を絞り込むための何らかの仮定が必要である。音声強調に関する従来アプローチでは、雑音に定常性、室内伝達系に時不変性などの制約を仮定し、その仮定の下で当該逆問題が定式化されるものが多かった。

例えば、古典的な雑音除去手法であるスペクトラルサブトラクション法では、雑音の定常性が仮定され、既知の非音声区間から得られる雑音のパワー(または振幅)スペクトルを観測信号のパワー(または振幅)スペクトルから減算することで音声信号のパワー(または振幅)スペクトルを推定する。また、逆フィルタリングに基づく残響除去法では、時不変かつ最小位相な室内伝達系が仮定され、音声に関する仮定やモデル(調波性、自己回帰モデル、自己相関関数コードブック等)を手がかりに最適な残響除去フィルタを推定する。

こうした従来手法は、仮定が成立する状況においては高い性能を発揮する一方で、実環境においては、雑音が例えば携帯電話の着信音及び背景音楽などのように非定常なものである場合や、僅かでも音源移動がある場合等では、雑音や室内伝達系に対する上述の仮定が成立しない場合があり、このような環境の下では高い性能を発揮できなくなるという問題があった。

こうした室内伝達系の変化に対する頑健な残響除去の手法として、音声スペクトログラムと室内インパルス応答のスペクトログラムの時間方向の畳み込みモデルに基づき、音声スペクトログラムのスパース性を手がかりに音声スペクトログラムを推定する方法が提案されている(非特許文献1)。

しかし、この手法では残響のような乗法的な雑音の混入プロセスしか想定されておらず、加法的な雑音の混入プロセスは陽に想定されていない(雑音が存在しない状況を仮定している)ため、観測信号に雑音が含まれる場合には高い性能を発揮できないという問題があった。

一方、非定常な雑音を抑圧するための手法として、スペクトログラムを非負値行列と見立て、これを二つの非負値行列に分解する非負値行列因子分解アプローチが近年注目されている(非特許文献2)。非負値行列因子分解を用いれば、観測スペクトログラムを構成する基底スペクトルと各時刻におけるそれらの結合係数を同時に得ることができる。これにより、例えばクリーンな音声サンプルからあらかじめ基底スペクトルを獲得(学習)しておき、雑音つき音声の観測スペクトログラムを基底行列と係数行列に分解する際、基底スペクトルセットの一部をあらかじめ学習した基底スペクトルのセットに固定しておくことで、音声とそれ以外(すなわち雑音)のスペクトログラムに分解する。

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.

しかし、上記の非特許文献2の手法では加法的な雑音の混入プロセスしか想定されておらず、残響のような乗法的な雑音の混入プロセスは陽に想定されていない(残響がない状況を仮定している)ため、観測信号に残響が含まれる場合には高い性能を発揮できないという問題があった。

本発明は、上記の事情を鑑みてなされたもので、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる音響信号解析装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る音響信号解析装置は、音響信号の時系列データを入力として、各時刻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]とに基づいて算出されるパワースペクトル密度系列モデルΦx k[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のパワースペクトル密度系列モデルΦx k[t]とに基づいて、前記係数行列の各要素G(s) m[t]と、前記基底行列の各要素B(n) q,kと、前記係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新するパラメータ更新部と、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、を含んで構成されている。

本発明に係る音響信号解析方法は、時間周波数解析部によって、音響信号の時系列データを入力として、各時刻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]とに基づいて算出されるパワースペクトル密度系列モデルΦx k[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のパワースペクトル密度系列モデルΦx k[t]とに基づいて、前記係数行列の各要素G(s) m[t]と、前記基底行列の各要素B(n) q,kと、前記係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新し、終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う。

本発明に係るプログラムは、上記の音響信号解析装置の各部としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の音響信号解析装置、方法、及びプログラムによれば、原音声信号sのパワースペクトル密度Φ(s) k[t]の基底行列及び係数行列と、室内インパルス応答のパワーHk[τ]と、雑音信号nのパワースペクトル密度Φ(n) k[t]とに基づいて算出されるパワースペクトル密度系列モデルΦx k[t]と、観測時間周波数成分yk[t]との距離が小さくなるように、パワースペクトル密度Φ(s) k[t]の係数行列と、パワースペクトル密度Φ(n) k[t]の基底行列及び係数行列と、室内インパルス応答のパワーHk[τ]を更新することを繰り返すことにより、原音声信号に雑音及び残響成分が重畳した音響信号の時系列データから原音声信号を精度よく分離して、原音声信号を強調させることができる、という効果が得られる。

マイクロホンの位置を変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分及び位相成分を示す図である。 本発明の実施の形態に係る音響信号解析装置の構成を示す概略図である。 本発明の実施の形態に係る音響信号解析装置における音響信号解析処理ルーチンの内容を示すフローチャートである。 (A)メル周波数ケプストラム係数歪みの評価結果を示すグラフ、(B)バークスペクトル歪みスコアの評価結果を示すグラフ、及び(C)信号対干渉比の評価結果を示すグラフである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

<発明の原理>

まず、本発明の原理について説明する。

<音声強調問題の定式化>

室内伝達系が時不変の場合、原音声信号をs[t]、室内インパルス応答をh[t]、雑音信号をn[t] とすると、音響信号としての観測信号y[t]は時間領域において、以下の(1)式のように表される。

ただし、tは離散時刻のインデックスである。右辺第1項Στs[t-τ]h[τ]は残響音声信号を表す。ここで、観測信号、原信号、室内インパルス応答、雑音信号の短時間フーリエ変換をそれぞれyk[t]、sk[t]、hk[t]、nk[t]とする。ただし、k=1,...,Kは周波数のインデックス、t=1,...,Tはフレーム(時刻)のインデックスを表す。時間領域における信号同士の畳み込みは、特定の条件の下でそれぞれの信号の各周波数における狭帯域信号同士の畳み込みで近似できることが知られている。もしこの近似が成立するような短時間フーリエ変換の分析条件を選べば、yk[t]は(2)式のように表される。

以上では時不変な系を考えたが、実際には音源が移動することがあるため、室内伝達系の時不変性の仮定は必ずしも成り立たない。ここで、時変な室内伝達系を仮定することは、(2)式において室内インパルス応答hk[τ]が時刻tごとに変化しうることを許容することに相当する。すなわち、以下の(3)式のように表される。

非定常雑音、時変残響環境下における音声強調問題は、(3)式をもとに観測時間周波数成分Y={yk[t]}1≦k≦K,1≦t≦TからS={sk[t]}1≦k≦K,1≦t≦Tを推定する問題として定式化される。(3)式では室内インパルス応答が時刻ごとに変化することを許容しているため、(2)式に比べて未知パラメータが圧倒的に多くなっており、同じ右辺を与えるsとhとnの組み合わせは無数に存在する。そこで以下では、その無数の組み合わせの中から正しい解を得るためにどのような仮定を置くべきか、について説明する。

<室内伝達系の「半時変性」>

室内インパルス応答が時間的に変化する主要因の1つは話者の移動である。なぜなら、話者の移動に伴い、話者からマイクロホンまでの(壁からの反射音の到来経路を含む)あらゆる音響経路の距離が一斉に変化するからである。

ここで、単一の音響経路のみに着目すると、伝達系のインパルス応答はデルタ関数となるが、経路の距離に変化があった場合、当該インパルス応答のフーリエ変換の位相成分には音響経路の距離変化に伴う到来時間の変化、振幅成分には音響経路の距離変化に伴う強度の変化がそれぞれ反映されることになる。音の強度は伝播距離に対して反比例的であるため、伝播距離がある程度大きい場合には伝播距離に変化があっても振幅成分はほとんど変化しないことになる。

一方、位相成分については、到来時間と周波数に比例した量だけ変化するため、到来時間の変化がたとえわずかであっても、特に高周波域においては著しく変化しうる。

よって、単一の音響経路の伝達系のインパルス応答の短時間フーリエ変換には、話者の移動に伴って振幅は変化しにくく位相は変化しやすい、という性質があることが予想される。

実際の室内インパルス応答はあらゆる音響経路の伝達系のインパルス応答の重ね合わせとなるので、一つ一つのインパルス応答が上述のような傾向があるならば、その重ね合わせもおおよそ同様な傾向があるのではないかと考えられる。

図1は、話者の代わりにマイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分と位相成分を比較したものである。

上記図1の上段は、室内インパルス応答の短時間フーリエ変換の振幅成分(A)及び位相成分(B)を示したものであり、上記図1の下段は、マイクロホンの位置を数センチ変えて測定した室内インパルス応答の短時間フーリエ変換の振幅成分(C)と位相成分(D)を示したものである。なお、室内インパルス応答のデータはRWCP実環境音声・音響データベースのマイクロホンアレーデータベースに収録されている、残響時間0.3[s]の残響可変室(パネル)で測定されたものである。

上記図1から、振幅成分はさほど変化していないのに対し位相成分は著しく変化していることが分かる。

以上より、発明者らは、話者の移動等、環境の軽微な変化については、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変と仮定した特殊な系により扱えるのではないかと考えた。以後、このような系を「半時変系」と呼ぶこととする。

<観測信号の生成モデル>

以上の2つの仮定をベースに、観測信号の生成モデルを定式化する。(3)式に半時変性の仮定を導入すると、観測信号の短時間フーリエ変換は(4)式のように表される。

ただし、φk,t[τ]とHk[t]は室内インパルス応答の短時間フーリエ変換の位相成分と振幅成分であり、Hk[t]=|hk,t[τ]|である。ここで、φk,t[τ]をk,n,mごとに独立に区間[0;2π)上の一様分布に従う確率変数とする。さらに、原音声信号、雑音については、それぞれ平均が0の複素正規分布(5)式及び(6)式に従う確率変数とする。

ただし、Φ(s) k[t]、Φ(n) k[t]はそれぞれ原音声信号および雑音信号の時刻t、k番目の周波数ビンにおけるパワースペクトル密度である。ここで、NC(z;0,λ)は複素正規分布の確率密度関数を表し、(7)式によって与えられる。

証明は省略するが、原音声信号と雑音のガウス性の仮定と、正規分布の再生性より、観測信号の短時間フーリエ変換yk[t]はk、tごとに独立に、(8)式に示す複素正規分布に従うことが示される。

次に、原音声信号のパワースペクトル密度系列Φ(s) k[t]に対し、非負値行列積表現を導入する。基底スペクトルをM個とし、m番目の基底スペクトルをB(s) m,kとする。また、時刻tにおける基底スペクトルB(s) m,kのゲインをG(s) m[t]とすると、原音声信号のパワースペクトル密度系列は(9)式及び(10)式で示される。

なお、基底スペクトルとそのゲインは非負であることに注意する.また、スケールの任意性を除くため、(11)式に示す条件を満たすものとする。

音声の基底スペクトルB(s) m,kは未知変数として観測信号から他の未知パラメータとしてブラインドで推定することもありえるが、ある程度の量のクリーン音声サンプルからあらかじめ学習しておくことを想定する。その意味で、本実施の形態に係る手法(以下、提案法と称する)はセミブラインドな音声強調手法である。基底の事前学習は、クリーン音声データのスペクトログラムに対してNMFを適用することで行うことができる。

<パラメータ推定アルゴリズム>

以上で立てた観測時間周波数成分Y={yk[t]}1≦k≦K,1≦t≦Tの確率密度関数は、未知パラメータの尤度関数に対応する。(8)式より、観測時間周波数成分Yが与えられたもとでの未知パラメータの対数尤度は、定数項・符号を無視すれば、観測スペクトログラム|yk[t]|2と(12)式で示されるパワースペクトル密度系列モデルとの距離を表す板倉斎藤距離(13)式と等しくなる。

従って、提案法は観測スペクトログラム|yk[t]|2とパワースペクトル密度系列モデルΦx k[t]との最適フィッティング問題に帰着する。ここで、雑音信号のパワースペクトル密度系列Φ(n) k[t]に関しても、(14)式のようにモデル化される。

なお、板倉斎藤距離I(H,G(s),B(n),G(n))は、距離尺度を表すβ-divergenceにおいてβを0にした場合に得られる距離に相当する。

(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))の上限関数となる。

すなわち、板倉斎藤距離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)式によって示される。

そして、この上限関数D0の更新則は(17)式〜(23)式で示される。


<システム構成>

次に、原音声信号に雑音及び残響成分が重畳された観測信号を解析して、観測信号に含まれる原音声信号を推定する音響信号解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。

図2に示すように、本発明の実施の形態に係る音響信号解析装置は、CPUと、RAMと、後述する音響信号解析処理ルーチンを実行するためのプログラムを記憶したROMとを備えたコンピュータで構成され、機能的には次に示すように構成されている。

音響信号解析装置100は、入力部10と、演算部20と、記憶部30と、出力部40とを備えている。

入力部10により、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが入力される。記憶部30は、入力部10により入力された観測信号の時系列データを記憶する。また、記憶部30は、後述する各処理での結果を記憶すると共に、本処理ルーチンで用いる各パラメータの初期値を記憶している。

演算部20は、時間周波数解析部21と、初期設定部22と、補助変数更新部23と、パラメータ更新部24と、終了判定部25と、信号変換部26とを備えている。

時間周波数解析部21は、例えばマイクロホンの時系列信号としての観測された観測信号y[t]を入力として、観測信号y[t]の観測時間周波数成分Y={yk[t]}(周波数k=1,・・・,K、時刻t=1,・・・,T)を計算する。また、計算した観測時間周波数成分Yを、記憶部30に記憶しておく。より詳細には、時間周波数解析部21は、例えばマイクロホンで観測された観測信号の時系列データを入力として、短時間フーリエ変換(Short-Time Fourier Transform;STFT)を用いて時間周波数解析を行うことにより観測時間周波数成分Yを計算する。

初期設定部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]の各パラメータの初期値は非負値となるように設定する。

そして、B(s) m,kの値には、一定量のクリーン音声サンプルから予め学習して得た値を設定する。なお、スケールの任意性を排除するため、B(s) m,kの値は上記(11)式を満たすように予め学習されているものとする。

補助変数更新部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に格納する。

また、補助変数更新部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に格納する。

更に、補助変数更新部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に格納する。

パラメータ更新部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に格納する。

また、パラメータ更新部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に格納する。

また、パラメータ更新部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に格納する。

更に、パラメータ更新部24は、(k、τ)の全ての組み合わせの各々について、記憶部30に記憶されているyk[t]、B(s) m,k、G(s) m[t]、ζk,m,t,τ、ηk[t]に基づいて、上記(20)式に従って、上限関数D0を小さくするように、室内インパルス応答のパワーHk[τ]を更新し、記憶部30に格納する。

終了判定部25は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、補助変数更新部23及びパラメータ更新部24の各処理を繰り返す。終了判定部25は、終了条件を満足したと判定した場合には、信号変換部26による処理に移行する。信号変換部26は、記憶部30に記憶されているB(s) m,k 及びG(s) m[t]に基づいて、雑音信号を取り除くことにより推定される原音声信号(以下、推定原音声信号と称する)に変換し、出力部40により、推定原音声信号を出力する。

なお、終了条件としては、繰り返し回数がL-1回目の上限関数(15)式の値と、繰り返し回数がL回目の上限関数(15)式の値との差が、予め定めた閾値よりも小さくなったことを用いればよい。あるいは、終了条件として、繰り返し回数が、予め定められた上限回数に到達したことを用いてもよい。

<音響信号解析装置の作用>

次に、本実施の形態に係る音響信号解析装置100の作用について説明する。まず、解析対象の信号として、原音声信号に雑音及び残響成分が重畳された観測信号の時系列データが音響信号解析装置100に入力され、記憶部30に格納される。そして、音響信号解析装置100において、図3に示す音響信号解析処理ルーチンが実行される。

まず、ステップS101において、記憶部30から、観測信号y[t]を読み込み、当該観測信号y[t]に対して、短時間フーリエ変換を用いた時間周波数分析を行い、観測信号y[t]の観測時間周波数成分Y={yk[t]}(k=1,・・・,K、t=1,・・・,T)を算出すると共に、得られた観測時間周波数成分Yを記憶部30に記憶する。

そして、ステップS102において、乱数を用いて、G(s) m[t]、Hk[τ]、B(n) q,k 、G(n) q[t]の各初期値を設定して、記憶部30に記憶すると共に、一定量のクリーン音声サンプルから予め学習して得た各周波数kのM個の基底スペクトルmからなる基底行列の各要素をB(s) m,kとして設定して、記憶部30に記憶する。

次にステップ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に格納する。

ステップ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に格納する。

ステップ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に格納する。

ステップ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に格納する。

そして、ステップ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に格納する。

次のステップ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より推定原音声信号を出力して音響信号解析処理ルーチンを終了する。

<実施結果>

次に、本実施の形態に係る手法の有用性を示す目的で、残響室内(残響時間1.3秒)での移動音声信号に、PHSの着信音、背景雑音、当該PHSの着信音及び背景雑音を組み合わせた音、粒子落下音、BGM音、他者の音声といった非定常雑音を重畳したものを観測信号として用いた場合の、提案法での評価実験の結果について説明する。

図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は値が大きくなる程、音声の明瞭度が良くなることを示す値である。

上記図4から分かるように、移動音声信号にBGM又は他者の音声を重畳した場合のMFCCDの評価結果を除いて、観測信号に含まれる音声の明瞭度に比べ提案法による推定原音声信号の明瞭度の方が向上するという結果を得られた。特に、移動音声信号にPHSの着信音を重畳した場合のSIRに関する評価結果では、4倍以上の改善が見られた。

以上説明したように、本発明の実施の形態に係る音響信号解析装置によれば、上記(15)式で示される上限関数の値が小さくなるように、係数行列の各要素G(s) m[t]と、基底行列の各要素B(n) q,kと、係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、補助変数ζk,m,t,τと、補助変数ξk,q,tと、補助変数ηk[t]とを繰り返し更新することにより、音声信号に雑音及び残響成分が重畳した観測信号から精度よく音声信号を推定して、音声信号の明瞭度を向上させることができる。

このように、本発明で提案する手法では、音源の移動等に関して、室内インパルス応答の短時間フーリエ変換の振幅成分を時不変、位相成分を時変とする観測信号の生成モデルを立てて、観測信号に含まれる原音声信号を高い精度で推定する

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上述の音響信号解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、WWWシステムを利用している場合であれば、ホームページ提供環境(あるいは表示環境)も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

また、本実施の形態に係る手法では、パワースペクトル密度系列モデルΦx k[t]と観測時間周波数成分yk[t]との距離の尺度として板倉斎藤距離I(H,G(s),B(n),G(n))を用いたが、これに限らず、β-divergenceの一種であるGeneralized Kullback-Leibler距離やFrobeniusノルムといった尺度を用いるようにしてもよい。

10 入力部
20 演算部
21 時間周波数解析部
22 初期設定部
23 補助変数更新部
24 パラメータ更新部
25 終了判定部
26 信号変換部
30 記憶部
40 出力部
100 音響信号解析装置

Claims (5)

  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]とに基づいて算出されるパワースペクトル密度系列モデルΦx k[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のパワースペクトル密度系列モデルΦx k[t]とに基づいて、前記係数行列の各要素G(s) m[t]と、前記基底行列の各要素B(n) q,kと、前記係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新するパラメータ更新部と、
    予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う終了判定部と、
    を含む音響信号解析装置。
  2. 前記パワースペクトル密度系列モデルΦx k[t]と、前記観測時間周波数成分yk[t]との距離の尺度としてβ−divergenceを用いた
    請求項1記載の音響信号解析装置。
  3. 補助変数更新部を更に含み、
    前記パワースペクトル密度系列モデルΦx k[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のパワースペクトル密度系列モデルΦx k[t]と、(k、m、t、τ)の全ての組み合わせについての補助変数ζk,m,t,τと、(k、q、t)の全ての組み合わせについての補助変数ξk,q,tと、各時刻t及び各周波数kについての補助変数ηk[t]を用いて表され、かつ、前記パワースペクトル密度系列モデルΦx k[t]と、前記観測時間周波数成分yk[t]とのβ−divergenceの上限関数である補助関数とし、
    前記補助変数更新部は、前記補助関数を小さくするように、前記基底行列の各要素B(s) m,kと、前記係数行列の各要素G(s) m[t]と、前記基底行列の各要素B(n) q,kと、前記係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]と、各時刻t及び各周波数kのパワースペクトル密度系列モデルΦx k[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]とに基づいて算出されるパワースペクトル密度系列モデルΦx k[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のパワースペクトル密度系列モデルΦx k[t]とに基づいて、前記係数行列の各要素G(s) m[t]と、前記基底行列の各要素B(n) q,kと、前記係数行列の各要素G(n) q[t]と、各時刻τの室内インパルス応答の各周波数kのパワーHk[τ]とを更新し、
    終了判定部によって、予め定められた終了条件を満たすまで、前記パラメータ更新部による更新を繰り返し行う
    音響信号解析方法。
  5. コンピュータを、請求項1〜請求項3の何れか1項に記載の音響信号解析装置の各部として機能させるためのプログラム。
JP2013181701A 2013-09-02 2013-09-02 音響信号解析装置、方法、及びプログラム Active JP6142402B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2013181701A JP6142402B2 (ja) 2013-09-02 2013-09-02 音響信号解析装置、方法、及びプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2013181701A JP6142402B2 (ja) 2013-09-02 2013-09-02 音響信号解析装置、方法、及びプログラム

Publications (2)

Publication Number Publication Date
JP2015049406A JP2015049406A (ja) 2015-03-16
JP6142402B2 true JP6142402B2 (ja) 2017-06-07

Family

ID=52699461

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2013181701A Active JP6142402B2 (ja) 2013-09-02 2013-09-02 音響信号解析装置、方法、及びプログラム

Country Status (1)

Country Link
JP (1) JP6142402B2 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6448567B2 (ja) * 2016-02-23 2019-01-09 日本電信電話株式会社 音響信号解析装置、音響信号解析方法、及びプログラム

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8015003B2 (en) * 2007-11-19 2011-09-06 Mitsubishi Electric Research Laboratories, Inc. Denoising acoustic signals using constrained non-negative matrix factorization
JP2012027196A (ja) * 2010-07-22 2012-02-09 Nippon Telegr & Teleph Corp <Ntt> 信号分析装置、方法、及びプログラム

Also Published As

Publication number Publication date
JP2015049406A (ja) 2015-03-16

Similar Documents

Publication Publication Date Title
Aarabi et al. Phase-based dual-microphone robust speech enhancement
Adler et al. Audio inpainting
Vincent et al. From blind to guided audio source separation: How models and side information can improve the separation of sound
Goh et al. Kalman-filtering speech enhancement method based on a voiced-unvoiced speech model
KR100486736B1 (ko) 두개의 센서를 이용한 목적원별 신호 분리방법 및 장치
AU2009278263B2 (en) Apparatus and method for processing an audio signal for speech enhancement using a feature extraction
Xu et al. An experimental study on speech enhancement based on deep neural networks
Matassini et al. Optimizing of recurrence plots for noise reduction
JP5587396B2 (ja) 信号分離のためのシステム、方法、および装置
Yegnanarayana et al. Enhancement of reverberant speech using LP residual signal
CN1168069C (zh) 识别系统和识别方法
US9111526B2 (en) Systems, method, apparatus, and computer-readable media for decomposition of a multichannel music signal
CN104053092B (zh) 用于双麦克风通信装置的降噪
Magi et al. Stabilised weighted linear prediction
Shrawankar et al. Techniques for feature extraction in speech recognition system: A comparative study
JP2009128906A (ja) 音響信号と雑音信号とを含む混成信号の雑音を除去するための方法およびシステム
JP2004502977A (ja) サブバンド指数平滑雑音消去システム
JP2005249816A (ja) 信号強調装置、方法及びプログラム、並びに音声認識装置、方法及びプログラム
US8548803B2 (en) System and method of processing a sound signal including transforming the sound signal into a frequency-chirp domain
EP1356461A1 (fr) Procede et dispositif de reduction de bruit
JP4469882B2 (ja) 音響信号処理方法及び装置
Ephraim et al. On second-order statistics and linear estimation of cepstral coefficients
Xiao et al. Deep beamforming networks for multi-channel speech recognition
JP4880036B2 (ja) 音源と室内音響の確率モデルに基づく音声残響除去のための方法及び装置
JP5227393B2 (ja) 残響除去装置、残響除去方法、残響除去プログラム、および記録媒体

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20150709

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20150709

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20160706

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20160809

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20170328

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20170418

R150 Certificate of patent or registration of utility model

Ref document number: 6142402

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150