JP5188300B2 - 基本周波数軌跡モデルパラメータ抽出装置、基本周波数軌跡モデルパラメータ抽出方法、プログラム及び記録媒体 - Google Patents

Description

本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定するための基本周波数軌跡モデルパラメータ抽出装置、基本周波数軌跡モデルパラメータ抽出方法、プログラム及び記録媒体に関する。

歌唱音響信号から推定される基本周波数（Ｆ０）軌跡には、歌唱者が歌おうとする音高目標値の時系列と歌唱力・歌唱スタイル・個人性・感情に基づく様々な動的変動が観測される。歌声は、多くのジャンルの音楽を特徴付ける重要な要素の１つであり、この歌声のＦ０軌跡に着目した様々な研究が現在盛んに行われている。特に、歌声から楽曲を検索するハミング検索では、歌唱された歌声のＦ０軌跡から、歌唱者の意図する音高列を正しく推定して、楽曲データベースの旋律と照合する必要がある。

従来、Ｆ０軌跡そのものをＤＰマッチングによって照合する方法が提案されている（例えば、非特許文献１、非特許文献２参照）。しかしながら、これら従来技術では、歌声の動的変動の影響を受けて検索性能が低下するという問題があった。

ところで、オーバーシュートや、ビブラートのようなＦ０軌跡の動的変動は、歌声知覚に影響を与え、歌声の自然性を保つためには必要不可欠な成分であることが知られている。そこで、これらの動的変動を制御する２次系Ｆ０制御モデルが提案され、自然性、かつ、明瞭性のある歌声合成技術が実現されている（例えば、非特許文献３参照）。つまり、楽譜に相当する階段状の信号に２次系のインパルス応答を畳み込むことが、歌声のＦ０軌跡の生成モデルとして有効であることが示された。
橋口博樹、西村拓一、張建新、滝田順子、岡隆一、"モデル依存傾斜制限型の連続ＤＰを用いた鼻歌入力による楽曲信号のスポッティング検索,"電子情報通信学会論文誌D-II, Vol. J84-D-II, No. 12, pp. 2479-2488, 2001. Adams, N. H. et al., "Time Series Alignment for Music Information Retrieval," In Proc. ISMIR 2004, 2004. Saitou, T., Unoki, M. and Akagi, M., "Development of an F0 control model based on F0 dynamic characteristics for singing-voice synthesis," Speech Communication, Vol. 46, pp. 405-417, 2005.

しかしながら、上述した非特許文献３による従来技術では、制御パラメータが手作業あるいは規則に基づいて決定されるものであり、Ｆ０軌跡から自動推定する方法は確立されていない。つまり、音声認識や音声合成で提案される学習アルゴリズムの枠組みが、上述した歌声合成技術では確立されていない。すなわち、入力となる階段状の信号および２次系の制御パラメータがいずれも未知の下で、観測されるＦ０軌跡だけから、それらを推定することは不良設定問題であり、その解法は提案されていない。

本発明は、このような事情を考慮してなされたものであり、その目的は、動的変動そのものを適切にモデル化し、Ｆ０軌跡のみから旋律を構成する音高列を正しく推定することができる基本周波数軌跡モデルパラメータ抽出装置、基本周波数軌跡モデルパラメータ抽出方法、プログラム及び記録媒体を提供することにある。

上述した課題を解決するために、本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置であって、入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出手段と、前記基本周波数抽出手段により抽出された観測基本周波数時系列をフレーム分割するフレーム分割手段と、前記基本周波数抽出手段により抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成手段と、前記初期音高目標値生成手段により生成された初期音高目標値時系列と前記フレーム分割手段によりフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新手段と、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新手段と、前記逆フィルタ出力値系列更新手段により生成された特性パラメータと前記音高目標値更新手段により生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新手段と前記音高目標値更新手段による処理を再度行なわせる収束判定手段と、前記収束判定手段により所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力手段とを備え、前記逆フィルタ出力値系列更新手段は、前記観測基本周波数時系列ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ 、前記音高目標値時系列ｍ＝（ｍ _ｑ１ ，…，ｍ _ｑＮ ） ^Ｔ とし、定数行列Ｂ、Ａを、下記式（２５）とし、Ｂ ^Ｔ はＢの転置行列を表わすものとし、ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ は、ｙ _１ からｙ _Ｎ を縦に並べたベクトルを表わし、行列のｉ，ｊ成分がｘ _ｉ，ｊ から構成される行列Ｘ＝（Ｘ _ｉ，ｊ ）と書くとき、Ｔｒ ^（ｆ） （Ｘ）を、下記式（２６）で表わし、Ｄ _ｆ を下記式（２７）で表した場合に下記式（２４）を解いて、α ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） を求め、求まったα ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） と前記定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ ^（ｆ） を、下記式（２８）により算出して出力し、これにより、前記逆フィルタ出力値系列を、下記式（２９）により求めて出力し、前記音高目標値更新手段は、前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、前記逆フィルタ出力値系列の成分である下記式（３６）と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である下記式（４２）を用いて、下記式（３７）、下記式（３８）、下記式（３９）、下記式（４０）、下記式（４１）の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新することを特徴とする基本周波数軌跡モデルパラメータ抽出装置である。

本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置であって、入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出手段と、前記基本周波数抽出手段により抽出された観測基本周波数時系列をフレーム分割するフレーム分割手段と、前記基本周波数抽出手段により抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成手段と、前記初期音高目標値生成手段により生成された初期音高目標値時系列と前記フレーム分割手段によりフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新手段と、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新手段と、前記逆フィルタ出力値系列更新手段により生成された特性パラメータと前記音高目標値更新手段により生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新手段と前記音高目標値更新手段による処理を再度行なわせる収束判定手段と、前記収束判定手段により所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力手段とを備え、前記逆フィルタ出力値系列更新手段は、行列ｍ ^（ｆ） を下記式（３２）で表わし、行列Ｕを下記式（３３）で表した場合に、下記式（３０）により計算されるベクトルｗ＝（ｗ _１ ，…，ｗ _Ｍ−１ ） ^Ｔ の要素ｗ _１ ，…，ｗ _Ｎ−１ を、下記式（３１）に代入して特徴パラメータ行列Ｗ ^（ｆ） を出力し、下記式（３４）により前記逆フィルタ出力値系列を求めて出力し、前記音高目標値更新手段は、前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、前記逆フィルタ出力値系列の成分である下記式（３６）と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である下記式（４２）を用いて、下記式（３７）、下記式（３８）、下記式（３９）、下記式（４０）、下記式（４１）、の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新することを特徴とする基本周波数軌跡モデルパラメータ抽出装置である。

本発明は、上記の発明において、前記逆フィルタ出力値系列更新手段は、Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、下記式（９８）のＪ_１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を求めることにより前記特徴パラメータ行列Ｗを決定することを特徴とする。

