(57)【要約】【課題】動的変動そのものを適切にモデル化し、Ｆ０軌跡のみから旋律を構成する音高列を正しく推定する。【解決手段】初期音高目標値生成部２は、基本周波数抽出部１が抽出した観測基本周波数時系列から初期音高目標値系列を生成する。逆フィルタ出力値系列更新部４は、初期音高目標値時系列とフレーム分割部３によりフレーム分割された観測基本周波数系列とから逆フィルタ出力値系列と特性パラメータとを生成する。音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから音高目標値時系列を生成する。収束判定部６は、特性パラメータと音高目標値時系列とが所定の規準を満たしているかどうかを判定し、所定の規準を満たすと判定された場合、基本周波数軌跡パラメータ出力部７は、特性パラメータと音高目標値時系列とを出力する。【選択図】図５
2010020258.tif 000002
【特許請求の範囲】
【請求項１】
基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置であって、
入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出手段と、
前記基本周波数抽出手段により抽出された観測基本周波数時系列をフレーム分割するフレーム分割手段と、
前記基本周波数抽出手段により抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成手段と、
前記初期音高目標値生成手段により生成された初期音高目標値時系列と前記フレーム分割手段によりフレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新手段と、
全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新手段と、
前記逆フィルタ出力値系列更新手段により生成された特性パラメータと前記音高目標値更新手段により生成された音高目標値時系列とが所定の規準を満たしているか否かを判定し、所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新手段と前記音高目標値更新手段による処理を再度行なわせる収束判定手段と、
前記収束判定手段により所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力手段と
を備えることを特徴とする基本周波数軌跡モデルパラメータ抽出装置。
【請求項２】
前記初期音高目標値生成手段は、
各時刻において、
音高目標値集合の中から適切なものを選択し、音高目標値時系列を生成する音高目標値時系列生成手段と、
前記音高目標値時系列が各時刻で取り得る有限個の音高目標値の集合である音高目標値集合を生成する音高目標値集合生成手段と
を備えることを特徴とする請求項１記載の基本周波数軌跡モデルパラメータ抽出装置。
【請求項３】
前記逆フィルタ出力値系列更新手段は、
１階数値微分フィルタおよび２階数値微分フィルタの係数を算出する数値微分フィルタ係数算出手段と、
前記数値微分フィルタ係数算出手段により算出された係数に基づいて立てられる特性パラメータに関する正規方程式に基づいて、特性パラメータの更新値と逆フィルタ出力値の更新値とを生成する逆フィルタ出力更新値生成手段と
を更に備えることを特徴とする請求項１記載の基本周波数軌跡モデルパラメータ抽出装置。
【請求項４】
前記逆フィルタ出力値系列更新手段は、
Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、
【数１】
2010020258.tif 000125
のＪ１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ１，ｗ２，・・・，ｗＮ−１を求めることにより前記特徴パラメータ行列Ｗを決定し、
前記特徴パラメータ行列Ｗと前記観測基本周波数時系列とから逆フィルタ出力値系列を算出する
ことを特徴とする請求項１に記載の基本周波数軌跡モデルパラメータ抽出装置。
【請求項５】
基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出方法であって、
入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、
前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、
前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、
前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、
全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、
前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、
前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、
前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップと
を含むことを特徴とする基本周波数軌跡モデルパラメータ抽出方法。
【請求項６】
前記初期音高目標値生成ステップは、
各時刻において、
音高目標値集合の中から適切なものを選択し、音高目標値時系列を生成する音高目標値時系列生成ステップと、
前記音高目標値時系列が各時刻で取り得る有限個の音高目標値の集合である音高目標値集合を生成する音高目標値集合生成ステップと
を含むことを特徴とする請求項５記載の基本周波数軌跡モデルパラメータ抽出方法。
【請求項７】
前記逆フィルタ出力値系列更新ステップは、
１階数値微分フィルタおよび２階数値微分フィルタの係数を算出する数値微分フィルタ係数算出ステップと、
前記算出された係数に基づいて立てられる特性パラメータに関する正規方程式に基づいて、特性パラメータ更新値と逆フィルタ出力値系列とを生成する逆フィルタ出力更新値生成ステップと
を更に含むことを特徴とする請求項５記載の基本周波数軌跡モデルパラメータ抽出方法。
【請求項８】
前記逆フィルタ出力値系列更新ステップは、
Ｔｏｅｐｌｉｔｚ型行列で示される特徴パラメータ行列Ｗにおいて、
【数２】
2010020258.tif 000126
のＪ１の値が最小となるような前記特徴パラメータ行列Ｗの各成分ｗ１，ｗ２，・・・，ｗＮ−１を求めることにより前記特徴パラメータ行列Ｗを決定し、
前記特徴パラメータ行列Ｗと前記観測基本周波数時系列とから逆フィルタ出力値系列を算出する
ことを特徴とする請求項５記載の基本周波数軌跡モデルパラメータ抽出方法。
【請求項９】
基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置のコンピュータに、
入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、
前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、
前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、
前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、
全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、
前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、
前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、
前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップと
を実行させることを特徴とするプログラム。
【請求項１０】
基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定する基本周波数軌跡モデルパラメータ抽出装置のコンピュータに、
入力される音響信号から観測基本周波数時系列を抽出する基本周波数抽出ステップと、
前記抽出された観測基本周波数時系列をフレーム分割するフレーム分割ステップと、
前記抽出された観測基本周波数時系列に基づいて、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する初期音高目標値系列を生成する初期音高目標値生成ステップと、
前記生成された初期音高目標値時系列と前記フレーム分割された観測基本周波数系列とに基づいて、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表す、フィルタの特性パラメータを更新し、逆フィルタ出力値系列と特性パラメータとを生成する逆フィルタ出力値系列更新ステップと、
全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とから、音高目標値時系列を生成する音高目標値更新ステップと、
前記生成された特性パラメータと前記生成された音高目標値時系列とが所定の規準を満たしているか否かを判定する判定ステップと、
前記所定の規準を満たしていないと判定された場合に、所定の規準を満たすまで、前記逆フィルタ出力値系列更新ステップと前記音高目標値更新ステップによる処理を再度行なわせる再帰ステップと、
前記所定の規準を満たすと判定された場合に、その特性パラメータと音高目標値時系列とを出力する基本周波数軌跡パラメータ出力ステップと
を実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、基本周波数軌跡から、音高目標値および基本周波数生成系のフィルタ特性を推定するための基本周波数軌跡モデルパラメータ抽出装置、基本周波数軌跡モデルパラメータ抽出方法、プログラム及び記録媒体に関する。
【背景技術】
【０００２】
歌唱音響信号から推定される基本周波数（Ｆ０）軌跡には、歌唱者が歌おうとする音高目標値の時系列と歌唱力・歌唱スタイル・個人性・感情に基づく様々な動的変動が観測される。歌声は、多くのジャンルの音楽を特徴付ける重要な要素の１つであり、この歌声のＦ０軌跡に着目した様々な研究が現在盛んに行われている。特に、歌声から楽曲を検索するハミング検索では、歌唱された歌声のＦ０軌跡から、歌唱者の意図する音高列を正しく推定して、楽曲データベースの旋律と照合する必要がある。
【０００３】
従来、Ｆ０軌跡そのものをＤＰマッチングによって照合する方法が提案されている（例えば、非特許文献１、非特許文献２参照）。しかしながら、これら従来技術では、歌声の動的変動の影響を受けて検索性能が低下するという問題があった。
【０００４】
ところで、オーバーシュートや、ビブラートのようなＦ０軌跡の動的変動は、歌声知覚に影響を与え、歌声の自然性を保つためには必要不可欠な成分であることが知られている。そこで、これらの動的変動を制御する２次系Ｆ０制御モデルが提案され、自然性、かつ、明瞭性のある歌声合成技術が実現されている（例えば、非特許文献３参照）。つまり、楽譜に相当する階段状の信号に２次系のインパルス応答を畳み込むことが、歌声のＦ０軌跡の生成モデルとして有効であることが示された。
【非特許文献１】橋口博樹、西村拓一、張建新、滝田順子、岡隆一、“モデル依存傾斜制限型の連続ＤＰを用いた鼻歌入力による楽曲信号のスポッティング検索,”電子情報通信学会論文誌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で表わすこととする。
【００２０】
Ａ−１．基本周波数抽出過程
音響信号を入力として観測基本周波数時系列を出力する。出力した観測基本周波数時系列を次式（１）と表す。
【００２１】
【数１】
2010020258.tif 000003
【００２２】
以後、ｙｎのことをサンプルｎにおける観測基本周波数と呼ぶことにする。また、サンプルｎは、時刻ｎ△に対応するものとする。ここで、△をサンプリング周期と呼ぶ。
【００２３】
Ａ−２．フレーム分割過程
基本周波数抽出過程１において出力された観測基本周波数時系列を適当なフレーム長およびフレームシフト長により複数のフレームに分ける。フレーム番号をｆとし、ｆ番目のフレームにおける観測基本周波数の集合を、次式（２）と置く。
【００２４】
【数２】
2010020258.tif 000004
【００２５】
また、集合У（ｆ）の要素を縦に並べたベクトルを、次式（３）と置く。
【００２６】
【数３】
2010020258.tif 000005
【００２７】
但し、Ｔはフレームシフト長、Ｍはフレーム長である。また、Fはフレーム数である。例えば、フレーム長Ｍが４点、フレームシフト長Ｔが２点の場合、次式（４）に示すようになる。
【００２８】
【数４】
2010020258.tif 000006
【００２９】
Ａ−３．初期音高目標値生成過程
基本周波数抽出過程１において出力された観測基本周波数時系列ｙ１，…，ｙＮを入力として、次式（５）で示される、音高目標値時系列の初期値を生成する。
【００３０】
【数５】
2010020258.tif 000007
【００３１】
但し、ｑｎを「状態」と呼び、各時刻ｎおいて１，２，…，Ｉ（Ｉは正の整数）の中のいずれかの値をとる。例えば、Ｉ＝３とすると、次式（６）のようになる。
【００３２】
【数６】
2010020258.tif 000008
【００３３】
ｍｉは、ｉに対応した実数値を表し、ｍ１，…，ｍｉ，…，ｍＩを「音高目標値集合」と呼ぶ。従って、ｍｑ１，…，ｍｑｎは、状態系列ｑ１，…，ｑＮに対応した実数値列を表し、これを「音高目標値時系列」と呼ぶ。例えば、ｍ１＝５０、ｍ２＝１５０、ｍ３＝１２５とし、ｑｎを数式（６）とすると、ｍｑ１，…，ｍｑｎは、次式（７）となる。
【００３４】
【数７】
2010020258.tif 000009
【００３５】
音高目標値時系列の初期値ｍｑ１，ｍｑ２，…，ｍｑＮは、具体的には以下のような２つのステップにより求める。
【００３６】
Ａ−３−１．音高目標値時系列生成過程（状態系列ｑ１，…，ｑＮの決定）
ｍ１，…，ｍＩを適当に設定し（例えば、１２平均律音階に対応する周波数値）、動的計画法に基づいて、数式（８）〜（１１）に示すように、以下のような再帰計算を行なう。
Ｓａ１．初期化：
【００３７】
【数８】
2010020258.tif 000010
【００３８】
Ｓａ２．再帰計算：
【００３９】
【数９】
2010020258.tif 000011
【００４０】
Ｓａ３．終了：
【００４１】
【数１０】
2010020258.tif 000012
【００４２】
Ｓａ４．バックトラック：
【００４３】
【数１１】
2010020258.tif 000013
【００４４】
Ａ−３−２．音高目標値集合生成過程（音高目標値集合ｍｌ，…，ｍＩの決定）
ステップＳａ１により決まったｑ１，…，ｑＮをそれぞれ次式（１２）と置き、次式（１３）により、数式（１４）を求める。
【００４５】
【数１２】
2010020258.tif 000014
【００４６】
【数１３】
2010020258.tif 000015
【００４７】
【数１４】
2010020258.tif 000016
【００４８】
以上より求まった、数式（１２）、（１４）を用いて、次式（１５）で示す、音高目標値時系列の初期値が求まる。
【００４９】
【数１５】
2010020258.tif 000017
【００５０】
但し、Pｉ，ｊ（ｉ，ｊ∈｛１，…，Ｉ｝）は、予め設定しておく定数であり、「状態ｉから状態ｊへの行きやすさ（難しさ）」を表す。例えば、Ｐｉ，ｊを大きめ、Ｐｉ，ｊ（ｉ≠ｊ）を小さめにとると、ｍｑ１，…，ｍｑｎは、移り変わりの少ない安定した階段状の系列として推定されやすくなる。逆に、Ｐｉ，ｊを小さめ、Ｐｉ，ｊ（ｉ≠ｊ）を大きめにとると、ｍｑ１，…，ｍｑｎは、同じ値に長く留まろうとしないようになり、移り変わりの激しい系列として推定されやすくなる。このように、定数Ｐｉ，ｊ（ｉ，ｊ∈｛１，…，Ｉ｝）は、音高目標値系列の移り変わりの激しさを調節するための定数である。
【００５１】
Ａ−４．逆フィルタ出力値更新過程
音高目標値時系列と観測基本周波数系列とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、特性パラメータの更新値と逆フィルタ出力値とを生成する。
【００５２】
Ａ−４−１．数値微分フィルタ係数算出過程
数値微分フィルタ係数算出過程Ａ−４−１では、以上のように、適当な１階数値微分係数｛ａｌ｝と２階数値微分係数｛ｂｌ｝を任意に設定する。
【００５３】
１階数値微分係数｛ａｌ｝は、任意のサンプリング時刻における観測基本周波数の１階微分を近似するために用いる係数である。具体的には、時刻ｎ△における観測基本周波数の１階微分は、結合係数ａ１−ｎ，…，ａＮ−ｎによるｙ１，…，ｙＮの線形結合によって次式（１６）で近似できる。
【００５４】
【数１６】
2010020258.tif 000018
【００５５】
この結合係数｛ａｌ｝選び方には任意性があり、例えば、数式（１６）が時刻ｎ△における観測基本周波数の１階微分の良い近似を与える結合係数の例として、次式（１７）が挙げられる（詳細は後述）。
【００５６】
【数１７】
2010020258.tif 000019
【００５７】
時刻ｎ△における観測基本周波数の１階微分を、前後の観測値ｙｎ＋１，ｙｎ−１の差分（中心差分）で近似する場合が１階数値微分係数のシンプルな例の１つであり、この場合、結合係数は、次式（１８）で表わされる。
【００５８】
【数１８】
2010020258.tif 000020
【００５９】
他に、後退差分で近似する場合には、次式（１９）となり、
【００６０】
【数１９】
2010020258.tif 000021
【００６１】
前進差分で近似する場合には、次式（２０）となる。
【００６２】
【数２０】
2010020258.tif 000022
【００６３】
同様に、各サンプリング時刻ｎ△における観測基本周波数の２階微分は、結合係数ｂ１−ｎ，…，ｂＮ−ｎによるｙ１，…，ｙＮの線形結合は、次式（２１）によって近似できる。
【００６４】
【数２１】
2010020258.tif 000023
【００６５】
この結合係数｛ｂｌ｝の選び方には、やはり任意性があり、例えば、数式（２１）が時刻ｎ△における観測基本周波数の２階微分の良い近似を与える結合係数の例として、次式（２２）が挙げられる。
【００６６】
【数２２】
2010020258.tif 000024
【００６７】
前進差分商と後退差分商の差分商とで近似する場合には、次式（２３）となる。
【００６８】
【数２３】
2010020258.tif 000025
【００６９】
他にも２階微分の近似を与える結合係数は多数あるが、ここでは省略する。
【００７０】
Ａ−４−２．微分方程式逆フィルタ出力更新値生成過程
この微分方程式逆フィルタ出力更新値生成過程Ａ−４−２では、フレーム毎の特性パラメータ行列Ｗ（ｆ）及び逆フィルタ出力更新値を生成する。具体的には、前段で出力した音高目標値時系列ｍｑ１，ｍｑ２，…，ｍｑＮと、数値微分フィルタ係数算出過程４−１で設定した１階数値微分係数｛ａｌ｝と２階数値微分係数｛ｂｌ｝とを用いて、次式（２４）を解いて、α（ｆ）、β（ｆ）、γ（ｆ）を求める。
【００７１】
【数２４】
2010020258.tif 000026
但し、ベクトルｙ＝（ｙ１，…，ｙＮ）Ｔ、ベクトルｍ＝（ｍｑ１，…，ｍｑＮ）Ｔとし、定数行列Ｂ、Ａは、次式（２５）の通りである。また、ＢＴはＢの転置行列を表わすものとし、ｙ＝（ｙ１，…，ｙＮ）Ｔは、下記に示すようにｙ１からｙＮを縦に並べたベクトルを表わす。
【００７２】
【数２５】
2010020258.tif 000027
【００７３】
また、行列のｉ，ｊ成分がｘｉ，ｊから構成される行列Ｘ＝（Ｘｉ，ｊ）と書くと、Ｔｒ（ｆ）（Ｘ）は、次式（２６）で表わされる。但し、Ｄｆは、次式（２７）で表わされる。
【００７４】
【数２６】
2010020258.tif 000028
【００７５】
【数２７】
2010020258.tif 000029
【００７６】
求まったα（ｆ）、β（ｆ）、γ（ｆ）と定数行列Ｂ、Ａとを用い、特徴パラメータ行列Ｗ（ｆ）を、次式（２８）として出力する。これにより、逆フィルタ出力値系列は、次式（２９）により求めて出力する。
【００７７】
【数２８】
2010020258.tif 000030
【００７８】
【数２９】
2010020258.tif 000031
なお、逆フィルタ出力値更新過程（Ａ−４）は、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを求める処理である。このときＷを式（２８）で定義し、パラメータα，β，γを求めることで、下記の式（８３）のＪの値を最大にするようなＷを算出する方法が、「微分方程式逆フィルタ出力更新値生成過程（Ａ−４−２）」である。
【００７９】
Ａ−４−３．差分方程式逆フィルタ出力更新値生成過程
該差分方程式逆フィルタ出力更新値生成過程Ａ−４−３では、次式（３０）により計算されるベクトルｗ＝（ｗ１，…，ｗＭ−１）Ｔの要素ｗ１，…，ｗＮ−１を、次式（３１）のように代入して特徴パラメータ行列Ｗ（ｆ）を出力する。
【００８０】
【数３０】
2010020258.tif 000032
【００８１】
【数３１】
2010020258.tif 000033
【００８２】
但し、行列ｍ（ｆ）は、次式（３２）で表わされ、行列Ｕは、次式（３３）で表わされる。
【００８３】
【数３２】
2010020258.tif 000034
【００８４】
【数３３】
2010020258.tif 000035
【００８５】
これにより、逆フィルタ出力値系列を、次式（３４）で求めて出力する。
【００８６】
【数３４】
2010020258.tif 000036
【００８７】
但し、行列０Ｌは、成分がすべて０のＬ次元縦ベクトルとする。
なお、上述した通り、逆フィルタ出力値更新過程（Ａ−４）は、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを求める処理である。このときＷをＴｏｅｐｌｉｔｚ型行列である式（３１）で定義し、下記式（９８）のＪ１の値を最小にするＷの各成分ｗ１，ｗ２，・・・，ｗＮ−１を決定することにより、下記の式（８３）のＪの値を最大にするような特徴パラメータ行列Ｗを算出する方法が、「差分方程式逆フィルタ出力更新値生成過程（Ａ−４−３）」である。つまり、式（３０），式（３１）は、下記式（９８）のＪ１の値を最小にするＷの各成分ｗ１，ｗ２，・・・，ｗＮ−１を決定するための具体的な計算方法の一つであり、式（３０），式（３１）を再帰的に計算することによって最適なＷに近づけることができる。
【００８８】
Ａ−５．音高目標値更新過程
基本周波数抽出過程Ａ−１において出力された観測基本周波数時系列ｙ１，…，ｙＮと、逆フィルタ出力値系列を入力として、次式（３５）で表わされる音高目標値時系列を更新する。
【００８９】
【数３５】
2010020258.tif 000037
【００９０】
前段で求めたベクトル（＝逆フィルタ出力値系列）の成分、すなわち、数式（３６）で示される成分を用いて、音高目標値時系列ｍｑ１，ｍｑ２，…，ｍｑＮを以下の２つのステップにより求める。
【００９１】
【数３６】
2010020258.tif 000038
【００９２】
Ａ−５−１．音高目標値時系列更新過程（状態系列ｑ１，…，ｑＮの決定）
前段で出力されたｍ１，…，ｍＩを用い、動的計画法に基づいて、数式（３７）〜（４１）に示すように、以下のような再帰計算を行う。
Ｓｂ１．初期化：
【００９３】
【数３７】
2010020258.tif 000039
【００９４】
Ｓｂ２．再帰計算：
【００９５】
【数３８】
2010020258.tif 000040
【００９６】
【数３９】
2010020258.tif 000041
【００９７】
Ｓｂ３．終了：
【００９８】
【数４０】
2010020258.tif 000042
【００９９】
Ｓｂ４．バックトラック：
【０１００】
【数４１】
2010020258.tif 000043
【０１０１】
但し、εｋは、観測ｙｋを含んだフレーム番号の集合、すなわち、次式（４２）で表わされる。
【０１０２】
【数４２】
2010020258.tif 000044
【０１０３】
Ａ−５−２．音高目標値集合更新過程５−２（音高目標値集合ｍ１，…，ｍＩの決定）
上記ステップＳｂ１により求めたｑ１，…，ｑＮを用いて、集合Ｃｉ（ｉ＝１，…，Ｉ）を次式（４３）により更新する。
【０１０４】
【数４３】
2010020258.tif 000045
【０１０５】
これを用いて、次式（４４）により、ｍ１，…，ｍＩを求める。但し、｜Ｃｉ｜、｜εｎ｜は、それぞれ集合Ｃｉ、εｎの要素数を表す。以上より求めたｑ１，…，ｑＮとｍ１，…，ｍＩを用いて、音高目標値時系列ｍｑ１，ｍｑ２，…，ｍｑＮが求まる。
【０１０６】
【数４４】
2010020258.tif 000046
【０１０７】
Ａ−６．収束判定過程
反復計算が所定の回数を満たしたか否か、あるいは、反復計算においてパラメータの更新の変化率が所定値以下になったか否か、あるいは、目的関数値の変化率が所定値以下になったか否かを判定する。
【０１０８】
Ｂ．歌唱の基本周波数制御モデル
歌唱の基本周波数（Ｆ０）軌跡は、歌唱者が頭の中で思い描く旋律（ステップ状成分からなる階段関数）に、歌唱者の表現意図や身体的特性による動的変動成分が付加されて観測される。本発明は、観測されるＦ０軌跡から、歌唱者が歌おうとする旋律概形と歌唱者ごとに異なる動的変動成分とを分離する手法に関する。
【０１０９】
Ｂ−１．２階の微分方程式に基づくＦ０制御
Ｆ０軌跡の動的特性は、次式（４５）で示される、２次系の制御システムを用いて良く表現できる。
【０１１０】
【数４５】
2010020258.tif 000047
【０１１１】
ここで、歌唱者が頭の中で思い描く旋律概形を、ステップ状成分からなる階段関数ｕ（ｔ）と表現する。このｕ（ｔ）を入力としたときに、数式（４５）に従ってＦ０軌跡ｙ（ｔ）が観測されるものとする。はじめに、このシステムの伝達関数Ｇ（ｓ）を求め、その特性を確認する。まず、数式（４５）の両辺をラプラス変換する。Ｙ（ｓ）＝Ｌ［ｙ（ｔ）］とすると、関数ｙ（ｔ）のｎ階導関数ｙ（ｎ）（ｔ）のラプラス変換は、次式（４６）と書ける。
【０１１２】
【数４６】
2010020258.tif 000048
【０１１３】
ｙ（０）＝０、ｕ（０）＝０とし、数式（４５）の両辺をラプラス変換すると、次式（４７）で表される。
【０１１４】
【数４７】
2010020258.tif 000049
【０１１５】
伝達関数は、Ｇ（ｓ）＝Ｙ（ｓ）／Ｕ（ｓ）で表されるので、次式（４８）となる。
【０１１６】
【数４８】
2010020258.tif 000050
【０１１７】
一方、制動２次系伝達関数Ｈ（ｓ）は、一般的に、次式（４８）のように記述される。
【０１１８】
【数４９】
2010020258.tif 000051
【０１１９】
この伝達関数のステップ応答を図１に示す。ζ＝０では、定常振動となり、これは音高安定時間が持続した場合に観測される４〜７［Ｈｚ］の周期的な振動であるビブラートに対応すると考える。さらに、（｜ζ｜＜１）では減衰振動となり、これは音高遷移時に目的音高より大きく振れてしまうオーバーシュートに対応すると考える。このように、音高遷移における基本周波数の挙動を２次系の制御システムは適切に表現できる。また、パラメータΩ、ζ、Ｋと数式（４５）の微分方程式の係数α、β、γとの関係は、次式（５０）の通りである。
【０１２０】
【数５０】
2010020258.tif 000052
【０１２１】
Ｂ−２．差分方程式への変換
連続時間信号を扱う数式（４５）の微分方程式を、離散時間信号を扱う差分方程式に変換する。連続時間信号ｙ（ｔ）が、Ｎｙｑｕｉｓｔ周波数より低い帯域制限を受けた信号であると仮定すると、このｙ（ｔ）は、対応する離散時間信号ｙｎと連続時間信号ｓｉｎｃ（πｔ／△）の離散畳込みで、次式（５１）に示すように表現することができる。
【０１２２】
【数５１】
2010020258.tif 000053
【０１２３】
ここで、△はサンプリング周期とする。数式（５１）よりｙ（ｔ）の１階微分、２階微分は、それぞれ数式（５２）、（５３）となる。
【０１２４】
【数５２】
2010020258.tif 000054
【０１２５】
【数５３】
2010020258.tif 000055
【０１２６】
ここで、ｔ＝ｎ△における１階微分は、次式（５４）で表わされる。
【０１２７】
【数５４】
2010020258.tif 000056
【０１２８】
観測Ｎ点だけで上記微分の近似を行うと、次式（５５）となる。
【０１２９】
【数５５】
2010020258.tif 000057
【０１３０】
同様に、ｔ＝ｎ△における２階微分は、次式（５６）で表わされる。
【０１３１】
【数５６】
2010020258.tif 000058
【０１３２】
これを同様に観測Ｎ点だけで近似を行うと、次式（５７）となる。
【０１３３】
【数５７】
2010020258.tif 000059
【０１３４】
数式（５５）と数式（５７）を数式（４５）の微分方程式に代入すると、次式（５８）で表わされる。
【０１３５】
【数５８】
2010020258.tif 000060
【０１３６】
上記数式（５８）は、ｙ＝（ｙ１，ｙ２，…，ｙＮ）Ｔ、ｕ＝（ｕ１，ｕ２，…，ｕＮ）Ｔと置くと、次式（５９）に示すように簡潔に表現できる。
【０１３７】
【数５９】
2010020258.tif 000061
【０１３８】
但し、行列Ｂ、Ａは、次式（６０）の通りである。
【０１３９】
【数６０】
2010020258.tif 000062
【数６１】
2010020258.tif 000063
【０１４０】
数式（５９）の線形方程式のパラメータは、実際には、α、β、γの３つ（行列Ｂ、Ａは定数行列）であるが、次式（６１）のように置き、行列Ｗ＝（Ｗｉ，ｊ）の要素が、ＷがＴｏｅｐｌｉｔｚ型であるという拘束つきの自由パラメータであるような状況も併せて考えることにする（数式（５９）が成立するための必要条件は、ＷがＴｏｅｐｌｉｔｚ型であることによる）。すなわち、α、β、γを推定すべきパラメータとする問題（以降、「ケース１」）、及び、行列Ｗの要素Ｗｉ，ｊをＴｏｅｐｌｉｔｚ行列の拘束条件の下で推定すべきパラメータとする問題（以後、「ケース２」）を、以後同時並行的に検討する。
【０１４１】
Ｂ−３．階段関数ｕｎのモデル化
システムの入力となる旋律概形を表す階段関数ｕｎは、図２に示すような状態集合Ｓ＝｛Ｓ１，…，ＳＩ｝からなるＨＭＭ（隠れマルコフモデル）を利用して、以下のようにモデル化する。すなわち、１回の状態遷移によりサンプル点が１個生成されるモデルであり、状態によって出力の統計的な傾向が異なる。ここでは、式の見やすさのため、Ｓｉ＝ｉとし、次式（６２）で表わす。
【０１４２】
【数６２】
2010020258.tif 000064
【０１４３】
すなわち、数式（６３）、（６４）で表わされる。
【０１４４】
【数６３】
2010020258.tif 000065
【０１４５】
【数６４】
2010020258.tif 000066
【０１４６】
ｍＳｉは、状態Ｓｉにおける出力確率分布（正規分布）の平均を表す。一様なマルコフ連鎖を想定し、状態Ｓｊから状態Ｓｉへの遷移確率は、＊Ｐ（Ｓｉ｜Ｓｊ）と表す。ここで、遷移確率＊Ｐ（Ｓｉ｜Ｓｊ）は、定数とする。前述したＰｉ，ｊは、ｌｏｇ＊Ｐ（Ｓｉ｜Ｓｊ）のことである。同じ状態へ遷移する遷移確率＊Ｐ（Ｓｉ｜Ｓｊ）を自己遷移確率といい、これが大きいほど同じ状態に留まろうとする傾向が強くなる。ｑｎは、状態集合Ｓの要素（ＨＭＭの状態番号）の中のいずれかの値をとる。従って、状態系列ｑ１，…，ｑＮと、各時刻の状態における出力確率分布の平均ｍｑｎによって階段関数が決定される。
【０１４７】
Ｂ−４．Ｆ０制御モデルのパラメータの解釈
上述では、２階の微分方程式を利用したＦ０制御モデルと、入力となる階段階数とをモデル化した。以上より、ケース１では、次式（６５）が、ケース２では、次式（６６）が推定したいパラメータである。
【０１４８】
【数６５】
2010020258.tif 000067
【０１４９】
【数６６】
2010020258.tif 000068
【０１５０】
ここで、各パラメータメータの解釈を図３を参照して説明する。まず、微分方程式の係数α、β、γ、あるいは（Ｗｉ，ｊ）は、歌声の動的変動を表すパラメータである。音高が安定するときの振動であったり、音高遷移における連続的なダイナミクスを表現する。状態系列ｑ１，…，ｑＮは、音高が安定する長さを決定するパラメータである。これは、必ずしも楽譜に記される音符の長さに対応するわけではなく、歌唱者の意図や、歌唱スタイルに基づいて生成される運動指令の長さを表現したものであると、ここでは想定している。最後に、ＨＭＭの各状態の平均ｍＳ１，…，ｍＳＩは、歌唱者が意図する旋律の音高（音高目標値）に対応するパラメータである。これは、必ずしも楽譜に記される音符の音高（客観的に定まっている音高値）には対応しない。
【０１５１】
Ｃ．Ｆ０制御モデルのパラメータ最尤推定
微分方程式の係数α、β、γと、階段関数をモデル化する状態系列ｑ１，…，ｑＮ、各状態の正規分布の平均ｍＳｌ，…，ｍＳＩを、観測系列ｙ＝（ｙ１，…，ｙＮ）Ｔから最尤推定する方法について述べる。
【０１５２】
観測系列が線形差分方程式に理想的に従うならば、数式（５９）を解けば良いが、実際には、理想的な差分方程式からの誤差があると考えられる。そこで、次式（６７）と次式（６８）との間に、数式（６９）が成り立つと仮定する。
【０１５３】
【数６７】
2010020258.tif 000069
【０１５４】
【数６８】
2010020258.tif 000070
【０１５５】
【数６９】
2010020258.tif 000071
【０１５６】
但し、行列Ｗ、Ｂ、Ａ、ｍは、それぞれ、次式（７０）、（７１）、（７２）で表わされる。
【０１５７】
【数７０】
2010020258.tif 000072
【０１５８】
【数７１】
2010020258.tif 000073
【０１５９】
【数７２】
2010020258.tif 000074
【０１６０】
ここで、εの要素εｎは平均０、分散σ２の正規分布に従う互いにＧａｕｓｓ性白色雑音である。いま、観測系列ｙに対するパラメータΘの尤度は、数式（６９）より、次式（７３）となるから、多次元正規分布であることが分かり、その正規化係数は、次式（７４）とすればよく、結局、次式（７５）のような形となる。
【０１６１】
【数７３】
2010020258.tif 000075
【０１６２】
【数７４】
2010020258.tif 000076
【０１６３】
【数７５】
2010020258.tif 000077
【０１６４】
このとき、パラメータΘの対数尤度は、次式（７６）となる。
【０１６５】
【数７６】
2010020258.tif 000078
【０１６６】
パラメータΘの事後確率は、次式（７７）であるので、ここで、次式（７８）が成り立つとすると、次式（７９）と表すことができる。
【０１６７】
【数７７】
2010020258.tif 000079
【０１６８】
【数７８】
2010020258.tif 000080
【０１６９】
【数７９】
2010020258.tif 000081
【０１７０】
ここで、事前確率＊Ｐ（α，β，γ）と＊Ｐ（ｍＳ１，…，ｍＳＩ）は一様分布とし、＊Ｐ（ｑ１，…，ｑＮ）は、先に述べたようにＨＭＭによる一様なマルコフ連鎖を想定している。このため、次式（８０）で表わすことができる。
【０１７１】
【数８０】
2010020258.tif 000082
【０１７２】
遷移確率（Ｐ）（Ｓｉ｜Ｓｊ）は、事前に決定する定数である。以後、簡単のため、次式（８１）で表わすことにする。したがって、次式（８２）となる。
【０１７３】
【数８１】
2010020258.tif 000083
【０１７４】
【数８２】
2010020258.tif 000084
【０１７５】
以上より、数式（７６）と数式（８０）を、数式（７９）に代入し、定数項を除いた次式（８３）がパラメータΘに関して最大化したい目的関数である。
【０１７６】
【数８３】
2010020258.tif 000085
【０１７７】
しかし、数式（８３）を最大化するパラメータΘは、解析的に求めることができない。そこで、Θの各要素（微分方程式の係数、ガウス性雑音の分散、状態系列、状態の出力分布の平均）に関して、他の要素を固定した下で、数式（８３）を最大化するステップを、数式（８３）の値が収束するまで繰り返す。
【０１７８】
Ｃ−１．特徴パラメータ行列Ｗの更新
状態系列ｑ１，…，ｑＮ、ＨＭＭの各状態の正規分布の平均値ｍＳ１，…，ｍＳＩを固定したとき、数式（８３）を最大にするＷを求めたい。ここでは、これを実現する３つの解法について説明する。
【０１７９】
Ｃ−１−１．解法１（微分方程式逆フィルタを用いた方法）
ここでは、ケース１（自由パラメータはα、β、γ）を想定する。上記数式（８３）の右辺の第２項は、第３項に比べてＪへの寄与が無視できるほど小さいと仮定し、次式（８４）を最小化するα、β、γが、Ｊを最大化するものと近似的に見なす。
【０１８０】
【数８４】
2010020258.tif 000086
【０１８１】
Ｊ１をαに関して偏微分すると、次式（８５）を得る。
【０１８２】
【数８５】
2010020258.tif 000087
【０１８３】
これを０と置くと、次式（８６）となる。
【０１８４】
【数８６】
2010020258.tif 000088
【０１８５】
また、同様に、Ｊ１をβ及びγに関して偏微分して０と置くと、次式（８７）、（８８）を得る。
【０１８６】
【数８７】
2010020258.tif 000089
【０１８７】
【数８８】
2010020258.tif 000090
【０１８８】
以上より立てられる次式（８９）の正規方程式を解けばよい。
【０１８９】
【数８９】
2010020258.tif 000091
【０１９０】
すなわち、次式（９０）がα、β、γの更新値となる。
【０１９１】
【数９０】
2010020258.tif 000092
【０１９２】
また、以上により求めたα、β、γを用い、雑音の分散推定値が次式（９１）により求まる。
【０１９３】
【数９１】
2010020258.tif 000093
【０１９４】
Ｃ−１−２．解法２（逆フィルタを用いた方法）
ここでも、ケース１（自由パラメータはα、β、γ）を想定する。解放１のように、数式（８３）の第２項（ｌｏｇ｜Ｗ｜）の寄与を無視しないとすると、Ｊを最大化するα、β、γは解析的に求まらない。そこで、ここでは、α、β、γの更新値を勾配法により数値計算する方法について説明する。以下で与えられる、Ｊのα、β、γに関する勾配ベクトル▽Ｊにより、最急降下法、共役勾配法、準ニュートン法などが適用できる。
【０１９５】
【数９２】
2010020258.tif 000094
【０１９６】
【数９３】
2010020258.tif 000095
【０１９７】
【数９４】
2010020258.tif 000096
【０１９８】
【数９５】
2010020258.tif 000097
また、勾配法により、α、β、γを更新するごとに、雑音の分散推定値σ２を、次式（９６）により更新する。
【０１９９】
【数９６】
2010020258.tif 000098
【０２００】
Ｃ−１−３．解法３（差分方程式逆フィルタを用いた方法）
ここでは、ケース２（自由パラメータは行列Ｗの要素）を想定する。行列Ｗを、以下のように対角成分が、次式（９７）で示すように、全て１のＴｏｅｐｌｉｔｚ型の上三角行列と仮定し、α、β、γの代わりに、ｗ１，…，ｗＭ（Ｍ≦Ｎ−１）をパラメータとして求める。但し、ｗｍ＝０（ｍ≧Ｍ＋１）とする。
【０２０１】
【数９７】
2010020258.tif 000099
【０２０２】
すなわち、ここでは、行列Ｗに関して数式（７０）のような要素に関する拘束を仮定しない代わりに、数式（９７）のような構造の拘束を与える。
【０２０３】
このとき、｜Ｗ｜＝１であることに注意すると、ｌｏｇ｜Ｗ｜＝０であるので、次式（９８）を最小化するｗ１，…，ｗＭを求めればよい。
【０２０４】
【数９８】
2010020258.tif 000100
【０２０５】
行列ＷはＴｏｅｐｌｉｔｚ行列のため、ベクトルｗ＝（ｗ１，…，ｗＮ−１）Ｔとすると、次式（９９）で表わされる。
【０２０６】
【数９９】
2010020258.tif 000101
【０２０７】
但し、行列Ｕは、次式（１００）で表わされる。
【０２０８】
【数１００】
2010020258.tif 000102
【０２０９】
したがって、Ｊ１は、次式（１０１）で表わされる。
【０２１０】
【数１０１】
2010020258.tif 000103
【０２１１】
これを行列ｗに関して偏微分して０と置くと、次式（１０２）で表わされる正規方程式を得る。
【０２１２】
【数１０２】
2010020258.tif 000104
【０２１３】
これを解くと、次式（１０３）を得る。
【０２１４】
【数１０３】
2010020258.tif 000105
【０２１５】
以上により求まった行列ｗの要素を数式（９７）に基づいて、行列Ｗの中に代入すれば、Ｔｏｅｐｌｉｔｚ型の拘束条件の下での最適な行列Ｗを求めたことになる。また、以上により求まった行列Ｗを用い、雑音の分散推定値が次式（１０４）により求まる。
【０２１６】
【数１０４】
2010020258.tif 000106
【０２１７】
Ｃ−２．状態系列ｑ１，…，ｑＮの推定
特徴パラメータ行列Ｗ、ＨＭＭの各状態における正規分布の平均値ｍＳ１，…，ｍＳＩを固定したとき、状態系列ｑ１，…，ｑＮに関して数式（８３）の最大化を考える。つまり、数式（８３）から関係する項だけを取り出して、次式（１０５）と置き、これを最大にする最適な状態系列ｑ１，…，ｑＮを求める。但し、次式（１０６）が成立するものとする。
【０２１８】
【数１０５】
2010020258.tif 000107
【０２１９】
【数１０６】
2010020258.tif 000108
【０２２０】
この問題は、ビタビ（Ｖｉｔｅｒｂｉ）アルゴリズム（動的計画法）により効率的に解くことができる。まず、次のような、最初から時刻ｋに状態Ｓｉに至るまでの部分系列に関する最適な状態系列について、次式（１０７）のような量が定義されているとする。
【０２２１】
【数１０７】
2010020258.tif 000109
【０２２２】
この量は、漸化式で、次式（１０８）のようにして得られる。
【０２２３】
【数１０８】
2010020258.tif 000110
【０２２４】
これをｋ＝Ｎまで計算すれば、最適経路が求められる。
図４にアルゴリズムを示す。図において、まず、次式（１０９）に従って初期化を行う（ステップＳｃ１）。
【０２２５】
【数１０９】
2010020258.tif 000111
【０２２６】
次に、次式（１１０）に従って再帰処理を行う（ステップＳｃ２）。
【０２２７】
【数１１０】
2010020258.tif 000112
【０２２８】
次に、次式（１１１）に従って終了処理を行う（ステップＳｃ３）。
【０２２９】
【数１１１】
2010020258.tif 000113
【０２３０】
そして、次式（１１２）で示されるように、状態系列のバックトラックを行う（ステップＳｃ４）。
【０２３１】
【数１１２】
2010020258.tif 000114
【０２３２】
Ｃ−３．ＨＭＭの各状態の出現確率分布（正規分布）の平均ｍＳ１，…，ｍＳＩの推定
特徴パラメータ行列Ｗと状態系列ｑ１，…，ｑＮを固定したとき、数式（８３）が最大となるように、ＨＭＭの各状態における正規分布の平均ｍＳ１，…，ｍＳＩを更新する。つまり、数式（８３）からｍｓ１，…，ｍＳＩに関係する項だけを取り出して符号を反転した、次式（１１３）で示される値の最小化を考えればよい。
【０２３３】
【数１１３】
2010020258.tif 000115
【０２３４】
すなわち、次式（１１４）を解くと、次式（１１５）を得る。
【０２３５】
【数１１４】
2010020258.tif 000116
【０２３６】
【数１１５】
2010020258.tif 000117
【０２３７】
ここで、集合Ｃｉ＝｛ｎ｜ｑｎ＝Ｓｉ｝とし、｜Ｃｉ｜をその要素数とする。数式（１１５）より、各状態における正規分布の平均ｍＳｉが更新される。以上は、前述したビタビアルゴリズムによって求められた最適な状態系列を用いて、ＨＭＭの各状態のパラメータを学習する操作であるため、しばしばビタビ学習と呼ばれる。最終的に状態系列ｑ１，…，ｑＮとＨＭＭの各状態の正規分布の平均値ｍＳ１，…，ｍＳＩから、次式（１１６）で示される入力の階段階数ｕｎ、すなわち、行列ｍが求まる。
【０２３８】
【数１１６】
2010020258.tif 000118
【０２３９】
Ｃ−４．初期値設定
初期値設定では、Ｂ−１、Ｂ−２、Ｂ−３で説明した３段階のパラメータ推定を、パラメータΘの対数事後確率が収束するまで順番に繰り返す。しかし、以上の反復法によるパラメータ推定は、初期値を適切に設定しないと、局所解に収束してしまう。この問題に対処するため、推定手順の前半に２つの初期値設定、すなわち第１の初期値設定及び第２の初期値設定を行う。
【０２４０】
第１の初期値設定では、観測系列ｙに、Ｂ−２で説明したビタビアルゴリズムを適用する。ここでは、次式（１１７）に示す値を最小化する状態系列ｑ１，…，ｑＮを、前述したＢ−２と同様の方法で求め、これを次式（１１８）で示す初期状態系列とする。
【０２４１】
【数１１７】
2010020258.tif 000119
【０２４２】
【数１１８】
2010020258.tif 000120
【０２４３】
第２の初期値設定では、第１の初期値設定で求めた初期状態系列をもとにＨＭＭの各状態の正規分布の平均を求める。すなわち、次式（１１９）、（１２０）を解き、次式（１２１）を得て、これを次式（１２２）とする。
【０２４４】
【数１１９】
2010020258.tif 000121
【０２４５】
【数１２０】
2010020258.tif 000122
【０２４６】
【数１２１】
2010020258.tif 000123
【０２４７】
【数１２２】
2010020258.tif 000124
【０２４８】
以上のように、初期値設定の段階で、観測系列ｙから、ある程度の階段関数ｕｎの概形を推定することにより、局所解に収束してしまうことを防ぐことが可能となる。
【０２４９】
Ｄ．時変なＦ０制御モデルへの拡張
前述では、特徴パラメータ行列Ｗは、時不変なものとして推定を行ったが、観測系列フレームと呼ぶ区間に分割し（但し、区間は重複してもよい）、α、β、γ、あるいは特徴パラメータ行列Ｗを、フレーム毎に自由度もつパラメータと見なして推定するようにしてもよい。
【０２５０】
（実施例）
図５は、本実施形態による、基本周波数軌跡モデルパラメータ抽出装置の構成を示すブロック図である。図において、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列を抽出する。初期音高目標値生成部２は、抽出された観測基本周波数時系列を入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列を生成する。フレーム分割部３は、観測基本周波数時系列をフレーム分割する。
【０２５１】
逆フィルタ出力値系列更新部４は、初期音高目標値時系列とフレーム分割された観測基本周波数系列とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、観測基本周波数系列から逆フィルタ出力値系列と特性パラメータとを生成する。音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列と観測基本周波数系列とを入力として、音高目標値時系列を生成する。
【０２５２】
収束判定部６は、特性パラメータと音高目標値時系列とが所定の規準を満たしているかどうかを判定し、満たしていない場合には、逆フィルタ出力値系列更新部４と音高目標値更新部５とに再処理させる。基本周波数軌跡パラメータ出力部７は、収束判定部６において所定の規準を満たすと判定された特性パラメータと音高目標値時系列とを出力する。
【０２５３】
図６は、本実施形態において、逆フィルタ出力値系列更新部４で、前述した逆フィルタを用いる第１動作例（Ｃ−１−２．解法２）を説明するためのフローチャートである。まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＡ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍｉｎｉｔを生成する（ステップＳＡ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）を出力する（ステップＳＡ３）。
【０２５４】
次に、逆フィルタ出力値系列更新部４は、初期音高目標値時系列ｍｉｎｉｔとフレーム分割された観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）とを入力として、フレーム毎の音高目標値時系列と観測基本周波数系列との関係を表すフィルタの特性パラメータを更新し、観測基本周波数系列から逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と特性パラメータＷ（１），…，Ｗ（Ｆ）とを生成する（ステップＳＡ４）。次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と観測基本周波数系列ｙ（１），…，ｙ（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＡ４）。
【０２５５】
収束判定部６は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＡ６）、満たしていない場合には（ステップＳＡ６のＮＧ）、ステップＳＡ４に戻り、逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＡ６の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＡ７）。
【０２５６】
次に、図７は、本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−１、Ａ−４−２で説明した微分方程式逆フィルタを用いる第１動作例（Ｃ−１−１．解法１）を説明するためのフローチャートである。
まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＢ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍｉｎｉｔを生成する（ステップＳＢ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）を出力する（ステップＳＢ３）。
【０２５７】
次に、逆フィルタ出力値系列更新部４は、前述したＡ−４−１で説明したように、適当な１階数値微分係数と２階数値微分係数を用いて、数値微分フィルタ係数Ａ、Ｂを算出し（ステップＳＢ４）、該数値微分フィルタ係数Ａ、Ｂと、初期音高目標値時系列ｍｉｎｉｔと、フレーム分割された観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）とを入力として、前述したＡ−４−２で説明したように、適当な１階数値微分係数と２階数値微分係数を用いて、フレーム毎の逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と特性パラメータＷ（１），…，Ｗ（Ｆ）とを生成する（ステップＳＢ４）。
【０２５８】
次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と観測基本周波数系列ｙ（１），…，ｙ（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＢ５）。収束判定部６は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＢ７）、満たしていない場合には（ステップＳＢ７のＮＧ）、ステップＳＢ５に戻り、微分方程式逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＢ７の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＢ８）。
【０２５９】
次に、図８は、本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−３で説明した差分方程式逆フィルタを用いる第２動作例（Ｃ−１−３．解法３）を説明するためのフローチャートである。まず、基本周波数抽出部１は、入力される音響信号から観測基本周波数時系列ｙを抽出する（ステップＳＣ１）。次に、初期音高目標値生成部２は、抽出された観測基本周波数時系列ｙを入力として、値が有限個の離散的な実数値の中からしか取り得ない拘束を有する音高目標値の初期値の時系列、すなわち、初期音高目標値系列ｍｉｎｉｔを生成する（ステップＳＣ２）。フレーム分割部３は、観測基本周波数時系列をフレーム分割し、観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）を出力する（ステップＳＣ３）。
【０２６０】
次に、逆フィルタ出力値系列更新部４は、初期音高目標値時系列ｍｉｎｉｔと、フレーム分割された観測基本周波数系列У（１），…，У（Ｆ）、ｙ（１），…，ｙ（Ｆ）とを入力として、前述したＡ−４−３で説明したように、フレーム毎の逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と特性パラメータＷ（１），…，Ｗ（Ｆ）とを生成する（ステップＳＣ４）。
【０２６１】
次に、音高目標値更新部５は、全てのフレームの逆フィルタ出力値系列ｕ＾（１），…，ｕ＾（Ｆ）と観測基本周波数系列ｙ（１），…，ｙ（Ｆ）とを入力として、音高目標値時系列ｍを生成する（ステップＳＣ５）。収束判定部６は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとが所定の規準を満たしているかどうかを判定し（ステップＳＣ６）、満たしていない場合には（ステップＳＣ６のＮＧ）、ステップＳＣ４に戻り、差分方程式逆フィルタ出力値系列更新と音高目標値更新とを再実行する。一方、収束判定部６において所定の規準を満たすと判定された場合には（ステップＳＣ６の収束）、基本周波数軌跡パラメータ出力部７は、特性パラメータＷ（１），…，Ｗ（Ｆ）と音高目標値時系列ｍとを出力する（ステップＳＣ７）。
【０２６２】
上述した実施形態によれば、歌声知覚に対しての妥当性が確認された２次系歌声生成モデルの制御パラメータを実データから教師なし学習する枠組みが確立されれば、歌声合成技術が飛躍的に向上させることができる。例えば、従来不可能であった歌唱スタイルの転写が実現される。歌唱者Ａの歌い方を制御パラメータから学習し、旋律ａを表す階段状の信号に畳み込むことによって、歌唱者Ａの歌い方による旋律ａの歌声を自動生成できる。
【０２６３】
また、上述した実施形態によれば、表現豊かな歌い方を制御パラメータによって予め学習することで、高品質な歌声合成が期待される。
【０２６４】
さらに、カラオケ採点システムへの応用も考えられる。これまで、歌唱力の自動評価手法が提案されているが、聴取実験に基づいて様々な評価尺度を提案し、歌の上手・下手の２クラス識別を行っているものの、歌い方を精密にモデル化するまでには至っていない。これに対して、上述した実施形態によれば、制御パラメータを評価尺度に利用することで、性能向上を図ることができる。
【０２６５】
また、上述した実施形態においては、基本周波数抽出部１、初期音高目標値生成部２、フレーム分割部３、逆フィルタ出力値系列更新部４、音高目標値更新部５、収束判定部６、基本周波数軌跡パラメータ出力部７は、プログラムの形式でコンピュータ読み取り可能な記録媒体に記憶されていてもよく、このプログラムをコンピュータが読み出して実行することによって、上記処理を行うようにしてもよい。すなわち、放送装置５における、各処理手段、処理部は、ＣＰＵ等の中央演算処理装置がＲＯＭやＲＡＭ等の主記憶装置に上記プログラムを読み出して、情報の加工・演算処理を実行することにより、実現するようにしてもよい。
【０２６６】
ここでコンピュータ読み取り可能な記録媒体とは、磁気ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、半導体メモリ等をいう。また、このコンピュータプログラムを通信回線によってコンピュータに配信し、この配信を受けたコンピュータが当該プログラムを実行するようにしても良い。
【図面の簡単な説明】
【０２６７】
【図１】２次系伝達関数のステップ応答を示す概念図である。
【図２】階段関数ｕｎのモデル化を説明するための概念図である。
【図３】Ｆ０制御パラメータの解釈を説明するための概念図である。
【図４】最適な状態系列を求めるためのビタビアルゴリズムの一例を示す概念図である。
【図５】本実施形態による、基本周波数軌跡モデルパラメータ抽出装置の構成を示すブロック図である。
【図６】本実施形態において、逆フィルタ出力値系列更新部４で、前述した逆フィルタを用いる第１動作例（Ｃ−１−２．解法２）を説明するためのフローチャートである。
【図７】本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−１、Ａ−４−２で説明した微分方程式逆フィルタを用いる第１動作例（Ｃ−１−１．解法１）を説明するためのフローチャートである。
【図８】本実施形態において、逆フィルタ出力値系列更新部４で、前述したＡ−４−３で説明した差分方程式逆フィルタを用いる第２動作例（Ｃ−１−３．解法３）を説明するためのフローチャートである。
【符号の説明】
【０２６８】
１  基本周波数抽出部
２  初期音高目標値生成部
３  フレーム分割部
４  逆フィルタ出力値系列更新部
５  音高目標値更新部
６  収束判定部
７  基本周波数軌跡パラメータ出力部
【図１】
2010020258.tif 000127
【図２】
2010020258.tif 000128
【図３】
2010020258.tif 000129
【図４】
2010020258.tif 000130
【図５】
2010020258.tif 000131
【図６】
2010020258.tif 000132
【図７】
2010020258.tif 000133
【図８】
2010020258.tif 000134