(57)【要約】【課題】音声信号から声道スペクトルを精度よく推定することができる。【解決手段】時間周波数展開部２４が、音声信号の時系列データを入力として、観測スペクトログラムを出力する。パラメータ推定部３６が、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、Ｋ個のフォルマント周波数のセットに対するフォルマント強度と結合係数で表し、分散を、Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の混合Ｇａｕｓｓ関数モデルと、観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、フォルマント強度と、フォルマント分散及び各時刻の結合係数とを推定する。【選択図】図５
2018141924.tif 000002
【特許請求の範囲】
【請求項１】
音声信号の時系列データを入力として、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを出力する時間周波数展開部と、
前記時間周波数展開部により出力された前記観測スペクトログラムに基づいて、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数で近似した混合Ｇａｕｓｓ関数モデルであって、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と前記結合係数で表し、分散を、前記Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の前記混合Ｇａｕｓｓ関数モデルと、前記観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを推定するパラメータ推定部と、
を含む声道スペクトル推定装置。
【請求項２】
前記規準は、前記Ｇａｕｓｓ分布関数の各々についての、Ｋ個のフォルマント周波数のセットに対するフォルマント強度の大きさに関するペナルティ項を更に含む請求項１記載の声道スペクトル推定装置。
【請求項３】
前記パラメータ推定部は、
前記規準の上界関数である補助関数を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを更新するパラメータ更新部と、
予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させる収束判定部と、
を含む請求項１又は２記載の声道スペクトル推定装置。
【請求項４】
時間周波数展開部が、音声信号の時系列データを入力として、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを出力し、
パラメータ推定部が、前記時間周波数展開部により出力された前記観測スペクトログラムに基づいて、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数で近似した混合Ｇａｕｓｓ関数モデルであって、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と前記結合係数で表し、分散を、前記Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の前記混合Ｇａｕｓｓ関数モデルと、前記観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを推定する
声道スペクトル推定方法。
【請求項５】
前記規準は、前記Ｇａｕｓｓ分布関数の各々についての、Ｋ個のフォルマント周波数のセットに対するフォルマント強度の大きさに関するペナルティ項を更に含む請求項４記載の声道スペクトル推定方法。
【請求項６】
前記パラメータ推定部が推定することでは、
パラメータ更新部が、前記規準の上界関数である補助関数を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを更新し、
収束判定部が、予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させる
ことを含む請求項４又は５記載の声道スペクトル推定方法。
【請求項７】
請求項１〜請求項３の何れか１項に記載の声道スペクトル推定装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、声道スペクトル推定装置、方法、及びプログラムに係り、特に、音声信号から声道スペクトルを推定する声道スペクトル推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音声のソースフィルタ理論（図１）によればスペクトル包絡は声道の共振特性（声道スペクトル）に対応しており、声質や音韻などに関する情報を豊富に含んでいる。通常発話では母音や音素の種類は限られていることから、各時刻の声道スペクトルを有限個のテンプレートの非負結合で近似できると仮定すると、声道スペクトログラムを二つの非負値行列の積で表現することができる。非負値行列を二つの非負値行列の積で近似することを非負値行列因子分解(Non-negative Matrix Factorization; NMF) といい、声道スペクトログラムのNMF による行列積表現はいくつかの応用場面において有用である。以下に二つの例を挙げる。
【０００３】
まず、この表現により声道スペクトログラムを時間に依存しない話者の声質に関する因子と時間依存する発話内容や発話スタイルに関する因子とに分解することができるため、後者の因子を固定したまま前者の因子のみを変換してスペクトログラムを再構成することで、音声の発話内容を変えずに声質のみを変えることができる（非特許文献１、２）。このような技術を声質変換という。
【０００４】
もう一つは、音声分析合成のための声道スペクトル推定への応用である。音声合成や音声変換をはじめ音声処理全般において、音声信号から声道スペクトルを推定する技術は多くの場面で用いられている。短区間ごとの音声信号を周期デルタ関数（パルス列）を入力とした線形時不変系の出力としてモデル化できるとすれば、この線形系の入力とインパルス応答がそれぞれ声帯音源信号と声道特性に対応する。この仮定は、周波数領域では周期デルタ関数で表される声帯音源スペクトルと声道スペクトルの積で音声スペクトルが表されることに相当する。従って、音声スペクトルは声道スペクトルを周期的に（基本周波数間隔で）サンプリングしたもの、と見なすことができる。代表的な声道スペクトル推定方法の一つとして広く知られるSTRAIGHT は、音声信号を基本周期の幅で切り出し、その切り出し信号のスペクトルを声道スペクトルの推定値とする方法である。これは周波数領域では、各調波成分のピークを滑らかに補間したものを声道スペクトルと見なしていることに相当する。しかしこの方法によって得られる声道スペクトル推定値は、定常な音声が対象であっても切り出しフレームのオフセットに依存して周期的に時間変化することが知られる。これは各調波成分が互いに干渉し合うからであり、調波成分間の干渉に起因するこの変動成分を除くよう改良された手法も提案されている。しかし、前述のとおり音声スペクトルは声道スペクトルを基本周波数(F0) 間隔でサンプリングしたものと見なせるため、音声のF0 が高いときほど声道スペクトル推定の手がかりは少なくなる。このことは、フレームごとに独立な処理に本質的な限界があることを示唆している。一方で、通常発話では母音や音素の種類は限られていることから、類似した声道スペクトルが複数の異なる時刻で現れることも手がかりとなる。複数のフレームが共通の声道スペクトルを持つと仮定できそれらのフレームでF0 が異なれば、実際に観測可能な声道スペクトルのサンプル点が単一のフレームの場合よりも増えるため、原理的には声道スペクトルの「超解像」推定が可能になる。この考え方に基づき、音声スペクトログラムの中の基本周波数成分や高調波成分にのみ適合するようにNMF を行うことで複数フレームの音声スペクトルを手がかりに高精度に声道スペクトルを推定する手法が提案されている（非特許文献３）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】R. Takashima、 R. Aihara、 T. Takiguchi、 and Y. Ariki、 “Exemplar-based voice conversion usin sparse representation in noisy environments、" IEICE Transactions on Information and Systems、 vol. E96-A、 no. 10、pp. 1946-1953、 2013.
【非特許文献２】Z. Wu、 T. Virtanen、 T. Kinnunen、 E.S. Chng、 and H. Li、“Exemplar-based voice conversion using non-negative spectrogram deconvolution、" Proc. 8th ISCA Speech Synthesis Workshop、 pp. 201-206、 2013.
【非特許文献３】中村、 亀岡、 "非負値行列因子分解に基づく欠損データ補間による声道スペクトル推定法の検討、" 音講論(春)、3-P-33， 393-396， 2016.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
以上のように、NMF による声道スペクトログラムの行列積表現は音声信号処理の各種重要課題に対し有用な解決を与えるポテンシャルを有している。しかしながら、この表現では声道スペクトルをテンプレートの線形結合で表現することしかできないため、図２ に示すような声道スペクトログラムに現れる連続的なピーク軌跡を表現することが難しい（図３）。この声道スペクトルに現れるいくつかのピークをフォルマントといい、それぞれの周波数は声道の共振周波数に対応しており、母音の種類や質を特徴付ける量と考えられている。
【０００７】
本発明では、上記事情を鑑みて成されたものであり、音声信号から声道スペクトルを精度よく推定することができる声道スペクトル推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、本発明に係る声道スペクトル推定装置は、音声信号の時系列データを入力として、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを出力する時間周波数展開部と、前記時間周波数展開部により出力された前記観測スペクトログラムに基づいて、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数で近似した混合Ｇａｕｓｓ関数モデルであって、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と前記結合係数で表し、分散を、前記Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の前記混合Ｇａｕｓｓ関数モデルと、前記観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを推定するパラメータ推定部と、を含んで構成されている。
【０００９】
本発明に係る声道スペクトル推定方法は、時間周波数展開部が、音声信号の時系列データを入力として、各時刻及び各周波数の信号の成分を表す観測スペクトログラムを出力し、パラメータ推定部が、前記時間周波数展開部により出力された前記観測スペクトログラムに基づいて、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数で近似した混合Ｇａｕｓｓ関数モデルであって、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と前記結合係数で表し、分散を、前記Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の前記混合Ｇａｕｓｓ関数モデルと、前記観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを推定する。
【００１０】
また、本発明のプログラムは、コンピュータを、上記の声道スペクトル推定装置を構成する各部として機能させるためのプログラムである。
【発明の効果】
【００１１】
以上説明したように、本発明の声道スペクトル推定装置、方法、及びプログラムによれば、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数で近似した混合Ｇａｕｓｓ関数モデルであって、各ピークに対応するＧａｕｓｓ分布関数の各々に対し、平均を、Ｋ個のフォルマント周波数のセットと結合係数で表し、重みを、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と前記結合係数で表し、分散を、前記Ｋ個のフォルマント周波数のセットに対するフォルマント分散と結合係数で表したときの、各時刻の前記混合Ｇａｕｓｓ関数モデルと、前記観測スペクトログラムとの誤差の大きさを表す規準を小さくするように、前記Ｇａｕｓｓ分布関数の各々についての、前記Ｋ個のフォルマント周波数のセット及び各時刻の結合係数と、前記Ｋ個のフォルマント周波数のセットに対するフォルマント強度と、Ｋ個のフォルマント周波数のセットに対するフォルマント分散及び各時刻の結合係数とを推定することにより、音声信号から声道スペクトルを精度よく推定することができる。
【図面の簡単な説明】
【００１２】
【図１】ソースフィルタ理論を説明するための図である。
【図２】STRAIGHT 分析によって推定された声道スペクトログラムを示す図である。
【図３】推定された声道スペクトログラムを近似した非負値行列積を示す図である。
【図４】GMM によるスペクトル包絡の表現を示す図である。
【図５】本発明の実施の形態に係る声道スペクトル推定装置の機能的構成を示すブロック図である。
【図６】本発明の実施の形態に係る声道スペクトル推定装置におけるパラメータ推定処理ルーチンを示すフローチャート図である。
【図７】STRAIGHT 分析によって推定された声道スペクトログラムを示す図である。
【図８】図７の声道スペクトログラムに対し、本実施の形態に係る手法により推定した推定結果を示す図である。
【図９】STRAIGHT 分析によって推定された声道スペクトログラムを示す図である。
【図１０】図９の声道スペクトログラムに対し、本実施の形態に係る手法により推定した推定結果を示す図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１４】
＜本実施の形態の概要＞
本発明の実施の形態では、スペクトルテンプレートの線形結合による従来の声道スペクトログラム表現の代わりに、フォルマント周波数セットとフォルマント周波数セットのテンプレートを考え、それらの線形結合による各時刻のフォルマント周波数セットおよびフォルマント周波数セットの表現を内部にもつ新しい声道スペクトログラム表現を提案し、観測音声スペクトログラムから各テンプレートと各時刻におけるそれらの結合係数を推定する分析アルゴリズムを提案する。具体的には、スペクトル包絡の各ピークをGauss 分布関数で近似した混合Gauss 関数モデル（Gaussian Mixture Model; GMM）を考え、各Gauss 関数のピーク周波数と重みをテンプレートの線形結合で表すモデルを構築し、このスペクトル包絡モデルを観測音声スペクトルにできるだけフィットするようにフォルマント周波数セットおよびフォルマント周波数セットのテンプレートとそれらの結合係数を推定する最適化アルゴリズムを提案する。
【００１５】
＜声道スペクトログラムモデル＞
母音の種類や質はフォルマント周波数セットによって特徴付けられる。図２からも分かるように、発話中フォルマント周波数は連続的に時間変化する傾向にある。これは発声器官の物理的制約による。また、通常の発話で使用する母音の種類は限られているため、各時刻のフォルマント周波数セットは各母音に対応した有限個のテンプレートの凸結合で表せると仮定することができる。これは、各時刻のフォルマント周波数セットを各列ベクトルとした行列を二つの行列積で表現することに相当する。以後、これをフォルマント周波数行列積表現と呼ぶ。以下では、以上のフォルマント周波数行列積表現を組み込んだ声道スペクトログラムモデルを構築する。
【００１６】
スペクトル包絡の各ピークをGauss 分布関数で近似した混合Gauss 関数モデル（Gaussian Mixture Model; GMM）F(ω，t)
【００１７】
【数１】
2018141924.tif 000003
2018141924.tif 000004
【００１８】
を考える（図４）。ただし、ωは周波数、t は時刻、I +1 はGauss 関数
2018141924.tif 000005
の個数、
2018141924.tif 000006
はそれぞれi 番目のGauss 関数の平均、分散、重みである。また、
2018141924.tif 000007
は
【００１９】
【数２】
2018141924.tif 000008
【００２０】
を満たすものとする。
2018141924.tif 000009
に2 を乗じているのは、
2018141924.tif 000010
とするためである。なお、
2018141924.tif 000011
のときは
2018141924.tif 000012
である。よって、式(3) の条件では
2018141924.tif 000013
である。従って、β(t) は時刻t におけるスペクトルのスケールを表す。
【００２１】
k 番目の母音テンプレートのフォルマント周波数を
2018141924.tif 000014
フォルマント強度を
2018141924.tif 000015
とする。前述の議論に従い、
2018141924.tif 000016
を
【００２２】
【数３】
2018141924.tif 000017
2018141924.tif 000018
【００２３】
のようにK 個のフォルマント周波数セットおよびフォルマント強度セットのテンプレートの凸結合で表現する。ただし、
【００２４】
【数４】
2018141924.tif 000019
2018141924.tif 000020
【００２５】
である。
2018141924.tif 000021
についても同様に、フォルマント分散
2018141924.tif 000022
を用いて
【００２６】
【数５】
2018141924.tif 000023
【００２７】
と表現する。以上の
2018141924.tif 000024
が推定すべき未知パラメータである。ここでいずれの変数も非負でなければならない点に注意が必要である。
【００２８】
フォルマント周波数やフォルマント強度は通常連続的に時間変化する。そこで、
2018141924.tif 000025
が滑らかな関数となることを保証するため、時刻tj を中心に局在する滑らかな非負関数（Gauss 関数やHanning関数など）
2018141924.tif 000026
を用い
2018141924.tif 000027
を
【００２９】
【数６】
2018141924.tif 000028
【００３０】
と表現し、
2018141924.tif 000029
の代わりに
2018141924.tif 000030
を推定すべきパラメータとすることもできる。ただし、
【００３１】
【数７】
2018141924.tif 000031
2018141924.tif 000032
とする。
【００３２】
＜補助関数法によるパラメータ推定アルゴリズム＞
＜最適化問題の定式化＞
観測スペクトログラムを
2018141924.tif 000033
とし、
2018141924.tif 000034
と
2018141924.tif 000035
ができるだけ近くなるように未知パラメータ
2018141924.tif 000036
2018141924.tif 000037
を求めるための最適化アルゴリズムについて述べる。ここでは
2018141924.tif 000038
と
2018141924.tif 000039
の近さは例えばI ダイバージェンス
【００３３】
【数８】
2018141924.tif 000040
2018141924.tif 000041
【００３４】
で測ることができる。ただし、
2018141924.tif 000042
は上記未知パラメータの集合とする。また、
2018141924.tif 000043
が小さくなりすぎなくする目的で
【００３５】
【数９】
2018141924.tif 000044
【００３６】
のようなペナルティ関数（対称Dirichlet 分布の対数の負）を考える。対称Dirichlet 分布は変数がすべて等しいときに最大となるため、式(14) は
2018141924.tif 000045
のときに最小となる。以下では、
2018141924.tif 000046
を式(9) で表現する場合を考える（
2018141924.tif 000047
をデルタ関数
2018141924.tif 000048
に置き換えれば
2018141924.tif 000049
を変数とする場合と等価である）。以上より、最適化問題は以下となる。
【００３７】
【数１０】
2018141924.tif 000050
【００３８】
以上の最適化問題の停留点は射影勾配法などを使って探索することもできるが、ここでは、補助関数法の原理に基づき停留点への収束性が保証された最適化アルゴリズムを導出する。
【００３９】
＜補助関数法＞
2018141924.tif 000051
を
2018141924.tif 000052
に関して最小化したい目的関数とすると、
2018141924.tif 000053
を満たす関数
2018141924.tif 000054
を補助関数、αを補助変数と呼ぶ。このような補助関数を設計できれば、 
2018141924.tif 000055
と
2018141924.tif 000056
を交互に繰り返すことで、目的関数
2018141924.tif 000057
の停留点を得ることができる。この最適化手法を補助関数法と呼ぶ。
【００４０】
＜補助関数の設計＞
まず、負の対数関数は凸関数よりJensen の不等式より、
【００４１】
【数１１】
2018141924.tif 000058
【００４２】
が言える。ただし、
2018141924.tif 000059
は
2018141924.tif 000060
を満たす非負の変数で、等号は
【００４３】
【数１２】
2018141924.tif 000061
【００４４】
のとき成立する。同様に、
【００４５】
【数１３】
2018141924.tif 000062
【００４６】
が言える。ただし、
2018141924.tif 000063
は
2018141924.tif 000064
を満たす非負の変数で、等号は
【００４７】
【数１４】
2018141924.tif 000065
【００４８】
のとき成立する。次に、正の対数関数は凹関数より、
【００４９】
【数１５】
2018141924.tif 000066
【００５０】
が言える。ただし、
2018141924.tif 000067
は非負の変数で、等号は
【００５１】
【数１６】
2018141924.tif 000068
【００５２】
のとき成立する。二次関数は凸関数のため再びJensen の不等式より
【００５３】
【数１７】
2018141924.tif 000069
【００５４】
が言える。ただし、
2018141924.tif 000070
は
2018141924.tif 000071
を満たす非負の変数で、等号は
【００５５】
【数１８】
2018141924.tif 000072
【００５６】
のとき成立する。最後に、逆数関数は凸関数のためJensen の不等式より
【００５７】
【数１９】
2018141924.tif 000073
【００５８】
が言える。ただし、
2018141924.tif 000074
は
2018141924.tif 000075
を満たす非負の変数で、等号は
【００５９】
【数２０】
2018141924.tif 000076
【００６０】
のとき成立する。
【００６１】
以上をまとめると、
【００６２】
【数２１】
2018141924.tif 000077
【００６３】
が言え、右辺は
2018141924.tif 000078
の補助関数としての要件を満たす。ただし、
2018141924.tif 000079
は補助変数
2018141924.tif 000080
の集合である。また、
2018141924.tif 000081
はパラメータ
2018141924.tif 000082
にも補助変数
2018141924.tif 000083
にも依らない定数項をまとめたものである。
【００６４】
＜パラメータ更新式＞
以上の補助関数
2018141924.tif 000084
を用いて各パラメータの更新式を導くことができる。
【００６５】
＜ｗi,kの更新式＞
【００６６】
【数２２】
2018141924.tif 000085
【００６７】
より、
2018141924.tif 000086
の更新式
【００６８】
【数２３】
2018141924.tif 000087
【００６９】
を得る。
【００７０】
＜ａi,kの更新式＞
2018141924.tif 000088
は式(7) の条件を満たす必要があるので、ラグランジアン
2018141924.tif 000089
の偏微分を
【００７１】
【数２４】
2018141924.tif 000090
【００７２】
と置くことにより、
【００７３】
【数２５】
2018141924.tif 000091
2018141924.tif 000092
2018141924.tif 000093
【００７４】
より、
2018141924.tif 000094
の更新式
【００７５】
【数２６】
2018141924.tif 000095
【００７６】
を得る。
【００７７】
＜ｌj,kの更新式＞
2018141924.tif 000096
は式(１１)の条件を満たす必要があるが、ここでは無制約下で
2018141924.tif 000097
を最小化する
2018141924.tif 000098
を求めたのちに正規化する方法をとることとする。
2018141924.tif 000099
の偏微分を
【００７８】
【数２７】
2018141924.tif 000100
【００７９】
と置くことにより、
【００８０】
【数２８】
2018141924.tif 000101
【００８１】
を得る。よって、
2018141924.tif 000102
を最小化する
2018141924.tif 000103
は
【００８２】
【数２９】
2018141924.tif 000104
【００８３】
となる。式(38) で
2018141924.tif 000105
を更新した後、
【００８４】
【数３０】
2018141924.tif 000106
【００８５】
により正規化する。
【００８６】
＜ｃi,kの更新式＞
【００８７】
【数３１】
2018141924.tif 000107
【００８８】
より、
2018141924.tif 000108
の更新式
【００８９】
【数３２】
2018141924.tif 000109
を得る。
【００９０】
＜ｄk（ｔ）の更新式＞
【００９１】
【数３３】
2018141924.tif 000110
【００９２】
より、
2018141924.tif 000111
の更新式
【００９３】
【数３４】
2018141924.tif 000112
を得る。
【００９４】
＜β（ｔ）の更新式＞
【００９５】
【数３５】
2018141924.tif 000113
【００９６】
より、
2018141924.tif 000114
の更新式
【００９７】
【数３６】
2018141924.tif 000115
【００９８】
を得る。
【００９９】
＜最適化アルゴリズム＞
補助関数法に基づく以上の最適化アルゴリズムは以下の１．〜３．のようにまとめられる。
【０１００】
1. 
2018141924.tif 000116
を初期設定する。
【０１０１】
2. 補助変数を更新する。
【数３７】
2018141924.tif 000117
2018141924.tif 000118
2018141924.tif 000119
2018141924.tif 000120
2018141924.tif 000121
【０１０２】
3. パラメータを更新する。
【数３８】
2018141924.tif 000122
2018141924.tif 000123
2018141924.tif 000124
2018141924.tif 000125
2018141924.tif 000126
2018141924.tif 000127
2018141924.tif 000128
【０１０３】
＜本発明の実施の形態に係る声道スペクトル推定装置の構成＞
次に、本発明の実施の形態に係る声道スペクトル推定装置の構成について説明する。図５に示すように、本発明の実施の形態に係る声道スペクトル推定装置１００は、ＣＰＵと、ＲＡＭと、後述するパラメータ推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することができる。この声道スペクトル推定装置１００は、機能的には図５に示すように入力部１０と、演算部２０と、出力部９０と、を含んで構成されている。
【０１０４】
入力部１０は、音声信号の時系列データを受け付ける。
【０１０５】
演算部２０は、時間周波数展開部２４と、パラメータ推定部３６とを含んで構成されている。
【０１０６】
時間周波数展開部２４は、音声信号の時系列データに基づいて、各時刻における各周波数の信号の成分を表す振幅スペクトログラム又はパワースペクトログラムである観測スペクトログラムＹを計算する。なお、本実施の形態においては、短時間フーリエ変換やウェーブレット変換などの時間周波数展開を行う。
【０１０７】
パラメータ推定部３６は、時間周波数展開部２４により出力された観測スペクトログラムＹに基づいて、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数Ｇi（ω、ｔ）で近似した混合Ｇａｕｓｓ関数モデルＦ（ω、ｔ）であって、各ピークに対応するＧａｕｓｓ分布関数Ｇi（ω、ｔ）の各々に対し、平均μi（ｔ）を、Ｋ個のフォルマント周波数のセットｗi,1,・・・, ｗi,kと結合係数ｈk（ｔ）で表し、重みαi（ｔ）を、Ｋ個のフォルマント周波数のセットに対するフォルマント強度ａi,1,・・・, ａi,kと結合係数ｈk（ｔ）で表し、分散σi2（ｔ）を、Ｋ個のフォルマント周波数のセットに対するフォルマント分散ｃi,1,・・・, ｃi,kと結合係数ｄk（ｔ）で表したときの、各時刻の混合Ｇａｕｓｓ関数モデルＦ（ω、ｔ）と、観測スペクトログラムＹ（ω、ｔ）との誤差の大きさを表す、上記（１５）式の規準を小さくするように、Ｇａｕｓｓ分布関数Ｇi（ω、ｔ）の各々についての、Ｋ個のフォルマント周波数のセットｗi,1,・・・, ｗi,k及び各時刻の結合係数ｈk（ｔ）を表すためのパラメータｌj,kと、Ｋ個のフォルマント周波数のセットに対するフォルマント強度ａi,1,・・・, ａi,kと、Ｋ個のフォルマント周波数のセットに対するフォルマント分散ｃi,1,・・・, ｃi,k及び各時刻の結合係数ｈk（ｔ）と、各時刻の重みβ（ｔ）とを推定する。
【０１０８】
具体的には、パラメータ推定部３６は、初期値設定部４０、補助変数更新部４２、パラメータ更新部４４、及び収束判定部４６を備えている。
【０１０９】
初期値設定部４０は、未知パラメータ
2018141924.tif 000129
の各々に初期値を設定する。
【０１１０】
補助変数更新部４２は、初期値である、又は前回更新したパラメータ
2018141924.tif 000130
に基づいて、上記式（４６）〜式（５０）に従って、補助変数
2018141924.tif 000131
を更新する。
【０１１１】
パラメータ更新部４４は、時間周波数展開部２４により出力された観測スペクトログラムＹと、補助変数更新部４２によって更新された補助変数
2018141924.tif 000132
と、初期値である、又は前回更新したパラメータ
2018141924.tif 000133
とに基づいて、上記式（５１）〜式（５７）に従って、パラメータ
2018141924.tif 000134
を更新する。
【０１１２】
収束判定部４６は、収束条件を満たすか否かを判定し、収束条件を満たすまで、補助変数更新部４２における更新処理と、パラメータ更新部４４における更新処理とを繰り返させ、最終的に推定されたパラメータを、出力部９０から出力する。
【０１１３】
収束条件としては、例えば、繰り返し回数が、上限回数に到達したことを用いることができる。あるいは、収束条件として、上記式（１５）の規準の値と前回の規準の値との差分が、予め定められた閾値以下であることを用いることができる。
【０１１４】
＜本発明の実施の形態に係る声道スペクトル推定装置の作用＞
次に、本発明の実施の形態に係る声道スペクトル推定装置１００の作用について説明する。まず、入力部１０において音声信号の時系列データを受け付けると、声道スペクトル推定装置１００は、図６に示すパラメータ推定処理ルーチンを実行する。
【０１１５】
まず、ステップＳ１００では、入力部１０において受け付けた音声信号の時系列データに基づいて、観測スペクトログラムＹを計算する。
【０１１６】
ステップＳ１０２では、未知パラメータ
2018141924.tif 000135
の各々に初期値を設定する。
【０１１７】
ステップＳ１０４では、初期値である、又は前回更新したパラメータ
2018141924.tif 000136
に基づいて、上記式（４６）〜式（５０）に従って、補助変数
2018141924.tif 000137
を更新する。
【０１１８】
次に、ステップＳ１０６では、上記ステップＳ１００で得られた観測スペクトログラムＹと、補助変数更新部４２によって更新された補助変数
2018141924.tif 000138
と、初期値である、又は前回更新したパラメータ
2018141924.tif 000139
とに基づいて、上記式（５１）〜式（５７）に従って、パラメータ
2018141924.tif 000140
を更新する。
【０１１９】
次に、ステップＳ１０８では、収束条件を満たすか否かを判定する。収束条件を満たした場合には、ステップＳ１１０へ移行し、収束条件を満たしていない場合には、ステップＳ１０４へ移行し、ステップＳ１０４〜ステップＳ１０６の処理を繰り返す。
【０１２０】
ステップＳ１１０では、上記ステップＳ１０６で最終的に更新されたパラメータ
2018141924.tif 000141
を出力部９０から出力して、パラメータ推定処理ルーチンを終了する。
【０１２１】
＜実験例＞
図７、９ に示す
2018141924.tif 000142
に対し、本発明の実施の形態に係る手法により推定した
2018141924.tif 000143
を図８、１０に示す。テンプレート数はK = 10 とした。
2018141924.tif 000144
のフォルマント周波数軌跡を適切に推定できていることが確認できた。
【０１２２】
以上説明したように、本発明の実施の形態に係る声道スペクトル推定装置によれば、各時刻のスペクトル包絡の各ピークをＧａｕｓｓ分布関数Ｇi（ω、ｔ）で近似した混合Ｇａｕｓｓ関数モデルＦ（ω、ｔ）であって、各ピークに対応するＧａｕｓｓ分布関数Ｇi（ω、ｔ）の各々に対し、平均μi（ｔ）を、Ｋ個のフォルマント周波数のセットｗi,1,・・・, ｗi,kと結合係数ｈk（ｔ）で表し、重みαi（ｔ）を、Ｋ個のフォルマント周波数のセットに対するフォルマント強度ａi,1,・・・, ａi,kと結合係数ｈk（ｔ）で表し、分散σi2（ｔ）を、Ｋ個のフォルマント周波数のセットに対するフォルマント分散ｃi,1,・・・, ｃi,kと結合係数ｄk（ｔ）で表したときの、各時刻の混合Ｇａｕｓｓ関数モデルＦ（ω、ｔ）と、観測スペクトログラムＹ（ω、ｔ）との誤差の大きさを表す規準を小さくするように、Ｇａｕｓｓ分布関数Ｇi（ω、ｔ）の各々についての、Ｋ個のフォルマント周波数のセットｗi,1,・・・, ｗi,k及び各時刻の結合係数ｈk（ｔ）を表すためのパラメータｌj,kと、Ｋ個のフォルマント周波数のセットに対するフォルマント強度ａi,1,・・・, ａi,kと、Ｋ個のフォルマント周波数のセットに対するフォルマント分散ｃi,1,・・・, ｃi,k及び各時刻の結合係数ｈk（ｔ）と、各時刻の重みβ（ｔ）とを推定することにより、音声信号から声道スペクトルを精度よく推定することができる。
【０１２３】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１２４】
また、更新するパラメータの順番には任意性があるため、上記の実施の形態の順番に限定されない。
【０１２５】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。
【符号の説明】
【０１２６】
１０   入力部
２０   演算部
２４   時間周波数展開部
３６   パラメータ推定部
４０   初期値設定部
４２   補助変数更新部
４４   パラメータ更新部
４６   収束判定部
９０   出力部
１００ 声道スペクトル推定装置
【図１】
2018141924.tif 000145
【図２】
2018141924.tif 000146
【図３】
2018141924.tif 000147
【図４】
2018141924.tif 000148
【図５】
2018141924.tif 000149
【図６】
2018141924.tif 000150
【図７】
2018141924.tif 000151
【図８】
2018141924.tif 000152
【図９】
2018141924.tif 000153
【図１０】
2018141924.tif 000154