また、上述した課題を解決するために、本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出方法であって、入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップとを含み、前記逆フィルタ出力値系列更新ステップにおいて、前記観測基本周波数時系列ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ 、前記音高目標値時系列ｍ＝（ｍ _ｑ１ ，…，ｍ _ｑＮ ） ^Ｔ とし、定数行列Ｂ、Ａを、下記式（２５）とし、Ｂ ^Ｔ はＢの転置行列を表わすものとし、ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ は、ｙ _１ からｙ _Ｎ を縦に並べたベクトルを表わし、行列のｉ，ｊ成分がｘ _ｉ，ｊ から構成される行列Ｘ＝（Ｘ _ｉ，ｊ ）と書くとき、Ｔｒ ^（ｆ） （Ｘ）を、下記式（２６）で表わし、Ｄ _ｆ を下記式（２７）で表した場合に、下記式（２４）を解いて、α ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） を求め、求まったα ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） と前記定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ ^（ｆ） を、下記式（２８）により算出して出力し、これにより、前記逆フィルタ出力値系列を、下記式（２９）により求めて出力し、前記音高目標値更新ステップにおいて、前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、前記逆フィルタ出力値系列の成分である下記式（３６）と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である下記式（４２）を用いて、下記式（３７）、下記式（３８）、下記式（３９）、下記式（４０）、下記式（４１）の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新することを特徴とする基本周波数軌跡モデルパラメータ抽出方法である。

また、上述した課題を解決するために、本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出方法であって、入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップとを含み、前記逆フィルタ出力値系列更新ステップにおいて、行列ｍ ^（ｆ） を下記式（３２）で表わし、行列Ｕを下記式（３３）で表した場合に、下記式（３０）により計算されるベクトルｗ＝（ｗ _１ ，…，ｗ _Ｍ−１ ） ^Ｔ の要素ｗ _１ ，…，ｗ _Ｎ−１ を、下記式（３１）に代入して特徴パラメータ行列Ｗ ^（ｆ） を出力し、下記式（３４）により前記逆フィルタ出力値系列を求めて出力し、前記音高目標値更新ステップにおいて、前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、前記逆フィルタ出力値系列の成分である下記式（３６）と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である下記式（４２）を用いて、下記式（３７）、下記式（３８）、下記式（３９）、下記式（４０）、下記式（４１）の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新することを特徴とする基本周波数軌跡モデルパラメータ抽出方法である。

本発明は、上記の発明において、前記逆フィルタ出力値系列更新ステップは、Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、下記式（９８）のＪ_１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を求めることにより前記特徴パラメータ行列Ｗを決定することを特徴とする。

また、上述した課題を解決するために、本発明は、上述の基本周波数軌跡モデルパラメータ抽出方法の各ステップをコンピュータにより実行させるためのプログラムである。

また、上述した課題を解決するために、本発明は、上述のプログラムを記録したコンピュータ読み取り可能な記録媒体である。

この発明によれば、入力される音響信号から観測基本周波数時系列を抽出し、抽出された観測基本周波数時系列をフレーム分割し、抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成し、生成された初期音高目標値時系列とフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成し、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成し、生成された特性パラメータと生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、逆フィルタ出力値系列更新と前記音高目標値更新による処理を再度行なわせ、所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する。したがって、動的変動そのものを適切にモデル化することができ、Ｆ０軌跡のみから旋律を構成する音高列を正しく推定することができるという利点が得られる。つまり、Ｆ０軌跡に含まれる動的変動そのものを適切にモデル化できるため、高精度なハミング検索や高品質な歌声合成が可能となる。

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

Ａ．原理
まず、本発明の原理について説明する。なお、以下の説明において、記号Ｎはサンプル数を示し、記号≡は、左辺の式を右辺の式で定義することを意味し、記号∝は、左辺が右辺と比例することを意味する。また、確率（遷移確率）を記号^＊Pで表わすこととする。

Ａ−１．基本周波数抽出過程
音響信号を入力として観測基本周波数時系列を出力する。出力した観測基本周波数時系列を次式（１）と表す。

以後、ｙ_ｎのことをサンプルｎにおける観測基本周波数と呼ぶことにする。また、サンプルｎは、時刻ｎ△に対応するものとする。ここで、△をサンプリング周期と呼ぶ。

Ａ−２．フレーム分割過程
基本周波数抽出過程１において出力された観測基本周波数時系列を適当なフレーム長およびフレームシフト長により複数のフレームに分ける。フレーム番号をｆとし、ｆ番目のフレームにおける観測基本周波数の集合を、次式（２）と置く。

また、集合У^（ｆ）の要素を縦に並べたベクトルを、次式（３）と置く。

但し、Ｔはフレームシフト長、Ｍはフレーム長である。また、Fはフレーム数である。例えば、フレーム長Ｍが４点、フレームシフト長Ｔが２点の場合、次式（４）に示すようになる。

Ａ−３．初期音高目標値生成過程
基本周波数抽出過程１において出力された観測基本周波数時系列ｙ_１，…，ｙ_Ｎを入力として、次式（５）で示される、音高目標値時系列の初期値を生成する。

但し、ｑ_ｎを「状態」と呼び、各時刻ｎおいて１，２，…，Ｉ（Ｉは正の整数）の中のいずれかの値をとる。例えば、Ｉ＝３とすると、次式（６）のようになる。

ｍ_ｉは、ｉに対応した実数値を表し、ｍ_１，…，ｍ_ｉ，…，ｍ_Ｉを「音高目標値集合」と呼ぶ。従って、ｍ_ｑ１，…，ｍ_ｑｎは、状態系列ｑ_１，…，ｑ_Ｎに対応した実数値列を表し、これを「音高目標値時系列」と呼ぶ。例えば、ｍ_１＝５０、ｍ_２＝１５０、ｍ_３＝１２５とし、ｑ_ｎを数式（６）とすると、ｍ_ｑ１，…，ｍ_ｑｎは、次式（７）となる。

音高目標値時系列の初期値ｍ_ｑ１，ｍ_ｑ２，…，ｍ_ｑＮは、具体的には以下のような２つのステップにより求める。

Ａ−３−１．音高目標値時系列生成過程（状態系列ｑ_１，…，ｑ_Ｎの決定）
ｍ_１，…，ｍ_Ｉを適当に設定し（例えば、１２平均律音階に対応する周波数値）、動的計画法に基づいて、数式（８）〜（１１）に示すように、以下のような再帰計算を行なう。
Ｓａ１．初期化：

Ｓａ２．再帰計算：

Ｓａ３．終了：

Ｓａ４．バックトラック：

Ａ−３−２．音高目標値集合生成過程（音高目標値集合ｍ_ｌ，…，ｍ_Ｉの決定）
ステップＳａ１により決まったｑ１，…，ｑＮをそれぞれ次式（１２）と置き、次式（１３）により、数式（１４）を求める。

以上より求まった、数式（１２）、（１４）を用いて、次式（１５）で示す、音高目標値時系列の初期値が求まる。

