(57)【要約】【課題】一部欠損した観測スペクトログラムから時間領域信号を精度良く推定することができるようにする。【解決手段】信号推定部５０が、欠損していない部分における、観測スペクトログラムと、時間領域信号から求められるスペクトログラムとの誤差と、時間領域信号の基底スペクトル、及びアクティベーションパラメータの積と、時間領域信号から求められるスペクトログラムの絶対値との誤差、欠損していない部分における、基底スペクトル、及びアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び時間領域信号から抽出される各時刻のケプストラム特徴量の事前確率とを用いて表される規準を最適化するように、時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する。【選択図】図２
2018136419.tif 000002
【特許請求の範囲】
【請求項１】
一部が欠損した、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを入力として、時間領域信号を推定する信号推定装置であって、
欠損していない部分における、前記観測スペクトログラムと、前記時間領域信号から求められるスペクトログラムとの誤差と、
基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、
欠損していない部分における、基底スペクトル、及び各時刻におけるアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び
基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率の少なくとも一つとを用いて表される規準を最適化するように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する信号推定部
を含む信号推定装置。
【請求項２】
前記基底スペクトル、及び前記アクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、又は前記基底スペクトル、及び前記アクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差は、ユークリッド距離又はＫＬダイバージェンスで表される請求項１記載の信号推定装置。
【請求項３】
前記ケプストラム特徴量の事前確率は、予め定められたケプストラム特徴量の混合正規分布における、基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の確率である請求項１又は２記載の信号推定装置。
【請求項４】
前記信号推定部は、
前記規準の上界関数である補助関数を小さくするように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを更新するパラメータ更新部と、
予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させる収束判定部と、
を含む請求項１〜請求項３の何れか１項記載の信号推定装置。
【請求項５】
一部が欠損した、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを入力として、時間領域信号を推定する信号推定装置における信号推定方法であって、
信号推定部が、欠損していない部分における、前記観測スペクトログラムと、前記時間領域信号から求められるスペクトログラムとの誤差と、
基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、
欠損していない部分における、基底スペクトル、及び各時刻におけるアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び
基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率の少なくとも一つとを用いて表される規準を最適化するように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する
信号推定方法。
【請求項６】
前記基底スペクトル、及び前記アクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、又は前記基底スペクトル、及び前記アクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差は、ユークリッド距離又はＫＬダイバージェンスで表される請求項５記載の信号推定方法。
【請求項７】
前記ケプストラム特徴量の事前確率は、予め定められたケプストラム特徴量の混合正規分布における、基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の確率である請求項５又は６記載の信号推定方法。
【請求項８】
請求項１〜請求項４の何れか１項に記載の信号推定装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号推定装置、方法、及びプログラムに係り、特に、観測スペクトログラムから、時間領域信号を推定する信号推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音源分離は音声信号や音響信号が重畳した混合音から所望の目的音を抽出し分離する技術であり、雑音下での音声認識精度向上や他の音響信号処理の前処理として利用され、これまでに様々な音源分離手法が提案されてきた。
【０００３】
音源分離の有効な手法の一つに時間周波数マスクに基づく手法がある。この手法では、混合音のスペクトログラムの各時間周波数点においてどの音源が最も優勢であるかを識別することで、同一音源に属すると見なされた時間周波数点の成分のみを通過させるような時間周波数マスクを構成することが目的となる。近年深層ニューラルネットワークに基づくアプローチにより高い信号対歪比での音源分離が可能となってきているが、この手法により得られる各音源のスペクトログラムは他音源が優勢であると見なされた時間周波数点の成分が欠損したものになるため、音声認識など後段に用いられる様々な音響信号処理の性能劣化を招く。従って、このように欠損があるスペクトログラムからいかにして欠損した成分を復元するかが重要課題となる。
【０００４】
スペクトログラムの欠損成分を復元する手法として非負値行列因子分解(Non-negative Matrix Factorization; NMF) に基づく方法が提案されている（非特許文献１）。この方法は、非負値行列と見なしたスペクトログラムを、二つの非負値行列の積で非欠損箇所において近似することにより欠損箇所の成分を補完しようというものである。これは、スペクトログラムが低ランクな行列となっているという仮定を手がかりとしていることに相当する。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Paris Smaragdis, Bhiksha Raj, and Madhusudana Shashanka, “Missing data imputation for time-frequency representations of audio signals," Journal of Signal Processing Systems, Volume 65, Issue 3, pp. 361-370, 2011.
【非特許文献２】李莉, 亀岡弘和, 樋口卓哉, 猿渡洋, “ケプストラム距離正則化半教師ありNMF による音声強調," 日本音響学会2016 年春季研究発表会講演論文集, 1-P-27, pp. 721-724, Mar. 2016.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上記の従来技術の方法では欠損箇所における振幅成分しか推定されないため、信号を構成するためには別手法で位相成分も推定する必要がある。
【０００７】
短時間Fourier 変換やウェーブレット変換（定Q 変換）などに代表される時間周波数表現は時間領域信号の冗長表現であるため、各時間周波数点の振幅成分および位相成分は（局所的に）互いに依存関係にある。これは、非欠損箇所の振幅成分だけでなく位相成分にも欠損箇所の振幅・位相成分を推定する手がかりがあることを示唆している。また、音響特徴量の事前情報が得られる場合にはその統計分布も欠損成分を補完するための手がかりとなりうる。そこで、従来考慮されていなかったこれらの手がかりを用いることが考えられる。
【０００８】
本発明は、上記事情を鑑みてなされたものであり、一部欠損した観測スペクトログラムから時間領域信号を精度良く推定することができる信号推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記の目的を達成するために本発明に係る信号推定装置は、一部が欠損した、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを入力として、時間領域信号を推定する信号推定装置であって、欠損していない部分における、前記観測スペクトログラムと、前記時間領域信号から求められるスペクトログラムとの誤差と、基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、欠損していない部分における、基底スペクトル、及び各時刻におけるアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率の少なくとも一つとを用いて表される規準を最適化するように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する信号推定部を含んで構成されている。
【００１０】
本発明に係る信号推定方法は、一部が欠損した、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを入力として、時間領域信号を推定する信号推定装置における信号推定方法であって、信号推定部が、欠損していない部分における、前記観測スペクトログラムと、前記時間領域信号から求められるスペクトログラムとの誤差と、基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、欠損していない部分における、基底スペクトル、及び各時刻におけるアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率の少なくとも一つとを用いて表される規準を最適化するように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する。
【００１１】
本発明に係るプログラムは、上記の信号推定装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１２】
以上説明したように、本発明の信号推定装置、方法、及びプログラムによれば、欠損していない部分における、前記観測スペクトログラムと、前記時間領域信号から求められるスペクトログラムとの誤差と、前記時間領域信号の基底スペクトル、及びアクティベーションパラメータの積と、前記時間領域信号から求められるスペクトログラムの絶対値との誤差、欠損していない部分における、基底スペクトル、及びアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率の少なくとも一つとを用いて表される規準を最適化するように、前記時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定することにより、一部欠損した観測スペクトログラムから時間領域信号を精度良く推定することができる、という効果が得られる。
【図面の簡単な説明】
【００１３】
【図１】巡回行列となるように選択する方法を説明するための図である。
【図２】本発明の実施の形態に係る信号推定装置の構成を示す概略図である。
【図３】本発明の実施の形態に係る信号推定装置の学習部の構成を示す概略図である。
【図４】本発明の実施の形態に係る信号推定装置の信号推定部の構成を示す概略図である。
【図５】本発明の実施の形態に係る信号推定装置における信号推定処理ルーチンの内容を示すフローチャートである。
【図６】実験結果を示す図である。
【発明を実施するための形態】
【００１４】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音源分離などにより一部欠損したスペクトログラムから時間領域信号を構成する技術である。
【００１５】
＜本発明の実施の形態の概要＞
本発明の実施の形態は、上記非特許文献１で用いられる振幅スペクトログラムの低ランク構造の仮定とともに短時間Fourier変換やウェーブレット変換（定Q変換）などの時間周波数表現が時間領域信号の冗長表現であることより導かれる各時間周波数点の振幅成分および位相成分の局所的な制約条件、及び目的音源のケプストラム特徴量の事前情報（統計分布）のいずれかまたは両方を手がかりとして一部の時間周波数成分が欠損したスペクトログラムから直接時間領域信号を構成する。
【００１６】
＜問題の定式化＞
一部の時間周波数点の成分が欠損した複素スペクトログラムを
2018136419.tif 000003
とし、その各要素をYk,mと表す。ただし、k∈{1,...,K}及びm∈{1,...,M}はそれぞれ周波数ビン、時刻を表すインデックスである。非欠損箇所の集合を
2018136419.tif 000004
とすると、
2018136419.tif 000005
は欠損箇所において
【００１７】
【数１】
2018136419.tif 000006
【００１８】
となっているものとする。この複素スペクトログラムに対応する時間領域信号
2018136419.tif 000007
（全サンプル点数がN）を推定することが目的で、
2018136419.tif 000008
の時間周波数表現（短時間Fourier 変換やウェーブレット変換など）は一般に
2018136419.tif 000009
と表される。ただし
2018136419.tif 000010
は時刻tmに局在する中心周波数がωｋの窓掛けされた複素正弦波である。上記非特許文献１では、振幅スペクトログラムを二つの非負値行列の積（低ランクな非負値行列）で近似できると仮定し、｜Yｋ，ｍ｜を非欠損箇所
2018136419.tif 000011
において
【００１９】
【数２】
2018136419.tif 000012
【００２０】
で近似することで、Xk,mを欠損成分の推定値とする方法が提案されている。
【００２１】
振幅スペクトログラムを式(2)で表現することを非負値行列因子分解(Non-negative Matrix Fac-torization, NMF) といい、各時刻の振幅スペクトルXk,mをL個の基底スペクトルの非負係数U1，m,...,UL,mによる線形和で表現することに相当する。
【００２２】
NMFに基づく上記の手法は、対象とする振幅スペクトログラム全体が低ランクな行列で近似できるという、振幅スペクトログラムの大域的な構造に関する仮定を手がかりに欠損成分を補完しようとするものである。一方で、
2018136419.tif 000013
の複素スペクトログラムは
2018136419.tif 000014
と表され、N < KMの場合には
2018136419.tif 000015
の冗長表現になることから、非欠損箇所のYk,mの振幅成分だけでなく位相成分も周辺の欠損箇所の振幅・位相成分を推定する手がかりとすることができる。また、対象スペクトログラムのケプストラム特徴量の事前情報が得られる場合にはその統計分布も欠損成分を補完するための手がかりとなりうる。以上の3種類の手がかりを用いて、一部欠損したスペクトログラムから時間領域信号を直接構成する問題は
【００２３】
【数３】
2018136419.tif 000016
【００２４】
を最小化する最適化問題として定式化することができる。ただし、λ1、λ2、λ3は正則化パラメータ、
2018136419.tif 000017
は推定したい未知パラメータの集合であり、
2018136419.tif 000018
である。また、
2018136419.tif 000019
は観測スペクトログラムで欠損していない時間周波数点の集合（非欠損領域）である。
【００２５】
2018136419.tif 000020
は誤差の大きさを表す関数であり、ここではEuclid 距離（二乗誤差）
【００２６】
【数４】
2018136419.tif 000021
【００２７】
またはKullback-Leibler (KL) ダイバージェンス
【００２８】
【数５】
2018136419.tif 000022
【００２９】
を用いる。式(3) の第一項は観測複素スペクトログラムと
2018136419.tif 000023
の複素スペクトログラムの非欠損領域における二乗誤差、第二項は式(2) と
2018136419.tif 000024
の振幅スペクトログラムの誤差の大きさ、第三項は式(2) と観測振幅スペクトログラムの誤差の大きさを表す。第四項は以下で与えられるケプストラム距離規準の正則化項である。
【００３０】
【数６】
2018136419.tif 000025
2018136419.tif 000026
【００３１】
ただし
2018136419.tif 000027
はメル周波数ケプストラム係数（Mel-Frequency Cepstrum Coefficients: MFCC）であり、fr,kはr番目のメルフィルタバンク係数、
2018136419.tif 000028
は離散コサイン変換係数である。式(6) はパラメータ
2018136419.tif 000029
の混合正規分布から生成される確率の対数を表す。ただし、
2018136419.tif 000030
はp番目の正規分布の重み、平均及び分散を表す。
【００３２】
＜パラメータ推定アルゴリズム＞
2018136419.tif 000031
を最小化するパラメータ
2018136419.tif 000032
を解析的に得ることは難しいが、補助関数法と呼ぶ原理に基づき、
2018136419.tif 000033
の停留点への収束が保証された反復更新アルゴリズムを導くことができる。
【００３３】
＜補助関数法＞
F(θ)をθに関して最小化したい目的関数とすると、
2018136419.tif 000034
を満たす関数F+(θ,α)を補助関数、αを補助変数と呼ぶ。このような補助関数を設計できれば、
2018136419.tif 000035
と
2018136419.tif 000036
を交互に繰り返すことで、目的関数F(θ)の停留点を得ることができる。この最適化手法を補助関数法と呼ぶ。以下では各パラメータに対する更新式を得るための式(3)の補助関数を設計し、更新式を導出する。
【００３４】
＜ｓの更新式＞
式(3) の中で
2018136419.tif 000037
に関係する項は第一項と第二項である。まず、第二項の補助関数を設計する。
【００３５】
2018136419.tif 000038
がEuclid 距離（二乗誤差）のとき式(3)の第二項は
【００３６】
【数７】
2018136419.tif 000039
【００３７】
となる。この関数は
2018136419.tif 000040
に関して微分不可能な点をもつが、以下の補題を用いることで微分可能な補助関数を設計することができる。
【００３８】
（補題1．）複素数z及び｜ａ｜＝１を満たす複素数aに対して、以下が成立する。
【００３９】
【数８】
2018136419.tif 000041
【００４０】
ただし等号成立はa=z/|ｚ|である。
【００４１】
補題1より
【００４２】
【数９】
2018136419.tif 000042
【００４３】
が成り立つ。この不等式の右辺は補助変数
2018136419.tif 000043
をもつ
【００４４】
2018136419.tif 000044
の補助関数であり、等号は
2018136419.tif 000045
のとき成立する。一方、
2018136419.tif 000046
がKLダイバージェンスのとき式(3)の第二項は
【００４５】
【数１０】
2018136419.tif 000047
【００４６】
となる。この関数も
2018136419.tif 000048
に関して微分不可能な点をもつが、以下に従い微分可能な補助関数を設計することができる。まず、対数関数は凹関数であり、凹関数は任意の点の接線より大きくなることはないため、
【００４７】
【数１１】
2018136419.tif 000049
【００４８】
が成り立つ。ただし、
【００４９】
【数１２】
2018136419.tif 000050
【００５０】
である。また、等号は
2018136419.tif 000051
のとき成立する。次に、
2018136419.tif 000052
に関して
【００５１】
【数１３】
2018136419.tif 000053
2018136419.tif 000054
【００５２】
が成り立つことを利用し、Dk,mの符号に応じていずれかを用いることで
【００５３】
【数１４】
2018136419.tif 000055
【００５４】
が言える。ただし、Fk,m、Gk,mは
【００５５】
【数１５】
2018136419.tif 000056
2018136419.tif 000057
【００５６】
である。式(16)の不等式の右辺は補助変数
2018136419.tif 000058
をもつ
2018136419.tif 000059
の補助関数であり、等号は
2018136419.tif 000060
のとき成立する。
【００５７】
いずれの補助関数も
2018136419.tif 000061
に関して微分可能かつ凸なので、勾配法などを用いて大域最適となる点を探索することができる。
【００５８】
2018136419.tif 000062
がEuclid距離の場合の式(3)の補助関数は
【００５９】
【数１６】
2018136419.tif 000063
【００６０】
となる。よって、
2018136419.tif 000064
に関する偏微分は
【００６１】
【数１７】
2018136419.tif 000065
【００６２】
となる。ただし、rk,mは
【００６３】
【数１８】
2018136419.tif 000066
【００６４】
である。一方、
2018136419.tif 000067
がKLダイバージェンスの場合の式(3)の補助関数は
【００６５】
【数１９】
2018136419.tif 000068
【００６６】
となる。よって、
2018136419.tif 000069
に関する偏微分は
【００６７】
【数２０】
2018136419.tif 000070
【００６８】
となる。なお、
2018136419.tif 000071
の形になっている項は高速Fourier変換(Fast Fourier Trans-form, FFT) を用いて効率的に計算することができる。
【００６９】
また、
2018136419.tif 000072
がEuclid 距離の場合は
2018136419.tif 000073
が巡回行列となることを利用してさらに以下の補助関数を設計することでアルゴリズムの効率化を図ることができる。今、
2018136419.tif 000074
を
【００７０】
【数２１】
2018136419.tif 000075
【００７１】
と定義すると、
2018136419.tif 000076
より
【００７２】
【数２２】
2018136419.tif 000077
2018136419.tif 000078
【００７３】
が成り立つ。よって、この不等式の右辺は式(3)の第一項の補助関数として用いることができる（補助変数は
2018136419.tif 000079
）。等号は
2018136419.tif 000080
のとき成立する。従って、式(10) と合わせて
【００７４】
【数２３】
2018136419.tif 000081
【００７５】
を
2018136419.tif 000082
の補助関数とすることができる。
【００７６】
2018136419.tif 000083
の
2018136419.tif 000084
に関する偏微分
【００７７】
【数２４】
2018136419.tif 000085
【００７８】
を0と置くことにより、大域最適解
【００７９】
【数２５】
2018136419.tif 000086
【００８０】
を得る。式(28) は
2018136419.tif 000087
の逆行列を含むが、
2018136419.tif 000088
を
2018136419.tif 000089
が巡回行列となるように選ぶことができる（図１）ため、その場合は離散Fourier変換行列
2018136419.tif 000090
で
2018136419.tif 000091
のように対角化することが可能である。これにより、
2018136419.tif 000092
となり、
2018136419.tif 000093
は対角行列より効率的に計算することができる。なお、
2018136419.tif 000094
を
2018136419.tif 000095
が単位行列となるように
2018136419.tif 000096
を選ぶこともできる。
【００８１】
例えば、
2018136419.tif 000097
がHanning関数の平方根を窓関数とするSTFTに相当するように
2018136419.tif 000098
を選べば
2018136419.tif 000099
は単位行列となることが示せる。
【００８２】
以上より、
2018136419.tif 000100
の更新方法は以下のようにまとめられる。
【００８３】
2018136419.tif 000101
がEuclid距離の場合、以下の方法１又は方法２により、更新する。
【００８４】
（方法1）
1.
2018136419.tif 000102
を以下に従い更新する。
【００８５】
【数２６】
2018136419.tif 000103
【００８６】
2. 式(20) を用いて
2018136419.tif 000104
を勾配法により更新する。
【００８７】
（方法2：
2018136419.tif 000105
が巡回行列の場合）
1.
2018136419.tif 000106
を以下に従い更新する。
【００８８】
【数２７】
2018136419.tif 000107
2018136419.tif 000108
【００８９】
2.
2018136419.tif 000109
を式(28) に従い更新する。
【００９０】
2018136419.tif 000110
がKLダイバージェンスの場合、以下の方法１により、更新する。
【００９１】
（方法1）
1.
2018136419.tif 000111
を以下に従い更新する。
【００９２】
【数２８】
2018136419.tif 000112
2018136419.tif 000113
2018136419.tif 000114
【００９３】
2. 式(23) を用いて
2018136419.tif 000115
を勾配法により更新する。
【００９４】
＜ＨとＵの更新式＞
式(3) の中で
2018136419.tif 000116
と
2018136419.tif 000117
に関係する項は第二、第三、第四項であるが、
2018136419.tif 000118
と
2018136419.tif 000119
は式(3) の第三項を規準として
2018136419.tif 000120
のみから推定することもできる。以下では第三項だけでなく第二、第四項を含めた場合の
2018136419.tif 000121
と
2018136419.tif 000122
の更新方法について述べる。
【００９５】
まずλ3=0の場合（第二項と第三項のみ）を考える。第二、第三項の補助関数は従来のNMF と同様Jensenの不等式を用いて導くことができる。
【００９６】
2018136419.tif 000123
がEuclid距離（二乗誤差）のとき、二次関数は凸関数のためJensenの不等式
【００９７】
【数２９】
2018136419.tif 000124
2018136419.tif 000125
【００９８】
が成り立つ。ただし、βk,l,mは
2018136419.tif 000126
を満たす非負の変数であり、等号は
2018136419.tif 000127
2018136419.tif 000128
のとき成立する。よって、式(35)の右辺は式(3) の第二項の補助関数として用いることができる。このとき、
2018136419.tif 000129
が補助変数となる。また、式(3)の第三項についても同様に
【００９９】
【数３０】
2018136419.tif 000130
【０１００】
が成り立ち、右辺を補助関数とすることができる。式(35)、(36)より、式(3)の補助関数は
【０１０１】
【数３１】
2018136419.tif 000131
【０１０２】
と書ける。
【０１０３】
2018136419.tif 000132
に関する偏微分
【０１０４】
【数３２】
2018136419.tif 000133
2018136419.tif 000134
【０１０５】
を0と置くことにより、
2018136419.tif 000135
を最小とするHk,l,Ul,m
【０１０６】
【数３３】
2018136419.tif 000136
2018136419.tif 000137
【０１０７】
を得る。
【０１０８】
2018136419.tif 000138
がKLダイバージェンスのとき、負の対数関数は凸関数のためJensenの不等式
【０１０９】
【数３４】
2018136419.tif 000139
2018136419.tif 000140
【０１１０】
が成り立つ。
【０１１１】
ただし、γk,l,mは
2018136419.tif 000141
を満たす非負の変数であり、等号は
2018136419.tif 000142
2018136419.tif 000143
のとき成立する。よって、式(42)の右辺は式(3)の第二項の補助関数として用いることができる。このとき、
2018136419.tif 000144
が補助変数となる。また、式(3)の第三項についても同様に
【０１１２】
【数３５】
2018136419.tif 000145
【０１１３】
が成り立ち、右辺を補助関数とすることができる。式(42)、(43)より、式(3)の補助関数は
【０１１４】
【数３６】
2018136419.tif 000146
【０１１５】
と書ける。
【０１１６】
2018136419.tif 000147
のHk,l、Ul,mに関する偏微分
【０１１７】
【数３７】
2018136419.tif 000148
2018136419.tif 000149
【０１１８】
を0 と置くことにより
2018136419.tif 000150
を最小とする
2018136419.tif 000151
【０１１９】
【数３８】
2018136419.tif 000152
2018136419.tif 000153
【０１２０】
を得る。
【０１２１】
次にλ３≠0の場合を考える。上記非特許文献２より
【０１２２】
【数３９】
2018136419.tif 000154
【０１２３】
が成り立つため、右辺を式(3) の第四項の補助関数とすることができる。ただし、
【０１２４】
【数４０】
2018136419.tif 000155
2018136419.tif 000156
2018136419.tif 000157
【０１２５】
である。また、
2018136419.tif 000158
は
【０１２６】
【数４１】
2018136419.tif 000159
2018136419.tif 000160
【０１２７】
で定義され、δxは条件xを満たす場合に1、それ以外の場合に0となる指示関数である。式(49)の等号は
【０１２８】
【数４２】
2018136419.tif 000161
2018136419.tif 000162
2018136419.tif 000163
2018136419.tif 000164
【０１２９】
のとき成立する。ただしωp,p,r,mは
2018136419.tif 000165
を満たす任意の正の定数である。
【０１３０】
以上より、式(42)、(43) および式(49) より、各不等式の右辺の和
【０１３１】
【数４３】
2018136419.tif 000166
2018136419.tif 000167
【０１３２】
は、HとUの更新式を得るための補助関数として用いることができる。ただし、
2018136419.tif 000168
はすべての補助変数の集合とする。式(59) を
2018136419.tif 000169
に関して偏微分し、0と置くことにより
2018136419.tif 000170
と
2018136419.tif 000171
の更新式
【０１３３】
【数４４】
2018136419.tif 000172
2018136419.tif 000173
【０１３４】
を得る。ただし、
2018136419.tif 000174
は以下で定義される。
【０１３５】
【数４５】
2018136419.tif 000175
2018136419.tif 000176
2018136419.tif 000177
2018136419.tif 000178
2018136419.tif 000179
2018136419.tif 000180
【０１３６】
以上より、
2018136419.tif 000181
の更新方法は以下のようにまとめられる。
【０１３７】
（方法1）
2018136419.tif 000182
がEuclid 距離の場合
【０１３８】
1.
2018136419.tif 000183
を以下に従い更新する。
【０１３９】
【数４６】
2018136419.tif 000184
【０１４０】
2. 式(40)、 (41) にλ1= 0 を代入した式により
2018136419.tif 000185
を更新する。
【０１４１】
2018136419.tif 000186
がKLダイバージェンスの場合
【０１４２】
1.
2018136419.tif 000187
を以下に従い更新する。
【０１４３】
【数４７】
2018136419.tif 000188
【０１４４】
2. 式(47)、(48)にλ1=0を代入した式により
2018136419.tif 000189
を更新する。
【０１４５】
（方法2）
2018136419.tif 000190
がEuclid 距離の場合
【０１４６】
1.
2018136419.tif 000191
を以下に従い更新する。
【０１４７】
【数４８】
2018136419.tif 000192
【０１４８】
2. 式(40)、 (41) により
2018136419.tif 000193
を更新する。
【０１４９】
2018136419.tif 000194
がKLダイバージェンスの場合
【０１５０】
1.
2018136419.tif 000195
を以下に従い更新する。
【０１５１】
【数４９】
2018136419.tif 000196
【０１５２】
2. 式(47)、 (48) により
2018136419.tif 000197
を更新する。
【０１５３】
（方法3）
1.
2018136419.tif 000198
を以下に従い更新する。
【０１５４】
【数５０】
2018136419.tif 000199
2018136419.tif 000200
2018136419.tif 000201
2018136419.tif 000202
2018136419.tif 000203
2018136419.tif 000204
2018136419.tif 000205
【０１５５】
2. 式(60)、 (61) により
2018136419.tif 000206
を更新する。
【０１５６】
＜システム構成＞
【０１５７】
次に、一部欠損したスペクトログラムから時間領域信号を推定する信号推定装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【０１５８】
図２に示すように、本発明の実施の形態に係る信号推定装置は、ＣＰＵと、ＲＡＭと、後述する信号推定処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【０１５９】
図２に示すように、信号推定装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。
【０１６０】
入力部１０は、学習サンプルの目的音響信号の時系列データを受け付ける。また、入力部１０は、目的音響信号について、一部欠損した観測スペクトログラムを受け付ける。
【０１６１】
演算部２０は、学習部３０と、パラメータ記憶部４０と、信号推定部５０とを備えている。
【０１６２】
図３に示すように、学習部３０は、特徴量抽出部３２と、パラメータ学習部３６とを備えている。
【０１６３】
特徴量抽出部３２は、入力部１０によって受け付けた学習サンプルの目的音響信号の時系列データから、目的音響信号のメル周波数ケプストラム係数（Mel-Frequency Cepstrum Coefficients: MFCC）を抽出する。
【０１６４】
パラメータ学習部３６は、特徴量抽出部３２によって抽出された目的音響信号のメル周波数ケプストラム係数に基づいて、メル周波数ケプストラム係数のパラメータ
2018136419.tif 000207
を学習し、パラメータ記憶部４０に格納する。
【０１６５】
信号推定部５０は、目的音響信号についての、一部欠損し、かつ、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを入力として、パラメータ記憶部４０に格納された、メル周波数ケプストラム係数のパラメータを用いて、欠損していない部分における、観測スペクトログラムと、時間領域信号から求められるスペクトログラムとの誤差と、時間領域信号の各周波数における振幅スペクトルを表わす基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、時間領域信号から求められるスペクトログラムの絶対値との誤差と、欠損していない部分における、基底スペクトル、及び各時刻における音量を表すアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差と、基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率を用いて表される、上記（３）式に示す規準を最適化するように、時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定する。
【０１６６】
本実施の形態では、図４に示すように、信号推定部５０は、初期値設定部５２、補助変数更新部５４、パラメータ更新部５８、及び収束判定部６０を備えている。
【０１６７】
初期値設定部５２は、基底スペクトル、及び各時刻におけるアクティベーションパラメータの各々の初期値を設定する。例えば、上記式（２）における第三項から、
2018136419.tif 000208
を推定し、これを固定したのち、上記式（２）における第一項及び第二項から時間領域信号
2018136419.tif 000209
を推定して初期値とする。
【０１６８】
補助変数更新部５４は、時間領域信号
2018136419.tif 000210
に基づいて、上記式（２９）又は上記式（３０）、（３１）に従って、補助変数
2018136419.tif 000211
又は
2018136419.tif 000212
を更新する。あるいは、上記式（３２）〜（３４）に従って、補助変数
2018136419.tif 000213
を更新する。
【０１６９】
また、補助変数更新部５４は、
2018136419.tif 000214
に基づいて、上記式（６８）又は上記式（６９）に従って、補助変数
2018136419.tif 000215
又は
2018136419.tif 000216
を更新する。あるいは、上記式（７０）又は上記式（７１）に従って、補助変数
2018136419.tif 000217
又は
2018136419.tif 000218
を更新する。あるいは、
2018136419.tif 000219
と、パラメータ記憶部４０に格納されているメル周波数ケプストラム係数のパラメータ
2018136419.tif 000220
と、時間領域信号
2018136419.tif 000221
から抽出されるメル周波数ケプストラム係数
2018136419.tif 000222
とに基づいて、上記式（７２）〜式（７８）に従って、補助変数
2018136419.tif 000223
を更新する。
【０１７０】
パラメータ更新部５８は、上記式（３）に示す規準の上界関数である、上記式（１９）又は（２２）に示す補助関数を小さくするように、時間領域信号
2018136419.tif 000224
、基底スペクトル
2018136419.tif 000225
、及び各時刻におけるアクティベーションパラメータ
2018136419.tif 000226
を更新する。
【０１７１】
具体的には、パラメータ更新部５８は、補助変数
2018136419.tif 000227
と、観測スペクトログラム
2018136419.tif 000228
と、
2018136419.tif 000229
とに基づいて、上記式（２０）を用いて、勾配法により、時間領域信号
2018136419.tif 000230
を更新する。あるいは、補助変数
2018136419.tif 000231
と、観測スペクトログラム
2018136419.tif 000232
と、
2018136419.tif 000233
とに基づいて、上記式（２８）に従って、時間領域信号
2018136419.tif 000234
を更新する。あるいは、補助変数
2018136419.tif 000235
と、観測スペクトログラム
2018136419.tif 000236
と、
2018136419.tif 000237
とに基づいて、上記式（２３）を用いて、勾配法により、時間領域信号
2018136419.tif 000238
を更新する。
【０１７２】
また、パラメータ更新部５８は、補助変数
2018136419.tif 000239
と、観測スペクトログラム
2018136419.tif 000240
と、時間領域信号
2018136419.tif 000241
とに基づいて、上記式（４０）、式（４１）に従って、
2018136419.tif 000242
を更新する。あるいは、補助変数
2018136419.tif 000243
と、観測スペクトログラム
2018136419.tif 000244
と、時間領域信号
2018136419.tif 000245
とに基づいて、上記式（４７）、式（４８）に従って、
2018136419.tif 000246
を更新する。あるいは、補助変数
2018136419.tif 000247
と、観測スペクトログラム
2018136419.tif 000248
と、時間領域信号
2018136419.tif 000249
と、パラメータ記憶部４０に格納されているメル周波数ケプストラム係数のパラメータ
2018136419.tif 000250
とに基づいて、上記式（６０）、式（６１）に従って、
2018136419.tif 000251
を更新する。
【０１７３】
収束判定部６０は、予め定められた収束判定条件を満たすまで、補助変数更新部５４及びパラメータ更新部５８による各処理を繰り返させる。収束判定条件としては、例えば、予め定められた繰り返し回数に到達することである。
【０１７４】
収束判定条件を満たしたときに、最終的に得られた時間領域信号
2018136419.tif 000252
を、観測スペクトログラムに対して欠損成分を補完して構成された時間領域信号の推定結果として、出力部９０により出力する。
【０１７５】
＜信号推定装置の作用＞
【０１７６】
次に、本実施の形態に係る信号推定装置１００の作用について説明する。まず、学習サンプルの目的音響信号の時系列データが、信号推定装置１００に入力されると、信号推定装置１００において、目的音響信号のメル周波数ケプストラム係数（Mel-Frequency Cepstrum Coefficients: MFCC）が抽出され、抽出された目的音響信号のメル周波数ケプストラム係数に基づいて、メル周波数ケプストラム係数のパラメータ
2018136419.tif 000253
が学習され、パラメータ記憶部４０に格納される。
【０１７７】
そして、目的音響信号について、一部欠損した観測スペクトログラムが、信号推定装置１００に入力されると、信号推定装置１００において、図５に示す信号推定処理ルーチンが実行される。
【０１７８】
まず、ステップＳ１００において、基底スペクトル、及び各時刻におけるアクティベーションパラメータの各々の初期値を設定する。例えば、上記式（２）における第三項から、
2018136419.tif 000254
を推定し、これを固定したのち、上記式（２）における第一項及び第二項から時間領域信号
2018136419.tif 000255
を推定して初期値とする。
【０１７９】
そして、ステップＳ１０２において、時間領域信号
2018136419.tif 000256
に基づいて、上記式（２９）又は上記式（３０）、（３１）に従って、補助変数
2018136419.tif 000257
又は
2018136419.tif 000258
を更新する。あるいは、上記式（３２）〜（３４）に従って、補助変数
2018136419.tif 000259
を更新する。
【０１８０】
また、
2018136419.tif 000260
に基づいて、上記式（６８）又は上記式（６９）に従って、補助変数
2018136419.tif 000261
又は
2018136419.tif 000262
を更新する。あるいは、上記式（７０）又は上記式（７１）に従って、補助変数
2018136419.tif 000263
又は
2018136419.tif 000264
を更新する。あるいは、
2018136419.tif 000265
と、パラメータ記憶部４０に格納されているメル周波数ケプストラム係数のパラメータ
2018136419.tif 000266
と、
2018136419.tif 000267
の積のメル周波数ケプストラム係数
2018136419.tif 000268
とに基づいて、上記式（７２）〜式（７８）に従って、補助変数
2018136419.tif 000269
を更新する。
【０１８１】
そして、ステップＳ１０４において、補助変数
2018136419.tif 000270
と、観測スペクトログラム
2018136419.tif 000271
と、
2018136419.tif 000272
とに基づいて、上記式（２０）を用いて、勾配法により、時間領域信号
2018136419.tif 000273
を更新する。あるいは、補助変数
2018136419.tif 000274
と、観測スペクトログラム
2018136419.tif 000275
と、
2018136419.tif 000276
とに基づいて、上記式（２８）に従って、時間領域信号
2018136419.tif 000277
を更新する。あるいは、補助変数
2018136419.tif 000278
と、観測スペクトログラム
2018136419.tif 000279
と、
2018136419.tif 000280
とに基づいて、上記式（２３）を用いて、勾配法により、時間領域信号
2018136419.tif 000281
を更新する。
【０１８２】
そして、ステップＳ１０６において、補助変数
2018136419.tif 000282
と、観測スペクトログラム
2018136419.tif 000283
と、時間領域信号
2018136419.tif 000284
とに基づいて、上記式（４０）、式（４１）に従って、
2018136419.tif 000285
を更新する。あるいは、補助変数
2018136419.tif 000286
と、観測スペクトログラム
2018136419.tif 000287
と、時間領域信号
2018136419.tif 000288
とに基づいて、上記式（４７）、式（４８）に従って、
2018136419.tif 000289
を更新する。あるいは、補助変数
2018136419.tif 000290
と、観測スペクトログラム
2018136419.tif 000291
と、時間領域信号
2018136419.tif 000292
と、パラメータ記憶部４０に格納されているメル周波数ケプストラム係数のパラメータ
2018136419.tif 000293
とに基づいて、上記式（６０）、式（６１）に従って、
2018136419.tif 000294
を更新する。
【０１８３】
ステップＳ１０８において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１０２へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１１０へ進む。
【０１８４】
ステップＳ１１０では、上記ステップＳ１０４で最終的に得られた時間領域信号
2018136419.tif 000295
を、出力部９０により出力し、信号推定処理ルーチンを終了する。
【０１８５】
＜実験＞
ATR音声データベース503文の音声データと雑音データ（babble noise）を用いて上述した実施の形態の手法による雑音抑圧効果を検証する評価実験を行った。比較対象は従来のNMF法に加えて位相復元アルゴリズムを併用したものとし、信号対雑音比（SNR）およびケプストラム歪みを評価したテストデータは全て音響信号はサンプリング周波数16kHzのモノラル信号で、フレーム長32 ms、シフト長16msで短時間Fourier変換を行い、これに対して理想バイナリマスクによりマスキングを行うことで欠損成分の存在する観測スペクトログラムYk,mを算出した。学習においては男性1名の話者の計100文の音声を用いてMFCCのGMMパラメータの学習を行った。学習で得られたGMMのパラメータを固定し、正則化パラメータを各項が同等の大きさを持つよう逐次更新し、パラメータ
2018136419.tif 000296
の推定を行った。提案法アルゴリズムの初期段階では式(2)における第三項からパラメータ
2018136419.tif 000297
を推定し、これを固定したのち、第一項及び第二項からパラメータ
2018136419.tif 000298
を推定して初期値とした。
【０１８６】
以上の条件下で未処理（理想バイナリマスク）、正則化の有無を含めた提案法及び従来法によって得られたSNR及びケプストラム歪みを図６（Ａ）、（Ｂ）に示す。横軸は理想バイナリマスク作成の際の音声データと雑音データのSNRを示す。いずれの評価尺度においても提案法の方が高い性能が得られることが確認できた。
【０１８７】
以上説明したように、本発明の実施の形態に係る信号推定装置によれば、欠損していない部分における、観測スペクトログラムと、時間領域信号から求められるスペクトログラムとの誤差と、時間領域信号の基底スペクトル、及びアクティベーションパラメータの積と、時間領域信号から求められるスペクトログラムの絶対値との誤差、欠損していない部分における、基底スペクトル、及びアクティベーションパラメータの積と、観測スペクトログラムの絶対値との誤差、及び基底スペクトル及び各時刻におけるアクティベーションパラメータの積のケプストラム特徴量の事前確率とを用いて表される規準を最適化するように、時間領域信号、基底スペクトル、及び各時刻におけるアクティベーションパラメータを推定することにより、一部欠損した観測スペクトログラムから時間領域信号を精度良く推定することができる。
【０１８８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１８９】
例えば、上述の信号推定装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１９０】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１９１】
１０   入力部
２０   演算部
３０   学習部
３２   特徴量抽出部
３６   パラメータ学習部
４０   パラメータ記憶部
５０   信号推定部
５２   初期値設定部
５４   補助変数更新部
５８   パラメータ更新部
６０   収束判定部
９０   出力部
１００ 信号推定装置
【図１】
2018136419.tif 000299
【図２】
2018136419.tif 000300
【図３】
2018136419.tif 000301
【図４】
2018136419.tif 000302
【図５】
2018136419.tif 000303
【図６】
2018136419.tif 000304