(57)【特許請求の範囲】
【請求項１】
学習サンプルのソース音声の時系列データとターゲット音声の時系列データとからなるパラレルデータを入力として、前記ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の時系列データから抽出される、各時刻の基本周波数とに基づいて、前記ソース音声の各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータを、前記第１確率分布と前記第２確率分布とを用いて表される規準を大きくするように学習する学習部と、
予測対象のソース音声の時系列データを入力として、前記予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記学習部によって学習された前記第１確率分布のパラメータと前記第２確率分布のパラメータとに基づいて、前記規準を大きくするように、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数を予測する変換処理部と、
を含む基本周波数パターン予測装置。
【請求項２】
前記第１確率分布を、前記ソース音声の各時刻のスペクトル特徴量ベクトルと前記ターゲット音声の各時刻の基本周波数と前記基本周波数の動的成分との同時確率分布を表す混合正規分布とし、
前記第２確率分布を、
各時刻の基本周波数と、
隠れマルコフモデルの各時刻の状態からなる状態系列、又は各時刻における甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令のペアからなる指令関数との組み合わせの確率分布とした請求項１記載の基本周波数パターン予測装置。
【請求項３】
前記規準は、以下の式で表わされる請求項２記載の基本周波数パターン予測装置。
【数１】
0006472005.tif 000075
ただし、ｃは、前記ソース音声の各時刻のスペクトル特徴量ベクトルであり、ｙは、前記ターゲット音声の各時刻の基本周波数であり、ｍは、各時刻のスペクトル特徴量ベクトルを生成したらしい確率が最も高い正規分布のインデックスであり、ｓは、前記状態系列であり、γは、前記第１確率分布のパラメータであり、θは、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータであり、φは、状態遷移確率であり、ｑは、前記ターゲット音声の各時刻の基本周波数及び動的成分の結合ベクトルであり、ｕは、各時刻のフレーズ指令及びアクセント指令である。
【請求項４】
前記学習部は、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記第２確率分布から求められる、各時刻の基本周波数と、隠れマルコフモデルの各時刻の状態からなる状態系列との尤もらしさが大きくなるように、前記第２確率分布のパラメータとして、前記状態系列の状態遷移確率、及び各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータを学習するか、又は
各時刻の基本周波数が与えられたときの、各時刻のフレーズ指令及びアクセント指令のペアからなる指令関数及び前記振幅を表すパラメータの対数事後確率を目的関数として、前記目的関数を増加させるように、前記第２確率分布のパラメータとして、前記状態系列の状態遷移確率、及び前記振幅を表すパラメータを、前記第２確率分布のパラメータとして学習する請求項１〜請求項３の何れか１項記載の基本周波数パターン予測装置。
【請求項５】
前記学習部は、前記学習された前記第２確率分布のパラメータを固定して、ＥＭ（Expectation-Maximization）アルゴリズムにより、前記規準が大きくなるように、前記第１確率分布のパラメータを学習する請求項４記載の基本周波数パターン予測装置。
【請求項６】
前記変換処理部は、
前記予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記学習部によって学習された前記第１確率分布のパラメータと前記第２確率分布のパラメータとに基づいて、ＥＭアルゴリズムにより、前記規準を大きくするように、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数、各時刻のスペクトル特徴量ベクトルを生成したらしい確率が最も高い正規分布のインデックス、及び前記状態系列を推定することにより、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数を予測する請求項２記載の基本周波数パターン予測装置。
【請求項７】
学習部と、変換処理部とを含む基本周波数パターン予測装置における基本周波数パターン予測方法であって、
前記学習部が、学習サンプルのソース音声の時系列データとターゲット音声の時系列データとからなるパラレルデータを入力として、前記ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の時系列データから抽出される、各時刻の基本周波数とに基づいて、前記ソース音声の各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータを、前記第１確率分布と前記第２確率分布とを用いて表される規準を大きくするように学習し、
前記変換処理部が、予測対象のソース音声の時系列データを入力として、前記予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記学習部によって学習された前記第１確率分布のパラメータと前記第２確率分布のパラメータとに基づいて、前記規準を大きくするように、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数を予測する
基本周波数パターン予測方法。
【請求項８】
請求項１〜請求項６の何れか１項に記載の基本周波数パターン予測装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、基本周波数パターン予測装置、方法、及びプログラムに係り、特に、ソース音声から、ターゲット音声の基本周波数パターンを予測する基本周波数パターン予測装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
他者とのコミュニケーションにおいて音声は利便性に優れた手段ではあるが、時として物理的制約により様々な障壁が必然的にもたらされる。例えば、発声器官の内、わずか一か所でも正常に動作しなくなると、深刻な発声障害を患い、音声コミュニケーションに支障をきたす。また、音声生成という物理的行為は、秘匿性の高い意思伝達には不向きであるし、周囲の騒音に脆弱である。これらの障壁を無くすためには、身体的制約を超えて発声器官を動作させて音声を生成したり、適切な発音動作を指定して音声を生成したり、聴取困難なほど微かな音声発声時の発声器官動作から通常音声を生成するなど、物理的・身体的制約を超えた音声生成機能の拡張が必要である。
【０００３】
例えば、喉頭癌などで喉頭を失った喉頭摘出者に対して、残存器官を用いた代替発声法により生成される自然性に乏しい音声を、より自然な音声へと変換する発声補助技術が提案されている（非特許文献１〜非特許文献３を参照）。この他にも、非可聴つぶやき音声を自然な音声に変換する技術も提案されており、秘匿性に優れた通話技術としての応用が期待されている。上述の技術はいずれも音声のスペクトル特徴量系列から自然音声の基本周波数(F0) パターンを予測する問題を扱っている点で共通しており、学習処理と変換処理で構成される。学習処理では、対象音声（前者であれば電気音声、後者であれば非可聴つぶやき音声）と通常音声の同一発話データを用いる。まず各離散時刻（以後、フレーム）において、前後数フレームから得られる対象音声のスペクトル特徴量と、通常音声の対数F0とその動的成分（時間微分または時間差分）を抽出し、スペクトル距離尺度に基づく動的時間伸縮によりこれらを対応付けた結合ベクトルを得る。これをパラレルデータと呼ぶ。各フレームのパラレルデータを用い、対象音声のスペクトル特徴量と通常音声の対数F0 の静的・動的成分の結合確率密度関数を混合正規分布モデル（Gaussian Mixture Model; GMM）で表現する。GMM のパラメータはExpectation-Maximization アルゴリズムにより学習することができる。変換処理では、学習されたGMM を用いて、系列内変動を考慮した最尤系列変換法により、対象音声のスペクトル特徴量系列から通常音声のF0 パターンへと変換することができる。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Keigo Nakamura, Tomoki Toda, Hiroshi Saruwatari, Kiyohiro Shikano, "Speaking-aid systems using GMM-based voice conversion for electrolaryngeal speech," Speech Communication, vol. 54, no. 1, pp. 134-146, 2012.
【非特許文献２】Kou Tanaka, Tomoki Toda, Graham Neubig, Sakriani Sakti, Satoshi Nakamura, "A hybrid approach to electrolaryngeal speech enhancement based on noise reduction and statistical excitation generation," IEICE Transactions on Information and Systems, vol. E97-D, no. 6, pp. 1429-1437, Jun. 2014.
【非特許文献３】Kou Tanaka, Tomoki Toda, Graham Neubig, Sakriani Sakti, Satoshi Nakamura, "Direct F0 controlof an electrolarynx based on statistical excitation feature prediction and its evaluation through simulation," Proc. INTERSPEECH, pp. 31-35, Sep. 2014.
【非特許文献４】Hirokazu Kameoka, Jonathan Le Roux, Yasunori Ohishi, "A statistical model of speech F0 contours," ISCA Tutorial and Research Workshop on Statistical And Perceptual Audition (SAPA 2010), pp. 43-48, Sep. 2010.
【非特許文献５】Kota Yoshizato, Hirokazu Kameoka, Daisuke Saito, Shigeki Sagayama, "Hidden Markov convolutive mixture model for pitch contour analysis of speech," in Proc. The 13th Annual Conference of the International Speech Communication Association (Interspeech 2012), Sep. 2012.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
従来技術では、学習処理や変換処理において音声のF0 パターンの物理的な生成過程を考慮したモデルが用いられていなかったため、物理的に人間が発声しえないような不自然なF0 パターンを生成することが起こりえた。この問題に対し、F0パターンの物理的な生成過程を考慮した予測を行うことで、より自然なF0 パターンを生成できる可能性がある。
【０００６】
F0パターンは声帯に張力を与える甲状軟骨の運動によって生み出されており、非特許文献４、５ではその制御機構の確率モデルに基づき、フレーズ・アクセント指令と呼ぶ甲状軟骨の運動に関係するパラメータを推定する技術が提案されている。この技術では，フレーズ・アクセント指令の時系列の生成プロセスを隠れマルコフモデル(HMM) により表現した点がポイントの一つであり，HMM のトポロジーの設計や遷移確率の学習を通して、指令列に関する言語学的ないし先験的な知識をパラメータ推定に組み込むことが可能である。
【０００７】
本発明は、上記事情を鑑みてなされたものであり、F0 パターンの物理的な生成過程の制約を考慮しながらスペクトル特徴量系列に対応する最適なF0パターンを推定することができる基本周波数パターン予測装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の目的を達成するために本発明に係る基本周波数パターン予測装置は、学習サンプルのソース音声の時系列データとターゲット音声の時系列データとからなるパラレルデータを入力として、前記ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の時系列データから抽出される、各時刻の基本周波数とに基づいて、前記ソース音声の各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータを、前記第１確率分布と前記第２確率分布とを用いて表される規準を大きくするように学習する学習部と、予測対象のソース音声の時系列データを入力として、前記予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記学習部によって学習された前記第１確率分布のパラメータと前記第２確率分布のパラメータとに基づいて、前記規準を大きくするように、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数を予測する変換処理部と、を含んで構成されている。
【０００９】
本発明に係る基本周波数パターン予測方法は、学習部と、変換処理部とを含む基本周波数パターン予測装置における基本周波数パターン予測方法であって、前記学習部が、学習サンプルのソース音声の時系列データとターゲット音声の時系列データとからなるパラレルデータを入力として、前記ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の時系列データから抽出される、各時刻の基本周波数とに基づいて、前記ソース音声の各時刻のスペクトル特徴量ベクトルと、前記ターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータを、前記第１確率分布と前記第２確率分布とを用いて表される規準を大きくするように学習し、前記変換処理部が、予測対象のソース音声の時系列データを入力として、前記予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、前記学習部によって学習された前記第１確率分布のパラメータと前記第２確率分布のパラメータとに基づいて、前記規準を大きくするように、前記予測対象のソース音声に対応する前記ターゲット音声の各時刻の基本周波数を予測する。
【００１０】
本発明に係るプログラムは、上記の基本周波数パターン予測装置の各部としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００１１】
以上説明したように、本発明の基本周波数パターン予測装置、方法、及びプログラムによれば、ソース音声の各時刻のスペクトル特徴量ベクトルとターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータを、前記第１確率分布と前記第２確率分布とを用いて表される規準を大きくするように学習し、予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルから、第１確率分布と第２確率分布とを用いて表される規準を大きくするように、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数を予測することにより、F0 パターンの物理的な生成過程の制約を考慮しながらスペクトル特徴量系列に対応する最適なF0パターンを推定することができる、という効果が得られる。
【図面の簡単な説明】
【００１２】
【図１】ＨＭＭの状態遷移ネットワークの一例を説明するための図である。
【図２】ＨＭＭの状態遷移ネットワークの一例を説明するための図である。
【図３】ＨＭＭの状態遷移ネットワークの一例を説明するための図である。
【図４】状態の分割を説明するための図である。
【図５】本発明の第１の実施の形態に係る基本周波数パターン予測装置の構成を示す概略図である。
【図６】本発明の第１の実施の形態に係る基本周波数パターン予測装置の学習部の構成を示す概略図である。
【図７】本発明の第１の実施の形態に係る基本周波数パターン予測装置の変換処理部の構成を示す概略図である。
【図８】本発明の第１の実施の形態に係る基本周波数パターン予測装置における学習処理ルーチンの内容を示すフローチャートである。
【図９】本発明の第１の実施の形態に係る基本周波数パターン予測装置における基本周波数パターン予測処理ルーチンの内容を示すフローチャートである。
【図１０】本発明の第２の実施の形態に係る基本周波数パターン予測装置における学習処理ルーチンの内容を示すフローチャートである。
【図１１】本発明の第３の実施の形態に係る基本周波数パターン予測装置の学習部の構成を示す概略図である。
【図１２】本発明の第３の実施の形態に係る基本周波数パターン予測装置における学習処理ルーチンの内容を示すフローチャートである。
【図１３】実験結果を示す図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する技術は、音声の特徴量系列から基本周波数パターンを予測し、原音声の基本周波数パターンを予測した基本周波数パターンに置き換えることで音声の自然性を向上させることを目的とした音声処理技術である。
【００１４】
＜関連技術１：スペクトル特徴量系列からのF0 パターン予測方法＞
まず、スペクトル特徴量系列からのF0 パターン予測方法について説明する。
【００１５】
上記非特許文献１〜非特許文献３では、スペクトル特徴量系列からF0 パターンを予測する方法が提案されている。当該従来手法はスペクトル特徴量系列とF0 パターンの同時確率分布モデルのパラメータを学習する処理と学習した当該モデルを用いて所与のスペクトル特徴量系列からF0 パターンに変換する処理からなる。
【００１６】
＜学習処理＞
ソース音声（例えば電気音声）とターゲット音声（例えば自然音声）のパラレルデータが与えられているものとする。ソース音声のスペクトル特徴量ベクトルをc[k] とし，ターゲット音声の対数F0 とその動的成分（時間微分または時間差分）の結合ベクトル（F0 特徴量と呼ぶ。）をq[k] = (y[k];Δy[k])T とする。ここでk は離散時刻のインデックスである。音声特徴量c[k] としては例えば時刻k を中心とした前後数フレーム分のメルケプストラム（ベクトル）の系列を連結したベクトルに対し主成分分析により次元圧縮を行ったものを用いる。本手法ではc[k] とq[k] の同時確率分布を混合正規分布モデル（Gaussian Mixture Model; GMM）
【００１７】
【数１】
0006472005.tif 000002
【００１８】
でモデル化し、学習処理では所与のパラレルデータ｛c[k]; q[k]｝Kk=1 から当該GMM のパラメータ（各正規分布の重み，平均、分散共分散行列）を学習する。ただし、N(x;μ，Σ) は、xの確率密度関数が平均がμ、分散共分散行列がΣの正規分布で与えられることを意味する。
【００１９】
GMM のパラメータはExpectation-Maximization (EM) アルゴリズムにより推定することができる。学習したGMM パラメータをγとすると、条件付分布Ｐ(q[k]｜c[k]，γ) はスペクトル特徴量c[k] からF0 特徴量q[k] を予測するための分布と見ることができ、
【００２０】
【数２】
0006472005.tif 000003
【００２１】
のようにＰ(c[k]，q[k]｜γ) と同様にGMM で与えられる。ただし、
【００２２】
【数３】
0006472005.tif 000004
【００２３】
でありe(q｜c)mおよびD(q｜c)m は
【００２４】
【数４】
0006472005.tif 000005
【００２５】
で与えられる。
【００２６】
＜変換処理＞
変換処理では。所与のスペクトル特徴量系列
【００２７】
【数５】
0006472005.tif 000006
【００２８】
の下で，最尤のF0パターン
【００２９】
【数６】
0006472005.tif 000007
【００３０】
を以下の式（８）により求めることが目的である。
【００３１】
【数７】
0006472005.tif 000008
【００３２】
ただし、
【００３３】
【数８】
0006472005.tif 000009
【００３４】
であり、Ｗはｙとｑの関係を記述した変換行列（定数）である。ここで、P(ｑ｜ｃ，γ) は学習処理により学習したパラメータγで与えられるGMM であり、
【００３５】
【数９】
0006472005.tif 000010
【００３６】
で与えられる。ただし、ｍ= (ｍ1，...，ｍK) であり、ｍk は時刻k におけるGMMの成分インデックスを表す。ここで、P(q[k]｜c[k]，γ) が
【００３７】
【数１０】
0006472005.tif 000011
【００３８】
により近似できるとする。式(12)より、^mkは、データｃ[k]を生成したらしい確率が最も高い正規分布のインデックスを意味する。よって、式(9)よりP(ｑ｜ｃ，ｍ，＾γ）はすべてのkについて
【００３９】
【数１１】
0006472005.tif 000012
【００４０】
の積をとったもので与えられる。式(11)の近似によりP(ｑ｜ｃ，γ）は
【００４１】
【数１２】
0006472005.tif 000013
【００４２】
を連結したベクトルｅ(ｑ｜ｃ)を平均、
【００４３】
【数１３】
0006472005.tif 000014
【００４４】
を対角成分としたブロック対角行列Ｄ(q｜c) を分散共分散行列とした正規分布
【００４５】
【数１４】
0006472005.tif 000015
【００４６】
となる。これにｑ＝Ｗｙを代入し、ｙの分布となるように正規化すると、
【００４７】
【数１５】
0006472005.tif 000016
【００４８】
となるため、
【００４９】
【数１６】
0006472005.tif 000017
【００５０】
のようにｙの条件付き分布を得る。よって、式(8)の解は、
【００５１】
【数１７】
0006472005.tif 000018
【００５２】
となる。
【００５３】
＜関連技術２：F0パターン生成過程モデル＞
次に、F0パターン生成過程の確率モデルについて説明する。
【００５４】
音声のF0パターンの生成過程を記述したモデルに、藤崎の基本周波数（F0）パターン生成過程モデル（藤崎モデル）が知られている（非特許文献６）。
【００５５】
［非特許文献６］：H. Fujisaki, "In Vocal Physiology: Voice Production, Mechanisms and Functions," Raven Press, 1988.
【００５６】
藤崎モデルとは、甲状軟骨の運動によるF0 パターンの生成過程を説明した物理モデルである。藤崎モデルでは、甲状軟骨の二つの独立な運動（平行移動運動と回転運動）にそれぞれ伴う声帯の伸びの合計がF0の時間的変化をもたらすと解釈され、声帯の伸びとF0パターンの対数値y(t) が比例関係にあるという仮定に基づいてF0パターンがモデル化される。甲状軟骨の平行移動運動によって生じるF0パターンxp(t)をフレーズ成分、回転運動によって生じるF0 パターンxa(t) をアクセント成分と呼ぶ。藤崎モデルでは、音声のF0 パターンy(t) は、これらの成分に声帯の物理的制約によって決まるベースライン成分b を足し合わせたものとして、
【００５７】
【数１８】
0006472005.tif 000019
【００５８】
と表現される。これら二つの成分は二次の臨界制動系の出力と仮定され、
【００５９】
【数１９】
0006472005.tif 000020
【００６０】
と表される（＊ は時刻t に関する畳み込み演算）。ここでup(t) はフレーズ指令関数と呼ばれ、デルタ関数（フレーズ指令）の列からなり、ua(t) はアクセント指令関数と呼ばれ、矩形波（アクセント指令）の列からなる。これらの指令列には、発話の最初にはフレーズ指令が生起する、フレーズ指令は二連続で生起しない、異なる二つの指令は同時刻に生起しない、という制約条件がある。またαとβはそれぞれフレーズ制御機構、アクセント制御機構の固有角周波数であり、話者や発話内容によらず、おおよそα=3 rad/s、β=20 rad/s 程度であることが経験的に知られている。
【００６１】
＜関連技術３：F0パターン生成過程モデルパラメータ推定法＞
上述の藤崎モデルは以下のような確率モデルで記述することができる（非特許文献４、５、７参照）。
【００６２】
[非特許文献７]：石原達馬, 吉里幸太, 亀岡弘和, 齋藤大輔, 嵯峨山茂樹, \音声基本周波数の藤崎モデル指令列の統計的語彙モデル," 日本音響学会2013 年春季研究発表会講演論文集, 1-7-9, pp. 283-286, Mar. 2013.
【００６３】
まずフレーズ、アクセント指令関数のペアo[k] = (up[k], ua[k])T を出力するHMM を考える。ただし、k は離散時刻のインデックスを表す。状態出力分布は正規分布とし、各時刻の状態が与えられた下で
【００６４】
【数２０】
0006472005.tif 000021
【００６５】
により指令関数ペアo[k] が生成されるものとする。ここで｛sk｝Kk=1 はHMMの状態系列であり、平均ベクトルρ[k] はHMM の状態遷移の結果として定まる値である。具体的なHMM の構成の例を図１〜図３に示す。
【００６６】
図1に示すHMM の状態遷移ネットワークの例では、状態t = r0 においてμ(p)t [k] とμ(a)tはいずれも0 である。状態t = r0 からは状態p0 にのみ遷移することができ、状態t = p0においてμ(p)t [k] は非負値A(p)[k] をとり、μ(a)t は0 となる。状態t = p0 の次は状態r1 にのみ遷移することが許される。状態t = r0 同様、状態t = r1 においてμ(p)t [k] とμ(a)t はいずれも0 である。状態t = r1 からは状態a0，...，aN-1のいずれかにのみ遷移することができ、状態t = an においてμ(a)tは非負値A(a)n をとり、μ(p)t [k] は0 となる。状態t = anの次は状態r0 またはr1にのみ遷移することが許される。これよりμa[k] が矩形パルス列となることが保証される。
【００６７】
図2に示すHMMの状態遷移ネットワークの例では、状態t = r0 においてμ(p)t [k] とμ(a)tはいずれも0 である。状態t = r0 からは状態p0，...，pM-1のいずれかにのみ遷移することができ、状態t = pm においてμ(p)tは非負値A(p)m をとり、μ(a)t は0 となる。状態t = pm の次は状態r1にのみ遷移することが許される。状態t = r0 同様、状態t = r1 においてμ(p)t [k] とμ(a)t はいずれも0 である。状態t = r1 からは状態a0，...，aN-1のいずれかにのみ遷移することができ、状態t = an においてμ(a)tは非負値A(a)n をとり、μ(p)t は0 となる。状態t = anの次は状態r0 またはr1 にのみ遷移することが許される。これよりμa[k] が矩形パルス列となることが保証される。
【００６８】
図3に示すHMMの状態遷移ネットワークの例では、それぞれの終点と始点が連結された複数のLeft-to-Right 型HMM からなる。図１、２と同様、状態t = rlにおいてμ(p)t [k] とμ(a)t はいずれも0 である。また、状態t = pm においてμ(p)tは非負値A(p)m をとり、μ(a)t は0 となる。状態t = anにおいてμ(a)tは非負値A(a)nをとり、μ(p)
【００６９】
tは0となる。
【００７０】
ρ[k]は、図1 の例では、以下の式（２４）で表わされる。
【００７１】
【数２１】
0006472005.tif 000022
【００７２】
また、図2, 3 の例では、ρ[k]は、以下の式（２５）で表わされる。
【００７３】
【数２２】
0006472005.tif 000023
【００７４】
いずれの例においても、図４のようにそれぞれの状態を同じ出力分布をもついくつかの小状態に分割し、Left-to-Right 型の状態遷移経路を制約することで同一状態に停留する時間長の確率をパラメータ化することができる。図4 は状態an を分割した例である。例えばこの図のように全てのm≠0 に対してan,n′からan,n′+1 への状態遷移確率を1に設定することで，an,0からan,n′への遷移確率が状態an がn′ステップだけ持続する確率に対応し、アクセント指令の持続長の確率を設定したり学習したりできるようになる。同様にpm とrl も小状態に分割することで、フレーズ指令の持続長と指令間の間隔の長さの分布をパラメータ化することが可能になる。以後、状態集合を
【００７５】
【数２３】
0006472005.tif 000024
【００７６】
と表記する。上記のHMM の構成は次のように書ける。
【００７７】
【数２４】
0006472005.tif 000025
【００７８】
状態系列s = ｛sk｝Kk=1 が与えられたとき、このHMM はフレーズ指令関数up[k] とアクセント指令関数ua[k] のペアを出力する。式(18) と式(20) で示した通り、up[k] とua[k] にそれぞれgp[k] とga[k] が畳み込まれてフレーズ成分xp[k] とアクセント成分xa[k] が出力される。これを式で表すと、
【００７９】
【数２５】
0006472005.tif 000026
【００８０】
と書ける（＊ は離散時刻k に関する畳み込み演算）。このとき，F0 パターンx[k] は
【００８１】
【数２６】
0006472005.tif 000027
【００８２】
と三種類の成分の重ね合わせで書ける。ただしb は時刻によらないベースライン成分である。
【００８３】
また、実音声においては、いつも信頼のできるF0 の値が観測できるとは限らない。藤崎モデルのパラメータ推定を行うにあたっては、信頼のおける観測区間のF0 値のみを考慮に入れて、そうでない区間は無視することが望ましい。例えば音声の無声区間においては通常声帯の振動に伴う周期的な粗密波は観測されないので、仮に自動ピッチ抽出によって音声の無声区間から何らかの値がF0 の推定値として得られたとしても、その値を声帯から発せられる信号のF0 の値と見なすのは適当ではない。そこで、提案モデルに観測F0値の時刻k における不確かさの程度ｖ2n[k] を導入する。具体的には、観測F0 値y[k] を、真のF0 値x[k] とノイズ成分
【００８４】
【数２７】
0006472005.tif 000028
【００８５】
との重ね合わせで
【００８６】
【数２８】
0006472005.tif 000029
【００８７】
と表現することで、信頼のおける区間かどうかに関わらず全ての観測区間を統一的に扱える。
【００８８】
xn[k] を周辺化することで、出力値系列o = ｛o[k]｝Kk=1 が与えられたときのy = ｛y[k]｝Kk=1 の確率密度関数
【００８９】
【数２９】
0006472005.tif 000030
【００９０】
が得られる。状態系列s = ｛sk｝Kk=1 と指令の振幅を表すパラメータ
【００９１】
【数３０】
0006472005.tif 000031
【００９２】
および遷移確率行列φ = (φi,j)I×I が与えられたとき、出力値系列o は
【００９３】
【数３１】
0006472005.tif 000032
【００９４】
に従って生成される。また、P(s｜φ) は状態遷移確率の積として
【００９５】
【数３２】
0006472005.tif 000033
と書ける。ただし、
【００９６】
【数３３】
0006472005.tif 000034
【００９７】
は初期状態がs1である確率をあらわす。式(30)、(32) および式(33) よりP(y,o,s｜θ,φ) は
【００９８】
【数３４】
0006472005.tif 000035
【００９９】
と書ける。これをo に関して周辺化すると
【０１００】
【数３５】
0006472005.tif 000036
【０１０１】
が得られる。ただし、
【０１０２】
【数３６】
0006472005.tif 000037
0006472005.tif 000038
0006472005.tif 000039
【０１０３】
である。一方、s に関して周辺化すると
【０１０４】
【数３７】
0006472005.tif 000040
が得られる。ただし、Σs はあらゆる状態系列に関して和をとる操作を意味する。
【０１０５】
＜パラメータ推定アルゴリズム１＞
ｙとｏを完全データと見なすことで、式(35) を局所最大化するｓとθ をExpectation-Maximization アルゴリズムにより探索することができる。導出は省略するが、
【０１０６】
【数３８】
0006472005.tif 000041
【０１０７】
が大きくなるようにｓとθを更新するステップと、更新したｓとθを用いて
【０１０８】
【数３９】
0006472005.tif 000042
【０１０９】
とＲを
【０１１０】
【数４０】
0006472005.tif 000043
【０１１１】
により更新するステップを繰り返すことで式(35) を単調増加させることができる（詳細は、上記非特許文献４参照）。
【０１１２】
具体的には、以下の初期設定、Ｅステップ、及びＭステップが実行される。
【０１１３】
（初期設定）
ｓとθを初期設定する。
【０１１４】
（E ステップ）
フレーズ成分、アクセント成分、ベースライン成分の条件付き期待値
【０１１５】
【数４１】
0006472005.tif 000044
【０１１６】
と条件付き分散共分散行列Ｒを
【０１１７】
【数４２】
0006472005.tif 000045
により更新する。ただし、
【０１１８】
【数４３】
0006472005.tif 000046
0006472005.tif 000047
【０１１９】
である。また、Ｒにおける
【０１２０】
【数４４】
0006472005.tif 000048
【０１２１】
に対応するブロック対角成分を
【０１２２】
【数４５】
0006472005.tif 000049
【０１２３】
とする。
【０１２４】
すなわち、
【０１２５】
【数４６】
0006472005.tif 000050
【０１２６】
である（＊ は以後用いないブロック成分である）。
【０１２７】
（Ｍステップ）
Q(ｓ,θ) が最大となる状態系列ｓ= (ｓ1,...,ｓK) を探索する。γp とγaは対角行
【０１２８】
列なので、
【０１２９】
【数４７】
0006472005.tif 000051
0006472005.tif 000052
【０１３０】
はいずれも
【０１３１】
【数４８】
0006472005.tif 000053
【０１３２】
のようにｋごとの項の和の形で書ける。従って、Ｑ(ｓ，θ) はｓに依らない項を除けば
【０１３３】
【数４９】
0006472005.tif 000054
【０１３４】
と書ける。従って、Ｑ(ｓ,θ) を最大にする状態系列ｓ= (ｓ1,...,ｓK)はViterbiアルゴリズムにより求めることができる（詳細は上記非特許文献４参照）。ただし、[・]k,k は行列のk 行k 列の要素、[・]kはベクトルの第k 要素を表す。
【０１３５】
続いて、Ｑ(ｓ,θ)を最大にするようにθを更新する。Ｑ(ｓ,θ)を最大にするθは、Ｑ(ｓ,θ)の各変数に関する偏微分を0 と置くことにより得られる（詳細は上記非特許文献４参照）。
【０１３６】
また、推定された状態系列ｓから、状態遷移確率φが求められる。
【０１３７】
［第１の実施の形態］
＜本発明の実施の形態の概要＞
本発明の実施の形態の技術は、上述した関連技術１と同様、学習処理と変換処理からなるが、式(8) の代わりに関連技術1の確率分布と関連技術２の確率分布を"Expert"としたProduct-of-Experts （非特許文献９参照）を用いることにより，関連技術２のF0 パターン生成過程モデルにできるだけ即したF0 パターンをスペクトル特徴量から統計的に予測することを可能にする技術である。
【０１３８】
［非特許文献９］：G. E. Hinton, “Training Products of Experts by Minimizing Contrastive Divergence," Neural Computation, no. 14, no. 8, pp. 1771-1800, 2002.
【０１３９】
学習処理と変換処理では、後述するｃ、ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を共通の規準とする。
【０１４０】
学習処理ではパラレルデータの学習サンプル｛c[k]，q[k]｝Kk=1が与えられた下でｃ,ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)ができるだけ大きくなるようにγ、θ、φを学習する。本実施の形態では、γと、θ、φとを別々に学習する。例えば、学習サンプルのF0 パターン｛y[k]｝Kk=1が与えられた下で、上述した関連技術３のパラメータ推定アルゴリズム１を用いてＰ（ｙ,ｓ｜θ,φ）ができるだけ大きくなるようにθとφを学習する。そして、θ、φを固定し、Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)ができるだけ大きくなるようにγを学習する。
【０１４１】
なお、学習サンプルのフレーズ・アクセント指令データｏが与えられる下では、上述した関連技術３を用いてＰ(ｏ｜θ、φ) ができるだけ大きくなるようにθとφを学習するようにしてもよい。
【０１４２】
＜システム構成＞
次に、ソース音声のスペクトル特徴量系列から、ターゲット音声の基本周波数パターンを予測する基本周波数パターン予測装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【０１４３】
図５に示すように、本発明の第１の実施の形態に係る基本周波数パターン予測装置は、ＣＰＵと、ＲＡＭと、後述する学習処理ルーチン、及び基本周波数パターン予測処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【０１４４】
図５に示すように、基本周波数パターン予測装置１００は、入力部１０と、演算部２０と、出力部９０とを備えている。
【０１４５】
入力部１０は、学習サンプルのソース音声（例えば電気音声）の時系列データとターゲット音声（例えば自然音声）の時系列データとからなるパラレルデータを受け付ける。また、入力部１０は、予測対象のソース音声の時系列データを受け付ける。
【０１４６】
演算部２０は、学習部３０と、パラメータ記憶部４０と、変換処理部５０とを備えている。
【０１４７】
図６に示すように、学習部３０は、特徴量抽出部３２と、基本周波数系列抽出部３４と、第１モデルパラメータ学習部３６と、第２モデルパラメータ学習部３８とを備えている。
【０１４８】
特徴量抽出部３２は、入力部１０によって受け付けた学習サンプルのソース音声の時系列データから、ソース音声のスペクトグラム特徴量ベクトルc[k]を抽出する。ここでk は離散時刻のインデックスである。例えば、非特許文献１〜３と同様に、時刻k を中心とした前後数フレーム分のメルケプストラム（ベクトル）の系列を連結したベクトルに対し主成分分析により次元圧縮を行ったものをc[k] として用いる。
【０１４９】
基本周波数系列抽出部３４は、入力部１０によって受け付けた学習サンプルのターゲット音声の時系列データから、ターゲット音声の各時刻ｋにおける基本周波数y[k]を抽出し、y = (y[1]，...， y[K])Tとする。
【０１５０】
この基本周波数の抽出処理は、周知技術により実現でき、例えば、非特許文献８（H. Kameoka, "Statistical speech spectrum model incorporating all-pole vocal tract model and F0 contour generating process model," in Tech. Rep. IEICE, 2010, in Japanese.）に記載の手法を利用して、８ｍｓごとに基本周波数を抽出する。
【０１５１】
また、y とその動的成分（時間微分または時間差分）の結合ベクトル（F0 特徴量と呼ぶ。）をq[k] = (y[k]，Δy[k])T とする。
【０１５２】
以上より、｛c[k]，q[k]｝Kk=1 というデータが得られる。
【０１５３】
第２モデルパラメータ学習部３８は、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数y[k]に基づいて、各時刻ｋの基本周波数y[k]と、隠れマルコフモデルの各時刻の状態からなる状態系列ｓとの組み合わせの確率分布である第２確率分布のパラメータを学習する。
【０１５４】
具体的には、第２モデルパラメータ学習部３８は、上述した関連技術３のF0 パターン生成過程モデルパラメータ推定法のパラメータ推定アルゴリズム１に従って、F0パターン生成過程モデルのパラメータθ、φを学習する。
【０１５５】
もし学習サンプルのフレーズ指令系列とアクセント指令系列のペアのデータo = ｛ok｝Kk=1 が入手できるのであれば，o からθ、φを学習しても良い（HMM の通常の学習に相当）。学習したF0パターン生成過程モデルのパラメータを^θ、^φとする。
【０１５６】
第１モデルパラメータ学習部３６は、特徴量抽出部３２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、第２モデルパラメータ学習部３８によって学習されたパラメータθ、φとに基づいて、ｃ、ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)ができるだけ大きくなるように、ソース音声の各時刻のスペクトル特徴量ベクトルc[k]とターゲット音声の各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］との同時確率分布を表す混合正規分布である第１確率分布のパラメータγを学習する。
【０１５７】
＜学習・変換規準＞
ここで、学習処理と変換処理で共通となる規準について説明する。
【０１５８】
ｃ、ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)は以下の形で与えられる。
【０１５９】
【数５０】
0006472005.tif 000055
ただし、
0006472005.tif 000056
0006472005.tif 000057
【０１６０】
以上よりｃ、ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)は、GMMベースのF0パターン予測モデルとF0 パターン生成過程モデルの確率分布の積を正規化した分布と見なせ、ΛおよびΓは両モデルの寄与の大きさを表した行列である。いずれも任意の対角行列（定数）とする。他の変数については、上述した関連技術1と関連技術3 の記述と同様である。また、上記では分かりやすさのためGMM の成分インデックス系列とHMMの状態系列に依存する変数には上付き文字(m), (s) を付与している。
【０１６１】
＜学習・変換処理のための統一的アルゴリズム＞
次に、学習・変換処理のための統一的アルゴリズムについて説明する。
【０１６２】
学習処理も変換処理も同一規準を用いた最適化問題
【０１６３】
【数５１】
0006472005.tif 000058
【０１６４】
となり（固定する変数と推定する変数が異なるだけである）、ｑとｕを潜在変数（隠れ変数）と扱うことでExpectation-Maximization アルゴリズムを適用することができる。ｙ、ｃ、ｑ、ｕを完全データとすると，完全データに対する尤度関数は、以下のように書ける。
【０１６５】
【数５２】
0006472005.tif 000059
ただし、
0006472005.tif 000060
【０１６６】
式(92) を平方完成すると、Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)は、以下のように表される。
【０１６７】
【数５３】
0006472005.tif 000061
0006472005.tif 000062
【０１６８】
E ステップで、以下の式（１００）、（１０１）に従って期待値を算出する。
【０１６９】
【数５４】
0006472005.tif 000063
【０１７０】
Ｍステップで、
【０１７１】
【数５５】
0006472005.tif 000064
【０１７２】
が最大となるように各変数を更新することでｌｏｇＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくすることができるので、E ステップとM ステップを繰り返すことで当該最適化問題の局所最適解を得ることができる。
【０１７３】
以上説明した原理に従って、第１モデルパラメータ学習部３６は、第２モデルパラメータ学習部３６で学習されたパラメータθ、φを固定して、特徴量抽出部３２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、に基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、ｌｏｇＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、第１確率分布のＧＭＭのパラメータγを学習する。学習したGMM パラメータを^γとする。
【０１７４】
変換処理部５０は、予測対象のソース音声の時系列データを入力として、ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、第１モデルパラメータ学習部３６によって学習された第１確率分布のパラメータγと、第２モデルパラメータ学習部３８によって学習された第２確率分布のパラメータθ、φとに基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、学習処理と共通の規準であるＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、各時刻の基本周波数ｙと、各時刻ｋのスペクトル特徴量ベクトルを生成したらしい確率が最も高い正規分布のインデックスｍ、各時刻の状態からなる状態系列ｓを推定することにより、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数ｙを予測する。
【０１７５】
本実施の形態では、図７に示すように、変換処理部５０は、特徴量抽出部５２と、期待値算出部５４と、変数更新部５６と、収束判定部５８とを備えている。
【０１７６】
特徴量抽出部５２は、入力部１０によって受け付けた予測対象のソース音声の時系列データから、特徴量抽出部３２と同様に、ソース音声の各時刻ｋのスペクトグラム特徴量ベクトルc[k]を抽出する。
【０１７７】
期待値算出部５４は、第１モデルパラメータ学習部３６によって学習された第１確率分布のパラメータγと、第２モデルパラメータ学習部３８によって学習された第２確率分布のパラメータθ、φと、特徴量抽出部５２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、変数更新部５６によって前回更新された各時刻の基本周波数ｙ、各時刻ｋの正規分布のインデックスｍ、及び各時刻の状態からなる状態系列ｓとに基づいて、上記式（１００）、式（１０１）に従って、期待値を算出する。
【０１７８】
変数更新部５６は、期待値算出部５４によって算出された期待値に基づいて、上記式（１０２）が最大となるように、各時刻の基本周波数ｙと、各時刻ｋの正規分布のインデックスｍと、各時刻の状態からなる状態系列ｓとを更新する。
【０１７９】
収束判定部５８は、予め定められた収束判定条件を満たすまで、期待値算出部５４及び変数更新部５６による各処理を繰り返させる。収束判定条件としては、例えば、予め定められた繰り返し回数に到達することである。
【０１８０】
収束判定条件を満たしたときに、最終的に得られた各時刻ｋの基本周波数ｙ［ｋ］を、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数の予測結果として、出力部９０により出力する。
【０１８１】
＜基本周波数パターン予測装置の作用＞
次に、本実施の形態に係る基本周波数パターン予測装置１００の作用について説明する。まず、学習サンプルのソース音声の時系列データ及びターゲット音声の時系列データからなるパラレルデータが、基本周波数パターン予測装置１００に入力されると、基本周波数パターン予測装置１００において、図８に示す学習処理ルーチンが実行される。
【０１８２】
まず、ステップＳ１０１において、入力されたソース音声の時系列データを読み込み、
【０１８３】
各時刻ｋのスペクトル特徴量ベクトルｃ［ｋ］を抽出する。ステップＳ１０２において、入力されたターゲット音声の時系列データを読み込み、ターゲット音声の各時刻ｋにおける基本周波数y[k]を抽出し、また、基本周波数y[k]とその動的成分の結合ベクトルq[k]を抽出する。
【０１８４】
ステップＳ１０４では、状態系列ｓと、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータθとを初期設定する。
【０１８５】
そして、ステップＳ１０５において、上記式（４９）、式（５０）に従って、フレーズ成分、アクセント成分、ベースライン成分の条件付き期待値￣ｘと、条件付き分散共分散行列Ｒとを更新する。
【０１８６】
次のステップＳ１０６では、上記ステップＳ１０４で初期設定された、又は後述するステップＳ１０７で前回更新されたパラメータθと、上記ステップＳ１０５で更新されたフレーズ成分、アクセント成分、ベースライン成分の条件付き期待値￣ｘと、条件付き分散共分散行列Ｒとに基づいて、上記式（６３）式を用いて、Ｑ(ｓ,θ) を最大にする状態系列ｓ= (ｓ1,...,ｓK)をViterbiアルゴリズムにより求めて、状態系列ｓを更新する。
【０１８７】
ステップＳ１０７では、上記ステップＳ１０６で更新された状態系列ｓと、上記ステップＳ１０５で更新されたフレーズ成分、アクセント成分、ベースライン成分の条件付き期待値￣ｘと、条件付き分散共分散行列Ｒとに基づいて、Ｑ(ｓ,θ)の各変数に関する偏微分を0 と置くことにより、Ｑ(ｓ,θ) を最大にするパラメータθを求めて、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータθを更新する。
【０１８８】
ステップＳ１０８において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１０５へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１０９へ進む。
【０１８９】
ステップＳ１０９では、上記ステップＳ１０１で抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、上記ステップＳ１０２で抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、上記ステップＳ１０７で最終的に得られたパラメータθと、上記ステップＳ１０６で最終的に得られた状態系列ｓから求められる状態遷移確率φと、後述するステップＳ１１０で前回更新されたパラメータγとに基づいて、上記式（１００）、式（１０１）に従って、期待値を算出する。
【０１９０】
ステップＳ１１０では、上記ステップＳ１０９で算出された期待値に基づいて、上記式（１０２）が最大となるように、パラメータγを更新する。
【０１９１】
ステップＳ１１２において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１０９へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１１３において、上記ステップＳ１０７で最終的に得られたパラメータθと、上記ステップＳ１０６で最終的に得られた状態系列ｓから求められる状態遷移確率φと、上記ステップＳ１１０で最終的に得られたパラメータγとを、パラメータ記憶部４０に格納する。
【０１９２】
次に、予測対象のソース音声の時系列データが、基本周波数パターン予測装置１００に入力されると、基本周波数パターン予測装置１００において、図９に示す基本周波数パターン予測処理ルーチンが実行される。
【０１９３】
まず、ステップＳ１２１において、入力された予測対象のソース音声の時系列データを読み込み、各時刻ｋのスペクトル特徴量ベクトルｃ［ｋ］を抽出する。
【０１９４】
ステップＳ１２２において、パラメータ記憶部４０に記憶されたパラメータγと、上記ステップＳ１２１で抽出された各時刻のスペクトル特徴量ベクトルc[k]とに基づいて、上記式（１６）に従って、各時刻ｋの基本周波数ｙ［ｋ］を推定することにより、各時刻ｋの基本周波数ｙ［ｋ］を初期設定すると共に、各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］に初期値を設定する。また、各時刻ｋのスペクトル特徴量ベクトルc[k]を生成したらしい確率が最も高い正規分布のインデックス＾ｍkと、状態系列ｓとに、初期値を設定する。
【０１９５】
ステップＳ１２３では、上記ステップＳ１２１で抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、各時刻ｋの基本周波数ｙ［ｋ］、各時刻ｋの基本周波数ｙ［ｋ］から得られる各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、パラメータ記憶部４０に記憶されたパラメータθ、状態遷移確率φ、及びパラメータγとに基づいて、上記式（１００）、式（１０１）に従って、期待値を算出する。
【０１９６】
ステップＳ１２４では、上記ステップＳ１２３で算出された期待値に基づいて、上記式（１０２）が最大となるように、各時刻ｋの基本周波数ｙ［ｋ］、各時刻ｋの正規分布のインデックス＾ｍk、及び状態系列ｓを更新する。
【０１９７】
そして、ステップＳ１２５において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１２３へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１２６において、上記ステップＳ１２４で最終的に得られた各時刻ｋの基本周波数ｙ［ｋ］を、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数の予測結果として、出力部９０により出力し、基本周波数パターン予測処理ルーチンを終了する。
【０１９８】
以上説明したように、第１の実施の形態に係る基本周波数パターン予測装置によれば、基本周波数パターン生成過程をモデル化した第２確率分布のパラメータθ、φを学習し、ソース音声の各時刻のスペクトル特徴量ベクトルとターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータγを、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように学習し、予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルから、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数を予測することにより、F0 パターンの物理的な生成過程の制約を考慮しながらスペクトル特徴量系列に対応する最適なF0パターンを推定することができる
【０１９９】
［第２の実施の形態］
次に、本発明の第２の実施の形態に係る基本周波数パターン予測装置について説明する。なお、第１の実施の形態と同様の構成となる部分については同一符号を付して説明を省略する。
【０２００】
第２の実施の形態では、第２の確率分布が、Ｐ(ｙ，ｏ｜θ，φ)であり、パラメータθ、φを推定する方法が第１の実施の形態と異なっている。
【０２０１】
第２の実施の形態に係る基本周波数パターン予測装置の学習部３０の第２モデルパラメータ学習部３８によるパラメータを学習する原理について説明する。
【０２０２】
まず、関連技術３のF0パターン生成過程モデルパラメータ推定法のパラメータ推定アルゴリズムについて説明する。
【０２０３】
＜パラメータ推定アルゴリズム２＞
観測F0系列ｙが与えられたもとで、モデルパラメータθとｏの事後確率Ｐ(ｏ，θ｜ｙ)の局所最適解を求める反復アルゴリズムを以下に示す。状態系列ｓを隠れ変数とし、事後確率Ｐ(ｏ，θ｜ｙ) が、
【０２０４】
【数５６】
0006472005.tif 000065
【０２０５】
をｓについて周辺化することで得られる点に注意すると、Ｑ関数Ｑ(ｏ，θ，ｏ´，θ´)は
【０２０６】
【数５７】
0006472005.tif 000066
と置ける。ただし、
【０２０７】
【数５８】
0006472005.tif 000067
【０２０８】
は定数項を除いて等しいことを表す。また、ｇb[k] = δ[k] (クロネッカーのデルタ) である。よって、Ｐ(ｓk = ｔ｜ｙ，ｏ´，θ´)をForward-Backward アルゴリズムにより計算するステップ、ｏとθについてＱ(ｏ，θ，ｏ´，θ´)を増加させるステップを繰り返すことで、Ｐ(ｏ，θ｜ｙ) が局所最大となる解を得ることができる。ｏはフレーズ・アクセント指令系列のペアであるため、Ｑ(ｏ，θ，ｏ´，θ´)を増加させるステップにおいては、ｏの非負制約を考慮する必要がある。ｏの非負制約を満たしながらＱ(ｏ，θ，ｏ´，θ´)を増加させるような更新則は以下により導くことができる。まず、Ｑ(ｏ，θ，ｏ´，θ´)の下界はJensen の不等式より
【０２０９】
【数５９】
0006472005.tif 000068
【０２１０】
のように設計することができる。また、i，k，l は、
【０２１１】
【数６０】
0006472005.tif 000069
【０２１２】
を満たす任意の変数である。従ってQ関数の下界は、
【０２１３】
【数６１】
0006472005.tif 000070
0006472005.tif 000071
【０２１４】
と表される。この下界関数をλi,k,l≧0 に関して最大化するステップとｏ に関して最大化するステップを交互に繰り返せばＱ(ｏ，θ，ｏ´，θ´)を増加させることができる。いずれのステップの更新則も解析的に求めることができ、それぞれ
【０２１５】
【数６２】
0006472005.tif 000072
【０２１６】
で表される。以上の反復が収束したあと、続けてθを更新する。更新式は、図１の場合、
【０２１７】
【数６３】
0006472005.tif 000073
【０２１８】
図２、３の場合、
【０２１９】
【数６４】
0006472005.tif 000074
【０２２０】
である。これらの更新値をｏ´とθ´に代入したのちに、Ｐ(ｓk = ｔ｜ｙ，ｏ´，θ´)の更新を再度行い、以後同様の処理を繰り返すことで事後確率Ｐ(ｏ，θ｜ｙ)を単調増加させることができる。
【０２２１】
以上の反復アルゴリズムが収束した後、上述したパラメータ推定アルゴリズム１のViterbi アルゴリズムにより求まる最適なｓを状態系列の推定結果とする。
【０２２２】
また、推定された状態系列ｓから、状態遷移確率φが求められる。
【０２２３】
以上説明した原理に従って、第２モデルパラメータ学習部３８は、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数y[k]に基づいて、各時刻ｋの基本周波数y[k]と、各時刻ｋにおける甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令ｕp［ｋ］及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令ｕa［ｋ］のペアからなる指令関数ｏ［ｋ］との組み合わせの確率分布である第２確率分布のパラメータθ、φを学習する。
【０２２４】
第１モデルパラメータ学習部３６は、上記第１の実施の形態と同様に、第２モデルパラメータ学習部３６で学習されたパラメータθ、φを固定して、特徴量抽出部３２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、に基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、ｌｏｇＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、第１確率分布のＧＭＭのパラメータγを学習する。学習したGMM パラメータを^γとする。
【０２２５】
変換処理部５０は、第１の実施の形態と同様に、予測対象のソース音声の時系列データを入力として、ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、第１モデルパラメータ学習部３６によって学習された第１確率分布のパラメータγと、第２モデルパラメータ学習部３８によって学習された第２確率分布のパラメータθ、φとに基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、学習処理と共通の規準であるＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、各時刻の基本周波数ｙと、各時刻ｋのスペクトル特徴量ベクトルを生成したらしい確率が最も高い正規分布のインデックスｍ、各時刻の状態からなる状態系列ｓを推定することにより、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数ｙを予測する。
【０２２６】
＜基本周波数パターン予測装置の作用＞
【０２２７】
次に、第２の実施の形態に係る基本周波数パターン予測装置の作用について説明する。なお、第１の実施の形態と同様の処理については、同一符号を付して詳細な説明を省略する。
【０２２８】
まず、学習サンプルのソース音声の時系列データ及びターゲット音声の時系列データからなるパラレルデータが、基本周波数パターン予測装置に入力されると、基本周波数パターン予測装置において、図１０に示す学習処理ルーチンが実行される。
【０２２９】
まず、ステップＳ１０１において、入力されたソース音声の時系列データを読み込み、
【０２３０】
各時刻ｋのスペクトル特徴量ベクトルｃ［ｋ］を抽出する。ステップＳ１０２において、入力されたターゲット音声の時系列データを読み込み、ターゲット音声の各時刻ｋにおける基本周波数y[k]を抽出し、また、基本周波数y[k]とその動的成分の結合ベクトルq[k]を抽出する。
【０２３１】
ステップＳ２００では、指令系列ｏと、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータθとを初期設定する。また、ターゲット音声の時系列データに基づいて、有声区間、無声区間を特定し、各時刻ｋの基本周波数の不確かさの程度vn2 [ｋ]を推定する。
【０２３２】
そして、ステップＳ２０１において、上記ステップＳ２００で設定された指令系列oの初期値、または後述するステップＳ２０３で前回更新された指令系列oに基づいて、（k,t）の全ての組み合わせについて、事後確率Ｐ(sk=t|y,o′,θ′)を更新する。
【０２３３】
ステップＳ２０２では、上記ステップＳ２００で設定された指令系列oの初期値、または後述するステップＳ２０３で前回更新された指令系列oに基づいて、（ｋ、ｌ）の全ての組み合わせについて、上記の式（７１）に従って、補助変数λp,k,l、λa,k,l、λb,k,lを算出して更新する。
【０２３４】
次のステップＳ２０３では、上記ステップＳ１０２で抽出されたされた基本周波数系列ｙと、上記ステップＳ２００で算出された各時刻ｋの不確かさの程度vn2 [ｋ]と、上記ステップＳ２０１で更新された事後確率Ｐ(sk=t|y,o′,θ′)と、上記ステップＳ２０２で更新された補助変数λp,k,l、λa,k,l、λb,k,lとに基づいて、上記式（７２）に従って、非負値である各時刻ｌのフレーズ指令ｕp［ｌ］及びアクセント指令ｕa［ｌ］からなる指令系列oとベース成分ｕbとを更新する。
【０２３５】
次のステップＳ２０４では、収束条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、収束条件を満足していないと判断して、上記ステップＳ２０２へ戻る。一方、繰り返し回数ｓがＳに到達した場合には、収束条件を満足したと判断し、ステップＳ２０５で、上記ステップＳ２０３で更新された各時刻ｋのフレーズ指令ｕp［ｋ］及びアクセント指令ｕa［ｋ］と、上記ステップＳ２０１で更新された事後確率Ｐ(sk=t|y,o′,θ′)とに基づいて、上記式（７３）、式（７４）、又は式（７５）、式（７６）に従って、各時刻ｋのフレーズ指令の振幅A(p)[k]、及び各位置ｎのアクセント指令の振幅Aa(a)を更新することにより、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータθを更新する。
【０２３６】
ステップＳ２０６において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ２０１へ戻る。一方、収束判定条件を満たした場合には、ステップＳ２０７において、上記ステップＳ２０３で最終的に更新された指令系列ｏに基づいて、Viterbi アルゴリズムにより、状態系列ｓを推定する。また、推定された状態系列ｓから、状態遷移確率φを求める。
【０２３７】
そして、ステップＳ１０９では、上記ステップＳ１０１で抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、上記ステップＳ１０２で抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、上記ステップＳ２０５で最終的に得られたパラメータθと、上記ステップＳ２０７で得られた状態遷移確率φと、後述するステップＳ１１０で前回更新されたパラメータγとに基づいて、上記式（１００）、式（１０１）に従って、期待値を算出する。
【０２３８】
ステップＳ１１０では、上記ステップＳ１０９で算出された期待値に基づいて、上記式（１０２）が最大となるように、パラメータγを更新する。
【０２３９】
ステップＳ１１２において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ１０９へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１１３において、上記ステップＳ２０５で最終的に得られたパラメータθと、上記ステップＳ２０７で得られた状態遷移確率φと、上記ステップＳ１１０で最終的に得られたパラメータγとを、パラメータ記憶部４０に格納する。
【０２４０】
次に、予測対象のソース音声の時系列データが、基本周波数パターン予測装置１００に入力されると、基本周波数パターン予測装置１００において、上記図９に示す基本周波数パターン予測処理ルーチンが実行され、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数の予測結果が、出力部９０により出力される。
【０２４１】
以上説明したように、第２の実施の形態に係る基本周波数パターン予測装置によれば、基本周波数パターン生成過程をモデル化した第２確率分布のパラメータθ、φを学習し、ソース音声の各時刻のスペクトル特徴量ベクトルとターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータγを、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように学習し、予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルから、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数を予測することにより、F0 パターンの物理的な生成過程の制約を考慮しながらスペクトル特徴量系列に対応する最適なF0パターンを推定することができる
【０２４２】
［第３の実施の形態］
次に、本発明の第３の実施の形態に係る基本周波数パターン予測装置について説明する。なお、第１の実施の形態と同様の構成となる部分については同一符号を付して説明を省略する。
【０２４３】
第３の実施の形態では、パラメータγ、θ、φを推定する方法が第１の実施の形態と異なっている。
【０２４４】
第３の実施の形態では、学習サンプルのF0 パターン｛y[k]｝Kk=1が与えられた下で、Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)ができるだけ大きくなるように、パラメータγ、θ、φを同時に学習する。
【０２４５】
図１１に示すように、第３の実施の形態に係る基本周波数パターン予測装置の学習部３０は、特徴量抽出部３２と、基本周波数系列抽出部３４と、モデルパラメータ学習部３３６とを備えている。
【０２４６】
モデルパラメータ学習部３３６は、特徴量抽出部３２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］とに基づいて、ｃ、ｙの結合分布Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)ができるだけ大きくなるように、ソース音声の各時刻のスペクトル特徴量ベクトルc[k]とターゲット音声の各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］との同時確率分布を表す混合正規分布である第１確率分布のパラメータγ、及び各時刻ｋの基本周波数y[k]と、隠れマルコフモデルの各時刻の状態からなる状態系列ｓとの組み合わせの確率分布である第２確率分布のパラメータθ、φを学習する。
【０２４７】
具体的には、特徴量抽出部３２によって抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、基本周波数系列抽出部３４によって抽出された各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］とに基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、ｌｏｇＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、第１確率分布のＧＭＭのパラメータγ、及び第２確率分布のパラメータθ、φを学習する。
【０２４８】
変換処理部５０は、第１の実施の形態と同様に、予測対象のソース音声の時系列データを入力として、ソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルと、モデルパラメータ学習部３３６によって学習された第１確率分布のパラメータγ、及び第２確率分布のパラメータθ、φとに基づいて、ＥＭ（Expectation-Maximization）アルゴリズムにより、E ステップとM ステップを繰り返すことで、学習処理と共通の規準であるＰ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、各時刻の基本周波数ｙと、各時刻ｋのスペクトル特徴量ベクトルを生成したらしい確率が最も高い正規分布のインデックスｍ、各時刻の状態からなる状態系列ｓを推定することにより、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数ｙを予測する。
【０２４９】
＜基本周波数パターン予測装置の作用＞
次に、第３の実施の形態に係る基本周波数パターン予測装置の作用について説明する。なお、第１の実施の形態と同様の処理については、同一符号を付して詳細な説明を省略する。
【０２５０】
まず、学習サンプルのソース音声の時系列データ及びターゲット音声の時系列データからなるパラレルデータが、基本周波数パターン予測装置に入力されると、基本周波数パターン予測装置において、図１２に示す学習処理ルーチンが実行される。
【０２５１】
まず、ステップＳ１０１において、入力されたソース音声の時系列データを読み込み、各時刻ｋのスペクトル特徴量ベクトルｃ［ｋ］を抽出する。ステップＳ１０２において、入力されたターゲット音声の時系列データを読み込み、ターゲット音声の各時刻ｋにおける基本周波数y[k]を抽出し、また、基本周波数y[k]とその動的成分の結合ベクトルq[k]を抽出する。
【０２５２】
ステップＳ３００では、状態系列ｓと、状態系列ｓから求められるパラメータφと、各時刻における状態に応じたフレーズ指令の振幅及び各アクセント指令の振幅を表すパラメータθと、パラメータγとを初期設定する。
【０２５３】
そして、ステップＳ３０１において、上記ステップＳ１０１で抽出された各時刻ｋのスペクトル特徴量ベクトルc[k]と、上記ステップＳ１０２で抽出された各時刻ｋの基本周波数ｙ［ｋ］、各時刻ｋの基本周波数ｙ［ｋ］から得られる各時刻ｋの基本周波数の結合ベクトルｑ［ｋ］と、初期設定された、又は後述するステップＳ３０２で前回更新されたパラメータθ、状態遷移確率φ、及びパラメータγとに基づいて、上記式（１００）、式（１０１）に従って、期待値を算出する。
【０２５４】
ステップＳ３０２では、上記ステップＳ３０１で算出された期待値に基づいて、上記式（１０２）が最大となるように、パラメータθ、状態系列ｓ、状態系列ｓから求められる状態遷移確率φ、パラメータγを更新する。
【０２５５】
そして、ステップＳ３０３において、予め定められた収束判定条件を満たしたか否かを判定し、収束判定条件を満たしていない場合には、上記ステップＳ３０１へ戻る。一方、収束判定条件を満たした場合には、ステップＳ１１３において、上記ステップＳ３０２で最終的に得られたパラメータθと、状態系列ｓから求められる状態遷移確率φと、パラメータγとを、パラメータ記憶部４０に格納する。
【０２５６】
次に、予測対象のソース音声の時系列データが、基本周波数パターン予測装置１００に入力されると、基本周波数パターン予測装置１００において、上記図９に示す基本周波数パターン予測処理ルーチンが実行され、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数の予測結果が、出力部９０により出力される。
【０２５７】
以上説明したように、第３の実施の形態に係る基本周波数パターン予測装置によれば、ソース音声の各時刻のスペクトル特徴量ベクトルとターゲット音声の各時刻の基本周波数との間の関係をモデル化した第１確率分布のパラメータγ、及び基本周波数パターン生成過程をモデル化した第２確率分布のパラメータθ、φを、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように学習し、予測対象のソース音声の時系列データから抽出される各時刻のスペクトル特徴量ベクトルから、第１確率分布と第２確率分布とを用いて表される規準Ｐ(ｃ,ｙ,ｍ,ｓ｜γ,θ,φ)を大きくするように、予測対象のソース音声に対応するターゲット音声の各時刻の基本周波数を予測することにより、F0 パターンの物理的な生成過程の制約を考慮しながらスペクトル特徴量系列に対応する最適なF0パターンを推定することができる
【０２５８】
＜実験＞
音声信号からスペクトル特徴量系列とF0 パターンおよびフレーズ・アクセント指令を抽出し、スペクトル特徴量系列とフレーズ・アクセント指令系列のペアデータを用いて学習処理により上記のモデルパラメータ（GMM のパラメータ）を学習したのちに、変換処理によりスペクトル特徴量系列をフレーズ・アセント指令系列に変換する実験を行い、変換されたフレーズ・アセント指令系列が元のF0 パターンをどの程度復元できているかを確認した。図１３にその結果の例を示す。点線が、音声信号から推定されたＦ0パターンであり、破線が、スペクトル特徴量系列から変換されたフレーズ・アセント指令系列から得られたF0 パターンである。スペクトル特徴量にはF0の情報が多く含まれていないにもかかわらず概ね元のF0 パターンを復元できていることが確認できた。
【０２５９】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０２６０】
例えば、上述の基本周波数パターン予測装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０２６１】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０２６２】
１０   入力部
２０   演算部
３０   学習部
３２   特徴量抽出部
３４   基本周波数系列抽出部
３６   第１モデルパラメータ学習部
３８   第２モデルパラメータ学習部
４０   パラメータ記憶部
５０   変換処理部
５２   特徴量抽出部
５４   期待値算出部
５６   変数更新部
５８   収束判定部
９０   出力部
１００ 基本周波数パターン予測装置
３３６ モデルパラメータ学習部
【図１】
0006472005.tif 000076
【図２】
0006472005.tif 000077
【図３】
0006472005.tif 000078
【図４】
0006472005.tif 000079
【図５】
0006472005.tif 000080
【図６】
0006472005.tif 000081
【図７】
0006472005.tif 000082
【図８】
0006472005.tif 000083
【図９】
0006472005.tif 000084
【図１０】
0006472005.tif 000085
【図１１】
0006472005.tif 000086
【図１２】
0006472005.tif 000087
【図１３】
0006472005.tif 000088