但し、P_ｉ，ｊ（ｉ，ｊ∈｛１，…，Ｉ｝）は、予め設定しておく定数であり、「状態ｉから状態ｊへの行きやすさ（難しさ）」を表す。例えば、Ｐ_ｉ，ｊを大きめ、Ｐ_ｉ，ｊ（ｉ≠ｊ）を小さめにとると、ｍ_ｑ１，…，ｍ_ｑｎは、移り変わりの少ない安定した階段状の系列として推定されやすくなる。逆に、Ｐ_ｉ，ｊを小さめ、Ｐ_ｉ，ｊ（ｉ≠ｊ）を大きめにとると、ｍ_ｑ１，…，ｍ_ｑｎは、同じ値に長く留まろうとしないようになり、移り変わりの激しい系列として推定されやすくなる。このように、定数Ｐ_ｉ，ｊ（ｉ，ｊ∈｛１，…，Ｉ｝）は、音高目標値系列の移り変わりの激しさを調節するための定数である。

Ａ−４．逆フィルタ出力値更新過程
音高目標値時系列と観測基本周波数系列とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、特性パラメータの更新値と逆フィルタ出力値とを生成する。

Ａ−４−１．数値微分フィルタ係数算出過程
数値微分フィルタ係数算出過程Ａ−４−１では、以上のように、適当な１階数値微分係数｛ａ_ｌ｝と２階数値微分係数｛ｂ_ｌ｝を任意に設定する。

１階数値微分係数｛ａ_ｌ｝は、任意のサンプリング時刻における観測基本周波数の１階微分を近似するために用いる係数である。具体的には、時刻ｎ△における観測基本周波数の１階微分は、結合係数ａ_１−ｎ，…，ａ_Ｎ−ｎによるｙ_１，…，ｙ_Ｎの線形結合によって次式（１６）で近似できる。

この結合係数｛ａ_ｌ｝選び方には任意性があり、例えば、数式（１６）が時刻ｎ△における観測基本周波数の１階微分の良い近似を与える結合係数の例として、次式（１７）が挙げられる（詳細は後述）。

時刻ｎ△における観測基本周波数の１階微分を、前後の観測値ｙ_ｎ＋１，ｙ_ｎ−１の差分（中心差分）で近似する場合が１階数値微分係数のシンプルな例の１つであり、この場合、結合係数は、次式（１８）で表わされる。

他に、後退差分で近似する場合には、次式（１９）となり、

前進差分で近似する場合には、次式（２０）となる。

同様に、各サンプリング時刻ｎ△における観測基本周波数の２階微分は、結合係数ｂ_１−ｎ，…，ｂ_Ｎ−ｎによるｙ_１，…，ｙ_Ｎの線形結合は、次式（２１）によって近似できる。

この結合係数｛ｂ_ｌ｝の選び方には、やはり任意性があり、例えば、数式（２１）が時刻ｎ△における観測基本周波数の２階微分の良い近似を与える結合係数の例として、次式（２２）が挙げられる。

前進差分商と後退差分商の差分商とで近似する場合には、次式（２３）となる。

他にも２階微分の近似を与える結合係数は多数あるが、ここでは省略する。

Ａ−４−２．微分方程式逆フィルタ出力更新値生成過程
この微分方程式逆フィルタ出力更新値生成過程Ａ−４−２では、フレーム毎の特性パラメータ行列Ｗ^（ｆ）及び逆フィルタ出力更新値を生成する。具体的には、前段で出力した音高目標値時系列ｍ_ｑ１，ｍ_ｑ２，…，ｍ_ｑＮと、数値微分フィルタ係数算出過程４−１で設定した１階数値微分係数｛ａ_ｌ｝と２階数値微分係数｛ｂ_ｌ｝とを用いて、次式（２４）を解いて、α^（ｆ）、β^（ｆ）、γ^（ｆ）を求める。

但し、ベクトルｙ＝（ｙ_１，…，ｙ_Ｎ）^Ｔ、ベクトルｍ＝（ｍ_ｑ１，…，ｍ_ｑＮ）^Ｔとし、定数行列Ｂ、Ａは、次式（２５）の通りである。また、Ｂ^ＴはＢの転置行列を表わすものとし、ｙ＝（ｙ_１，…，ｙ_Ｎ）^Ｔは、下記に示すようにｙ_１からｙ_Ｎを縦に並べたベクトルを表わす。

また、行列のｉ，ｊ成分がｘ_ｉ，ｊから構成される行列Ｘ＝（Ｘ_ｉ，ｊ）と書くと、Ｔｒ^（ｆ）（Ｘ）は、次式（２６）で表わされる。但し、Ｄ_ｆは、次式（２７）で表わされる。

求まったα^（ｆ）、β^（ｆ）、γ^（ｆ）と定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ^（ｆ）を、次式（２８）として出力する。これにより、逆フィルタ出力値系列は、次式（２９）により求めて出力する。

なお、逆フィルタ出力値更新過程（Ａ−４）は、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを求める処理である。このときＷを式（２８）で定義し、パラメータα，β，γを求めることで、下記の式（８３）のＪの値を最大にするようなＷを算出する方法が、「微分方程式逆フィルタ出力更新値生成過程（Ａ−４−２）」である。

Ａ−４−３．差分方程式逆フィルタ出力更新値生成過程
該差分方程式逆フィルタ出力更新値生成過程Ａ−４−３では、次式（３０）により計算されるベクトルｗ＝（ｗ_１，…，ｗ_Ｍ−１）^Ｔの要素ｗ_１，…，ｗ_Ｎ−１を、次式（３１）のように代入して特徴パラメータ行列Ｗ^（ｆ）を出力する。

但し、行列ｍ（ｆ）は、次式（３２）で表わされ、行列Ｕは、次式（３３）で表わされる。

これにより、逆フィルタ出力値系列を、次式（３４）で求めて出力する。

但し、行列０_Ｌは、成分がすべて０のＬ次元縦ベクトルとする。
なお、上述した通り、逆フィルタ出力値更新過程（Ａ−４）は、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを求める処理である。このときＷをＴｏｅｐｌｉｔｚ型行列である式（３１）で定義し、下記式（９８）のＪ_１の値を最小にするＷの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を決定することにより、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを算出する方法が、「差分方程式逆フィルタ出力更新値生成過程（Ａ−４−３）」である。つまり、式（３０），式（３１）は、下記式（９８）のＪ_１の値を最小にするＷの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を決定するための具体的な計算方法の一つであり、式（３０），式（３１）を再帰的に計算することによって最適なＷに近づけることができる。

Ａ−５．音高目標値更新過程
基本周波数抽出過程Ａ−１において出力された観測基本周波数時系列ｙ１，…，ｙＮと、逆フィルタ出力値系列を入力として、次式（３５）で表わされる音高目標値時系列を更新する。

前段で求めたベクトル（＝逆フィルタ出力値系列）の成分、すなわち、数式（３６）で示される成分を用いて、音高目標値時系列ｍ_ｑ１，ｍ_ｑ２，…，ｍ_ｑＮを以下の２つのステップにより求める。

Ａ−５−１．音高目標値時系列更新過程（状態系列ｑ_１，…，ｑ_Ｎの決定）
前段で出力されたｍ_１，…，ｍ_Ｉを用い、動的計画法に基づいて、数式（３７）〜（４１）に示すように、以下のような再帰計算を行う。
Ｓｂ１．初期化：

Ｓｂ２．再帰計算：

Ｓｂ３．終了：

Ｓｂ４．バックトラック：

但し、ε_ｋは、観測ｙ_ｋを含んだフレーム番号の集合、すなわち、次式（４２）で表わされる。

Ａ−５−２．音高目標値集合更新過程５−２（音高目標値集合ｍ_１，…，ｍ_Ｉの決定）
上記ステップＳｂ１により求めたｑ_１，…，ｑ_Ｎを用いて、集合Ｃ_ｉ（ｉ＝１，…，Ｉ）を次式（４３）により更新する。

これを用いて、次式（４４）により、ｍ_１，…，ｍ_Ｉを求める。但し、｜Ｃ_ｉ｜、｜ε_ｎ｜は、それぞれ集合Ｃ_ｉ、ε_ｎの要素数を表す。以上より求めたｑ_１，…，ｑ_Ｎとｍ_１，…，ｍ_Ｉを用いて、音高目標値時系列ｍ_ｑ１，ｍ_ｑ２，…，ｍ_ｑＮが求まる。

Ａ−６．収束判定過程
反復計算が所定の回数を満たしたか否か、あるいは、反復計算においてパラメータの更新の変化率が所定値以下になったか否か、あるいは、目的関数値の変化率が所定値以下になったか否かを判定する。

Ｂ．歌唱の基本周波数制御モデル
歌唱の基本周波数（Ｆ０）軌跡は、歌唱者が頭の中で思い描く旋律（ステップ状成分からなる階段関数）に、歌唱者の表現意図や身体的特性による動的変動成分が付加されて観測される。本発明は、観測されるＦ０軌跡から、歌唱者が歌おうとする旋律概形と歌唱者ごとに異なる動的変動成分とを分離する手法に関する。

Ｂ−１．２階の微分方程式に基づくＦ０制御
Ｆ０軌跡の動的特性は、次式（４５）で示される、２次系の制御システムを用いて良く表現できる。

ここで、歌唱者が頭の中で思い描く旋律概形を、ステップ状成分からなる階段関数ｕ（ｔ）と表現する。このｕ（ｔ）を入力としたときに、数式（４５）に従ってＦ０軌跡ｙ（ｔ）が観測されるものとする。はじめに、このシステムの伝達関数Ｇ（ｓ）を求め、その特性を確認する。まず、数式（４５）の両辺をラプラス変換する。Ｙ（ｓ）＝Ｌ［ｙ（ｔ）］とすると、関数ｙ（ｔ）のｎ階導関数ｙ^（ｎ）（ｔ）のラプラス変換は、次式（４６）と書ける。

ｙ（０）＝０、ｕ（０）＝０とし、数式（４５）の両辺をラプラス変換すると、次式（４７）で表される。

伝達関数は、Ｇ（ｓ）＝Ｙ（ｓ）／Ｕ（ｓ）で表されるので、次式（４８）となる。

一方、制動２次系伝達関数Ｈ（ｓ）は、一般的に、次式（４８）のように記述される。

この伝達関数のステップ応答を図１に示す。ζ＝０では、定常振動となり、これは音高安定時間が持続した場合に観測される４〜７［Ｈｚ］の周期的な振動であるビブラートに対応すると考える。さらに、（｜ζ｜＜１）では減衰振動となり、これは音高遷移時に目的音高より大きく振れてしまうオーバーシュートに対応すると考える。このように、音高遷移における基本周波数の挙動を２次系の制御システムは適切に表現できる。また、パラメータΩ、ζ、Ｋと数式（４５）の微分方程式の係数α、β、γとの関係は、次式（５０）の通りである。

Ｂ−２．差分方程式への変換
連続時間信号を扱う数式（４５）の微分方程式を、離散時間信号を扱う差分方程式に変換する。連続時間信号ｙ（ｔ）が、Ｎｙｑｕｉｓｔ周波数より低い帯域制限を受けた信号であると仮定すると、このｙ（ｔ）は、対応する離散時間信号ｙ_ｎと連続時間信号ｓｉｎｃ（πｔ／△）の離散畳込みで、次式（５１）に示すように表現することができる。

ここで、△はサンプリング周期とする。数式（５１）よりｙ（ｔ）の１階微分、２階微分は、それぞれ数式（５２）、（５３）となる。

ここで、ｔ＝ｎ△における１階微分は、次式（５４）で表わされる。

観測Ｎ点だけで上記微分の近似を行うと、次式（５５）となる。

同様に、ｔ＝ｎ△における２階微分は、次式（５６）で表わされる。

これを同様に観測Ｎ点だけで近似を行うと、次式（５７）となる。

数式（５５）と数式（５７）を数式（４５）の微分方程式に代入すると、次式（５８）で表わされる。

上記数式（５８）は、ｙ＝（ｙ_１，ｙ_２，…，ｙ_Ｎ）^Ｔ、ｕ＝（ｕ_１，ｕ_２，…，ｕ_Ｎ）^Ｔと置くと、次式（５９）に示すように簡潔に表現できる。

但し、行列Ｂ、Ａは、次式（６０）の通りである。

数式（５９）の線形方程式のパラメータは、実際には、α、β、γの３つ（行列Ｂ、Ａは定数行列）であるが、次式（６１）のように置き、行列Ｗ＝（Ｗ_ｉ，ｊ）の要素が、ＷがＴｏｅｐｌｉｔｚ型であるという拘束つきの自由パラメータであるような状況も併せて考えることにする（数式（５９）が成立するための必要条件は、ＷがＴｏｅｐｌｉｔｚ型であることによる）。すなわち、α、β、γを推定すべきパラメータとする問題（以降、「ケース１」）、及び、行列Ｗの要素Ｗ_ｉ，ｊをＴｏｅｐｌｉｔｚ行列の拘束条件の下で推定すべきパラメータとする問題（以後、「ケース２」）を、以後同時並行的に検討する。

Ｂ−３．階段関数ｕ_ｎのモデル化
システムの入力となる旋律概形を表す階段関数ｕｎは、図２に示すような状態集合Ｓ＝｛Ｓ_１，…，Ｓ_Ｉ｝からなるＨＭＭ（隠れマルコフモデル）を利用して、以下のようにモデル化する。すなわち、１回の状態遷移によりサンプル点が１個生成されるモデルであり、状態によって出力の統計的な傾向が異なる。ここでは、式の見やすさのため、Ｓ_ｉ＝ｉとし、次式（６２）で表わす。

すなわち、数式（６３）、（６４）で表わされる。

ｍ_Ｓｉは、状態Ｓ_ｉにおける出力確率分布（正規分布）の平均を表す。一様なマルコフ連鎖を想定し、状態Ｓ_ｊから状態Ｓ_ｉへの遷移確率は、^＊Ｐ（Ｓ_ｉ｜Ｓ_ｊ）と表す。ここで、遷移確率^＊Ｐ（Ｓ_ｉ｜Ｓ_ｊ）は、定数とする。前述したＰ_ｉ，ｊは、ｌｏｇ^＊Ｐ（Ｓ_ｉ｜Ｓ_ｊ）のことである。同じ状態へ遷移する遷移確率^＊Ｐ（Ｓ_ｉ｜Ｓ_ｊ）を自己遷移確率といい、これが大きいほど同じ状態に留まろうとする傾向が強くなる。ｑ_ｎは、状態集合Ｓの要素（ＨＭＭの状態番号）の中のいずれかの値をとる。従って、状態系列ｑ_１，…，ｑ_Ｎと、各時刻の状態における出力確率分布の平均ｍ_ｑｎによって階段関数が決定される。

Ｂ−４．Ｆ０制御モデルのパラメータの解釈
上述では、２階の微分方程式を利用したＦ０制御モデルと、入力となる階段階数とをモデル化した。以上より、ケース１では、次式（６５）が、ケース２では、次式（６６）が推定したいパラメータである。

ここで、各パラメータメータの解釈を図３を参照して説明する。まず、微分方程式の係数α、β、γ、あるいは（Ｗ_ｉ，ｊ）は、歌声の動的変動を表すパラメータである。音高が安定するときの振動であったり、音高遷移における連続的なダイナミクスを表現する。状態系列ｑ_１，…，ｑ_Ｎは、音高が安定する長さを決定するパラメータである。これは、必ずしも楽譜に記される音符の長さに対応するわけではなく、歌唱者の意図や、歌唱スタイルに基づいて生成される運動指令の長さを表現したものであると、ここでは想定している。最後に、ＨＭＭの各状態の平均ｍ_Ｓ１，…，ｍ_ＳＩは、歌唱者が意図する旋律の音高（音高目標値）に対応するパラメータである。これは、必ずしも楽譜に記される音符の音高（客観的に定まっている音高値）には対応しない。

Ｃ．Ｆ０制御モデルのパラメータ最尤推定
微分方程式の係数α、β、γと、階段関数をモデル化する状態系列ｑ_１，…，ｑ_Ｎ、各状態の正規分布の平均ｍ_Ｓｌ，…，ｍ_ＳＩを、観測系列ｙ＝（ｙ_１，…，ｙ_Ｎ）^Ｔから最尤推定する方法について述べる。

観測系列が線形差分方程式に理想的に従うならば、数式（５９）を解けば良いが、実際には、理想的な差分方程式からの誤差があると考えられる。そこで、次式（６７）と次式（６８）との間に、数式（６９）が成り立つと仮定する。

但し、行列Ｗ、Ｂ、Ａ、ｍは、それぞれ、次式（７０）、（７１）、（７２）で表わされる。

ここで、εの要素ε_ｎは平均０、分散σ^２の正規分布に従う互いにＧａｕｓｓ性白色雑音である。いま、観測系列ｙに対するパラメータΘの尤度は、数式（６９）より、次式（７３）となるから、多次元正規分布であることが分かり、その正規化係数は、次式（７４）とすればよく、結局、次式（７５）のような形となる。

このとき、パラメータΘの対数尤度は、次式（７６）となる。

パラメータΘの事後確率は、次式（７７）であるので、ここで、次式（７８）が成り立つとすると、次式（７９）と表すことができる。

ここで、事前確率^＊Ｐ（α，β，γ）と^＊Ｐ（ｍ_Ｓ１，…，ｍ_ＳＩ）は一様分布とし、^＊Ｐ（ｑ_１，…，ｑ_Ｎ）は、先に述べたようにＨＭＭによる一様なマルコフ連鎖を想定している。このため、次式（８０）で表わすことができる。

遷移確率（Ｐ）（Ｓｉ｜Ｓｊ）は、事前に決定する定数である。以後、簡単のため、次式（８１）で表わすことにする。したがって、次式（８２）となる。

以上より、数式（７６）と数式（８０）を、数式（７９）に代入し、定数項を除いた次式（８３）がパラメータΘに関して最大化したい目的関数である。

しかし、数式（８３）を最大化するパラメータΘは、解析的に求めることができない。そこで、Θの各要素（微分方程式の係数、ガウス性雑音の分散、状態系列、状態の出力分布の平均）に関して、他の要素を固定した下で、数式（８３）を最大化するステップを、数式（８３）の値が収束するまで繰り返す。

Ｃ−１．特徴パラメータ行列Ｗの更新
状態系列ｑ_１，…，ｑ_Ｎ、ＨＭＭの各状態の正規分布の平均値ｍ_Ｓ１，…，ｍ_ＳＩを固定したとき、数式（８３）を最大にするＷを求めたい。ここでは、これを実現する３つの解法について説明する。

Ｃ−１−１．解法１（微分方程式逆フィルタを用いた方法）
ここでは、ケース１（自由パラメータはα、β、γ）を想定する。上記数式（８３）の右辺の第２項は、第３項に比べてＪへの寄与が無視できるほど小さいと仮定し、次式（８４）を最小化するα、β、γが、Ｊを最大化するものと近似的に見なす。

Ｊ_１をαに関して偏微分すると、次式（８５）を得る。

これを０と置くと、次式（８６）となる。

また、同様に、Ｊ１をβ及びγに関して偏微分して０と置くと、次式（８７）、（８８）を得る。

以上より立てられる次式（８９）の正規方程式を解けばよい。

すなわち、次式（９０）がα、β、γの更新値となる。

また、以上により求めたα、β、γを用い、雑音の分散推定値が次式（９１）により求まる。

Ｃ−１−２．解法２（逆フィルタを用いた方法）
ここでも、ケース１（自由パラメータはα、β、γ）を想定する。解放１のように、数式（８３）の第２項（ｌｏｇ｜Ｗ｜）の寄与を無視しないとすると、Ｊを最大化するα、β、γは解析的に求まらない。そこで、ここでは、α、β、γの更新値を勾配法により数値計算する方法について説明する。以下で与えられる、Ｊのα、β、γに関する勾配ベクトル▽Ｊにより、最急降下法、共役勾配法、準ニュートン法などが適用できる。

また、勾配法により、α、β、γを更新するごとに、雑音の分散推定値σ_２を、次式（９６）により更新する。

Ｃ−１−３．解法３（差分方程式逆フィルタを用いた方法）
ここでは、ケース２（自由パラメータは行列Ｗの要素）を想定する。行列Ｗを、以下のように対角成分が、次式（９７）で示すように、全て１のＴｏｅｐｌｉｔｚ型の上三角行列と仮定し、α、β、γの代わりに、ｗ_１，…，ｗ_Ｍ（Ｍ≦Ｎ−１）をパラメータとして求める。但し、ｗ_ｍ＝０（ｍ≧Ｍ＋１）とする。

すなわち、ここでは、行列Ｗに関して数式（７０）のような要素に関する拘束を仮定しない代わりに、数式（９７）のような構造の拘束を与える。

このとき、｜Ｗ｜＝１であることに注意すると、ｌｏｇ｜Ｗ｜＝０であるので、次式（９８）を最小化するｗ_１，…，ｗ_Ｍを求めればよい。

行列ＷはＴｏｅｐｌｉｔｚ行列のため、ベクトルｗ＝（ｗ_１，…，ｗ_Ｎ−１）^Ｔとすると、次式（９９）で表わされる。

但し、行列Ｕは、次式（１００）で表わされる。

したがって、Ｊ_１は、次式（１０１）で表わされる。

これを行列ｗに関して偏微分して０と置くと、次式（１０２）で表わされる正規方程式を得る。

これを解くと、次式（１０３）を得る。

以上により求まった行列ｗの要素を数式（９７）に基づいて、行列Ｗの中に代入すれば、Ｔｏｅｐｌｉｔｚ型の拘束条件の下での最適な行列Ｗを求めたことになる。また、以上により求まった行列Ｗを用い、雑音の分散推定値が次式（１０４）により求まる。

Ｃ−２．状態系列ｑ_１，…，ｑ_Ｎの推定
特徴パラメータ行列Ｗ、ＨＭＭの各状態における正規分布の平均値ｍ_Ｓ１，…，ｍ_ＳＩを固定したとき、状態系列ｑ_１，…，ｑ_Ｎに関して数式（８３）の最大化を考える。つまり、数式（８３）から関係する項だけを取り出して、次式（１０５）と置き、これを最大にする最適な状態系列ｑ_１，…，ｑ_Ｎを求める。但し、次式（１０６）が成立するものとする。

この問題は、ビタビ（Ｖｉｔｅｒｂｉ）アルゴリズム（動的計画法）により効率的に解くことができる。まず、次のような、最初から時刻ｋに状態Ｓ_ｉに至るまでの部分系列に関する最適な状態系列について、次式（１０７）のような量が定義されているとする。

この量は、漸化式で、次式（１０８）のようにして得られる。

これをｋ＝Ｎまで計算すれば、最適経路が求められる。
図４にアルゴリズムを示す。図において、まず、次式（１０９）に従って初期化を行う（ステップＳｃ１）。

次に、次式（１１０）に従って再帰処理を行う（ステップＳｃ２）。

次に、次式（１１１）に従って終了処理を行う（ステップＳｃ３）。

そして、次式（１１２）で示されるように、状態系列のバックトラックを行う（ステップＳｃ４）。

Ｃ−３．ＨＭＭの各状態の出現確率分布（正規分布）の平均ｍ_Ｓ１，…，ｍ_ＳＩの推定
特徴パラメータ行列Ｗと状態系列ｑ_１，…，ｑ_Ｎを固定したとき、数式（８３）が最大となるように、ＨＭＭの各状態における正規分布の平均ｍ_Ｓ１，…，ｍ_ＳＩを更新する。つまり、数式（８３）からｍ_ｓ１，…，ｍ_ＳＩに関係する項だけを取り出して符号を反転した、次式（１１３）で示される値の最小化を考えればよい。

すなわち、次式（１１４）を解くと、次式（１１５）を得る。

ここで、集合Ｃ_ｉ＝｛ｎ｜ｑ_ｎ＝Ｓ_ｉ｝とし、｜Ｃ_ｉ｜をその要素数とする。数式（１１５）より、各状態における正規分布の平均ｍ_Ｓｉが更新される。以上は、前述したビタビアルゴリズムによって求められた最適な状態系列を用いて、ＨＭＭの各状態のパラメータを学習する操作であるため、しばしばビタビ学習と呼ばれる。最終的に状態系列ｑ_１，…，ｑ_ＮとＨＭＭの各状態の正規分布の平均値ｍ_Ｓ１，…，ｍ_ＳＩから、次式（１１６）で示される入力の階段階数ｕ_ｎ、すなわち、行列ｍが求まる。

Ｃ−４．初期値設定
初期値設定では、Ｂ−１、Ｂ−２、Ｂ−３で説明した３段階のパラメータ推定を、パラメータΘの対数事後確率が収束するまで順番に繰り返す。しかし、以上の反復法によるパラメータ推定は、初期値を適切に設定しないと、局所解に収束してしまう。この問題に対処するため、推定手順の前半に２つの初期値設定、すなわち第１の初期値設定及び第２の初期値設定を行う。

第１の初期値設定では、観測系列ｙに、Ｂ−２で説明したビタビアルゴリズムを適用する。ここでは、次式（１１７）に示す値を最小化する状態系列ｑ１，…，ｑＮを、前述したＢ−２と同様の方法で求め、これを次式（１１８）で示す初期状態系列とする。

第２の初期値設定では、第１の初期値設定で求めた初期状態系列をもとにＨＭＭの各状態の正規分布の平均を求める。すなわち、次式（１１９）、（１２０）を解き、次式（１２１）を得て、これを次式（１２２）とする。

以上のように、初期値設定の段階で、観測系列ｙから、ある程度の階段関数ｕｎの概形を推定することにより、局所解に収束してしまうことを防ぐことが可能となる。

Ｄ．時変なＦ０制御モデルへの拡張
前述では、特徴パラメータ行列Ｗは、時不変なものとして推定を行ったが、観測系列フレームと呼ぶ区間に分割し（但し、区間は重複してもよい）、α、β、γ、あるいは特徴パラメータ行列Ｗを、フレーム毎に自由度もつパラメータと見なして推定するようにしてもよい。

（実施例）
図５は、本実施形態による、基本周波数軌跡モデルパラメータ抽出装置の構成を示すブロック図である。図において、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列を抽出する。初期音高目標値生成部２は、抽出された観測基本周波数時系列を入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列を生成する。フレーム分割部３は、観測基本周波数時系列をフレーム分割する。

逆フィルタ出力値系列更新部４は、初期音高目標値時系列とフレーム分割された観測基本周波数系列とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、観測基本周波数系列から逆フィルタ出力値系列と特性パラメータとを生成する。音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とを入力として、音高目標値時系列を生成する。

収束判定部６は、特性パラメータと音高目標値時系列とが所定の規準を満たしているかどうかを判定し、満たしていない場合には、逆フィルタ出力値系列更新部４と音高目標値更新部５とに再処理させる。基本周波数軌跡パラメータ出力部７は、収束判定部６において所定の規準を満たすと判定された特性パラメータと音高目標値時系列とを出力する。

図６は、本実施形態において、逆フィルタ出力値系列更新部４で、前述した逆フィルタを用いる第１動作例（Ｃ−１−２．解法２）を説明するためのフローチャートである。まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＡ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍ^ｉｎｉｔを生成する（ステップＳＡ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）を出力する（ステップＳＡ３）。

次に、逆フィルタ出力値系列更新部４は、初期音高目標値時系列ｍ^ｉｎｉｔとフレーム分割された観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、観測基本周波数系列から逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と特性パラメータＷ^（１），…，Ｗ^（Ｆ）とを生成する（ステップＳＡ４）。次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と観測基本周波数系列ｙ^（１），…，ｙ^（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＡ４）。

収束判定部６は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＡ６）、満たしていない場合には（ステップＳＡ６のＮＧ）、ステップＳＡ４に戻り、逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＡ６の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＡ７）。

次に、図７は、本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−１、Ａ−４−２で説明した微分方程式逆フィルタを用いる第１動作例（Ｃ−１−１．解法１）を説明するためのフローチャートである。
まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＢ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍ^ｉｎｉｔを生成する（ステップＳＢ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）を出力する（ステップＳＢ３）。

次に、逆フィルタ出力値系列更新部４は、前述したＡ−４−１で説明したように、適当な１階数値微分係数と２階数値微分係数を用いて、数値微分フィルタ係数Ａ、Ｂを算出し（ステップＳＢ４）、該数値微分フィルタ係数Ａ、Ｂと、初期音高目標値時系列ｍ^ｉｎｉｔと、フレーム分割された観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）とを入力として、前述したＡ−４−２で説明したように、適当な１階数値微分係数と２階数値微分係数を用いて、フレーム毎の逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と特性パラメータＷ^（１），…，Ｗ^（Ｆ）とを生成する（ステップＳＢ４）。

次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と観測基本周波数系列ｙ^（１），…，ｙ^（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＢ５）。収束判定部６は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＢ７）、満たしていない場合には（ステップＳＢ７のＮＧ）、ステップＳＢ５に戻り、微分方程式逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＢ７の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＢ８）。

次に、図８は、本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−３で説明した差分方程式逆フィルタを用いる第２動作例（Ｃ−１−３．解法３）を説明するためのフローチャートである。まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＣ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍ^ｉｎｉｔを生成する（ステップＳＣ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）を出力する（ステップＳＣ３）。

次に、逆フィルタ出力値系列更新部４は、初期音高目標値時系列ｍ^ｉｎｉｔと、フレーム分割された観測基本周波数系列У^（１），…，У^（Ｆ）、ｙ^（１），…，ｙ^（Ｆ）とを入力として、前述したＡ−４−３で説明したように、フレーム毎の逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と特性パラメータＷ^（１），…，Ｗ^（Ｆ）とを生成する（ステップＳＣ４）。

次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾^（１），…，ｕ＾^（Ｆ）と観測基本周波数系列ｙ^（１），…，ｙ^（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＣ５）。収束判定部６は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＣ６）、満たしていない場合には（ステップＳＣ６のＮＧ）、ステップＳＣ４に戻り、差分方程式逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＣ６の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ^（１），…，Ｗ^（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＣ７）。

上述した実施形態によれば、歌声知覚に対しての妥当性が確認された２次系歌声生成モデルの制御パラメータを実データから教師なし学習する枠組みが確立されれば、歌声合成技術が飛躍的に向上させることができる。例えば、従来不可能であった歌唱スタイルの転写が実現される。歌唱者Ａの歌い方を制御パラメータから学習し、旋律ａを表す階段状の信号に畳み込むことによって、歌唱者Ａの歌い方による旋律ａの歌声を自動生成できる。

また、上述した実施形態によれば、表現豊かな歌い方を制御パラメータによって予め学習することで、高品質な歌声合成が期待される。

さらに、カラオケ採点システムへの応用も考えられる。これまで、歌唱力の自動評価手法が提案されているが、聴取実験に基づいて様々な評価尺度を提案し、歌の上手・下手の２クラス識別を行っているものの、歌い方を精密にモデル化するまでには至っていない。これに対して、上述した実施形態によれば、制御パラメータを評価尺度に利用することで、性能向上を図ることができる。

また、上述した実施形態においては、基本周波数抽出部１、初期音高目標値生成部２、フレーム分割部３、逆フィルタ出力値系列更新部４、音高目標値更新部５、収束判定部６、基本周波数軌跡パラメータ出力部７は、プログラムの形式でコンピュータ読み取り可能な記録媒体に記憶されていてもよく、このプログラムをコンピュータが読み出して実行することによって、上記処理を行うようにしてもよい。すなわち、放送装置５における、各処理手段、処理部は、ＣＰＵ等の中央演算処理装置がＲＯＭやＲＡＭ等の主記憶装置に上記プログラムを読み出して、情報の加工・演算処理を実行することにより、実現するようにしてもよい。

ここでコンピュータ読み取り可能な記録媒体とは、磁気ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、半導体メモリ等をいう。また、このコンピュータプログラムを通信回線によってコンピュータに配信し、この配信を受けたコンピュータが当該プログラムを実行するようにしても良い。

２次系伝達関数のステップ応答を示す概念図である。 階段関数ｕ_ｎのモデル化を説明するための概念図である。 Ｆ０制御パラメータの解釈を説明するための概念図である。 最適な状態系列を求めるためのビタビアルゴリズムの一例を示す概念図である。 本実施形態による、基本周波数軌跡モデルパラメータ抽出装置の構成を示すブロック図である。 本実施形態において、逆フィルタ出力値系列更新部４で、前述した逆フィルタを用いる第１動作例（Ｃ−１−２．解法２）を説明するためのフローチャートである。 本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−１、Ａ−４−２で説明した微分方程式逆フィルタを用いる第１動作例（Ｃ−１−１．解法１）を説明するためのフローチャートである。 本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−３で説明した差分方程式逆フィルタを用いる第２動作例（Ｃ−１−３．解法３）を説明するためのフローチャートである。

符号の説明

１ 基本周波数抽出部
２ 初期音高目標値生成部
３ フレーム分割部
４ 逆フィルタ出力値系列更新部
５ 音高目標値更新部
６ 収束判定部
７ 基本周波数軌跡パラメータ出力部

Claims (8)

1. 基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置であって、
   入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出手段と、
   前記基本周波数抽出手段により抽出された観測基本周波数時系列をフレーム分割するフレーム分割手段と、
   前記基本周波数抽出手段により抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成手段と、
   前記初期音高目標値生成手段により生成された初期音高目標値時系列と前記フレーム分割手段によりフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新手段と、
   全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新手段と、
   前記逆フィルタ出力値系列更新手段により生成された特性パラメータと前記音高目標値更新手段により生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新手段と前記音高目標値更新手段による処理を再度行なわせる収束判定手段と、
   前記収束判定手段により所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力手段と
   を備え、
   前記逆フィルタ出力値系列更新手段は、
   前記観測基本周波数時系列ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ 、前記音高目標値時系列ｍ＝（ｍ _ｑ１ ，…，ｍ _ｑＮ ） ^Ｔ とし、
   定数行列Ｂ、Ａを、
   

   

   とし、
   Ｂ ^Ｔ はＢの転置行列を表わすものとし、
   ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ は、ｙ _１ からｙ _Ｎ を縦に並べたベクトルを表わし、
   行列のｉ，ｊ成分がｘ _ｉ，ｊ から構成される行列Ｘ＝（Ｘ _ｉ，ｊ ）と書くとき、Ｔｒ ^（ｆ） （Ｘ）を、
   

   

   で表わし、Ｄ _ｆ を
   

   

   で表した場合に
   

   

   を解いて、α ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） を求め、
   求まったα ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） と前記定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ ^（ｆ） を、
   

   

   により算出して出力し、これにより、前記逆フィルタ出力値系列を、
   

   

   により求めて出力し、
   前記音高目標値更新手段は、
   前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、
   前記逆フィルタ出力値系列の成分である
   

   

   と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である
   

   

   を用いて、
   

   

   

   

   

   

   の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新する
   ことを特徴とする基本周波数軌跡モデルパラメータ抽出装置。
2. 基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置であって、
   入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出手段と、
   前記基本周波数抽出手段により抽出された観測基本周波数時系列をフレーム分割するフレーム分割手段と、
   前記基本周波数抽出手段により抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成手段と、
   前記初期音高目標値生成手段により生成された初期音高目標値時系列と前記フレーム分割手段によりフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新手段と、
   全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新手段と、
   前記逆フィルタ出力値系列更新手段により生成された特性パラメータと前記音高目標値更新手段により生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新手段と前記音高目標値更新手段による処理を再度行なわせる収束判定手段と、
   前記収束判定手段により所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力手段と
   を備え、
   前記逆フィルタ出力値系列更新手段は、
   行列ｍ ^（ｆ） を
   

   

   で表わし、行列Ｕを
   

   

   で表した場合に、
   

   

   により計算されるベクトルｗ＝（ｗ _１ ，…，ｗ _Ｍ−１ ） ^Ｔ の要素ｗ _１ ，…，ｗ _Ｎ−１ を、
   

   

   に代入して特徴パラメータ行列Ｗ ^（ｆ） を出力し、
   

   

   により前記逆フィルタ出力値系列を求めて出力し、
   前記音高目標値更新手段は、
   前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、
   前記逆フィルタ出力値系列の成分である
   

   

   と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である
   

   

   を用いて、
   

   

   

   

   

   

   の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新する
   ことを特徴とする基本周波数軌跡モデルパラメータ抽出装置。
3. 前記逆フィルタ出力値系列更新手段は、
   Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、
   

   

   のＪ_１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を求めることにより前記特徴パラメータ行列Ｗを決定する
   ことを特徴とする請求項２に記載の基本周波数軌跡モデルパラメータ抽出装置。
4. 基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出方法であって、
   入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、
   前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、
   前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、
   前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、
   全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、
   前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、
   前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、
   前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップとを含み、
   前記逆フィルタ出力値系列更新ステップにおいて、
   前記観測基本周波数時系列ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ 、前記音高目標値時系列ｍ＝（ｍ _ｑ１ ，…，ｍ _ｑＮ ） ^Ｔ とし、
   定数行列Ｂ、Ａを、
   

   

   とし、
   Ｂ ^Ｔ はＢの転置行列を表わすものとし、
   ｙ＝（ｙ _１ ，…，ｙ _Ｎ ） ^Ｔ は、ｙ _１ からｙ _Ｎ を縦に並べたベクトルを表わし、
   行列のｉ，ｊ成分がｘ _ｉ，ｊ から構成される行列Ｘ＝（Ｘ _ｉ，ｊ ）と書くとき、Ｔｒ ^（ｆ） （Ｘ）を、
   

   

   で表わし、Ｄ _ｆ を
   

   

   で表した場合に
   

   

   を解いて、α ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） を求め、
   求まったα ^（ｆ） 、β ^（ｆ） 、γ ^（ｆ） と前記定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ ^（ｆ） を、
   

   

   により算出して出力し、これにより、前記逆フィルタ出力値系列を、
   

   

   により求めて出力し、
   前記音高目標値更新ステップにおいて、
   前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、
   前記逆フィルタ出力値系列の成分である
   

   

   と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である
   

   

   を用いて、
   

   

   

   

   

   

   の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新する
   ことを特徴とする基本周波数軌跡モデルパラメータ抽出方法。
5. 基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出方法であって、
   入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、
   前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、
   前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、
   前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、
   全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、
   前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、
   前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、
   前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップとを含み、
   前記逆フィルタ出力値系列更新ステップにおいて、
   行列ｍ ^（ｆ） を
   

   

   で表わし、行列Ｕを
   

   

   で表した場合に、
   

   

   により計算されるベクトルｗ＝（ｗ _１ ，…，ｗ _Ｍ−１ ） ^Ｔ の要素ｗ _１ ，…，ｗ _Ｎ−１ を、
   

   

   に代入して特徴パラメータ行列Ｗ ^（ｆ） を出力し、
   

   

   により前記逆フィルタ出力値系列を求めて出力し、
   前記音高目標値更新ステップにおいて、
   前記基本周波数抽出において出力された観測基本周波数時系列ｙ _１ ，…，ｙ _Ｎ と、前記逆フィルタ出力値系列を入力して、
   前記逆フィルタ出力値系列の成分である
   

   

   と、前記観測基本周波数時系列の要素ｙを含んだフレーム番号の集合である
   

   

   を用いて、
   

   

   

   

   

   

   の各式による動的計画法に基づいて再帰計算を行い、前記音高目標値時系列ｍ _ｑ１ ，ｍ _ｑ２ ，…，ｍ _ｑＮ を更新する
   ことを特徴とする基本周波数軌跡モデルパラメータ抽出方法。
6. 前記逆フィルタ出力値系列更新ステップは、
   Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、
   

   

   のＪ_１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ_１，ｗ_２，・・・，ｗ_Ｎ−１を求めることにより前記特徴パラメータ行列Ｗを決定する
   ことを特徴とする請求項５記載の基本周波数軌跡モデルパラメータ抽出方法。
7. 請求項４〜６の何れか一項に記載の基本周波数軌跡モデルパラメータ抽出方法の各ステップをコンピュータにより実行させるためのプログラム。
8. 請求項７に記載のプログラムを記録したコンピュータ読み取り可能な記録媒体。

Images