(57)【特許請求の範囲】
【請求項１】
観測信号を、入力信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる信号生成系の出力信号と残差信号との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータと、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータと、前記残差信号を表すモデルを構成する残差信号パラメータとを推定することにより時系列信号の動特性特徴を抽出する信号解析装置であって、
前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成部と、
前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータとし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離部と、
前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られる関数を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新部と、
前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離部と前記モデルパラメータ更新部による処理を再度行わせるパラメータ収束判定部と、
前記パラメータ収束判定部により前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力部と、
を備えることを特徴とする信号解析装置。
【請求項２】
前記入力信号はステップ信号であり、
前記出力信号は多次元ガウス分布に従うものとして確率的にモデル化され、
前記残差信号は、ガウス性白色雑音として確率的にモデル化される
ことを特徴とする請求項１記載の信号解析装置。
【請求項３】
前記信号生成系のフィルタ特性は、差分法によって導出されるフィルタで表され、
前記フィルタ特性パラメータは、前記固有周波数の二乗に反比例するパラメータと、前記減衰率に比例し前記固有周波数に反比例するパラメータとである
ことを特徴とする請求項２記載の信号解析装置。
【請求項４】
前記信号生成系のフィルタ特性は、自己回帰過程に基づいて構成されるフィルタで表され、前記フィルタ特性パラメータは、自己回帰パラメータであることを特徴とする請求項２記載の信号解析装置。
【請求項５】
前記信号生成系のフィルタ特性は、複数の２次系フィルタの重み付き線形和によって構成されるフィルタで表され、前記フィルタ特性パラメータは、前記各２次系フィルタの重みであることを特徴とする請求項２記載の信号解析装置。
【請求項６】
前記モデルパラメータ更新部は、
補助変数から構成される前記目的関数の補助関数を、前記固有周波数の二乗に反比例するパラメータと前記補助関数を前記減衰率に比例し前記固有周波数に反比例するパラメータとでそれぞれ微分して得られる方程式からなる連立方程式をフィルタ特性パラメータについて解くことにより、フィルタ特性パラメータの値を更新するフィルタ特性パラメータ更新部と、
前記補助関数を、前記入力信号パラメータで微分して得られる方程式を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、
前記補助関数を、前記残差信号パラメータで微分して得られる方程式を解くことにより、残差信号パラメータを更新する残差信号パラメータとから構成されることを特徴とする請求項３記載の信号解析装置。
【請求項７】
前記モデルパラメータ更新部は、
前記目的関数を、前記フィルタ特性パラメータに含まれる各自己回帰パラメータでそれぞれ微分して得られる方程式からなる連立方程式をフィルタ特性パラメータについて解くことにより、フィルタ特性パラメータの値を更新するフィルタ特性パラメータ更新部と、
前記目的関数を、前記入力信号パラメータで微分して得られる方程式を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、
前記目的関数を、前記残差信号パラメータで微分して得られる方程式を解くことにより、残差信号パラメータを更新する残差信号パラメータとから構成されることを特徴とする請求項４記載の信号解析装置。
【請求項８】
前記モデルパラメータ更新部は、
補助変数から構成される前記目的関数の補助関数を、前記フィルタ特性パラメータである各２次系フィルタの重みでそれぞれ微分して得られる方程式からなる非線形連立方程式を、フィルタ特性パラメータについて解くことにより、フィルタ特性パラメータの値を更新するフィルタ特性パラメータ更新部と、
前記補助関数を、前記入力信号パラメータで微分して得られる方程式を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、
前記補助関数を、前記残差信号パラメータで微分して得られる方程式を解くことにより、残差信号パラメータを更新する残差信号パラメータとから構成されることを特徴とする請求項５記載の信号解析装置。
【請求項９】
前記信号分離部は、前記観測信号と前記モデルパラメータが与えられた時の前記出力信号および前記残差信号から構成される完全データの期待値と、前記完全データの自己相関とを用いて、前記観測信号を出力信号と残差信号とに分離することを特徴とする請求項２から８のいずれかに記載の信号解析装置。
【請求項１０】
観測信号を、入力信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる信号生成系の出力信号と残差信号との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータと、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータと、前記残差信号を表すモデルを構成する残差信号パラメータとを推定することにより時系列信号の動特性特徴を抽出する信号解析装置における信号解析方法であって、
前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成ステップと、
前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータとし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離ステップと、
前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られる関数を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新ステップと、
前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離ステップと前記モデルパラメータ更新ステップによる処理を再度行わせるパラメータ収束判定ステップと、
前記パラメータ収束判定ステップにより前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力ステップと、
を有することを特徴とする信号解析方法。
【請求項１１】
観測信号を、入力信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる信号生成系の出力信号と残差信号との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータと、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータと、前記残差信号を表すモデルを構成する残差信号パラメータとを推定することにより時系列信号の動特性特徴を抽出する信号解析装置上のコンピュータに信号解析を行わせる信号解析プログラムであって、
前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成ステップと、
前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータとし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離ステップと、
前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られる関数を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新ステップと、
前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離ステップと前記モデルパラメータ更新ステップによる処理を再度行わせるパラメータ収束判定ステップと、
前記パラメータ収束判定ステップにより前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力ステップと、
を前記コンピュータに行わせることを特徴とする信号解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時系列信号の動特性特徴を抽出する信号解析装置、信号解析方法及び信号解析プログラムに関する。
【背景技術】
【０００２】
歌声音響信号から抽出される基本周波数系列を例に挙げて、従来技術を説明する。この歌声の基本周波数系列には、歌唱者が歌おうとする音高目標値系列と歌唱力・歌唱スタイル・個人性・感情に基づく様々な動的変動成分（オーバーシュートやビブラートなど）が複雑に重ね合わされている。歌声は、多くのジャンルの音楽を特徴づける重要な要素の一つであり、この歌声の基本周波数系列に着目した様々な研究が現在盛んに行われている。
【０００３】
基本周波数系列に含まれる音高目標値系列を特徴抽出できれば、ハミング検索や自動採譜への応用が期待される。特に、歌声から楽曲を検索するハミング検索では、歌唱された歌声の基本周波数系列から、歌唱者の意図する音高目標値系列を正しく推定して、楽曲データベースの旋律と照合する必要がある。
【０００４】
一方で、オーバーシュートやビブラートのような基本周波数系列の動的変動成分は、歌声知覚、個人性知覚に影響を与える成分であることが知られている。したがって、歌唱スタイルの記述やそれを利用した類似歌声検索、歌唱力自動評価のための有用な尺度となりうる。また、より表情豊かかつ多様な歌声合成のためにも必要不可欠な成分である。そこで従来研究では、線形２次系を利用して歌声の基本周波数の動的変動成分を制御するモデルが提案されている（例えば、非特許文献１、２、３参照）。
【０００５】
これらの研究では、日本語の話声の基本周波数パターンを表現する藤崎モデルが参考にされた。藤崎モデルは、臨界制動２次系のインパルス応答とステップ応答を利用して、日本語の句頭から句末に向けて緩やかに下降するフレーズ成分と、語句に対応して急激に上昇下降するアクセント成分を表現し、これらを重畳することで、基本周波数系列を記述する。ただし、歌声の旋律に伴った急激な基本周波数の上昇・下降の制御及び、ビブラートのような周期的な振動は、臨界制動系では表現できない。そのため、歌声の基本周波数制御モデルでは２次系の伝達関数
【数１】
0005318042.tif 000002
における減衰率ζを調整することによって、指数減衰（ζ＞１）、減衰振動（０＜ζ＜１、オーバーシュートに対応する）、臨界制動（ζ＝１）、定常振動（ζ＝０、ビブラートに対応する）からなる様々な振動現象を表現する。
【０００６】
非特許文献３では、音階を表す階段状信号に式（１）のインパルス応答を畳み込んで得られる基本周波数系列を利用して、表情豊かな歌声合成音を実現した。しかしながら、これらの従来技術では、制御パラメータ（減衰率ζと固有周波数Ω）が手作業あるいは規則に基づいて決定されるものであった。これに対し、非特許文献４では、入力となる階段状信号および線形２次系の制御パラメータがいずれも未知の下で、観測される基本周波数系列だけから、それらを同時に推定する確率的なフレームワークが提案された。これは、音高目標値系列を表現する隠れマルコフモデル（ＨＭＭ）と、差分近似に基づく式（１）のパラメトリックな表現によって、最尤なモデルパラメータを反復推定する学習アルゴリズムである。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】N. Minematsu, B. Matsuoka, and K. Hirose, “Prosodic Modeling of Nagauta Singing and Its Evaluation," in Proc. SpeechProsody 2004, pp. 487-490, Mar. 2004
【非特許文献２】T. Saitou, M. Unoki, and M. Akagi, “Development of an F0 control Model Based on F0 Dynamic Characteristics for Singing-Voice Synthesis," Speech Communication, vol.46, pp. 405-417, 2005
【非特許文献３】T. Saitou, M. Goto, M. Unoki, and M. Akagi, “Speech-To-Singing Synthesis: Converting Speaking Voices to Singing Voices by Controlling Acoustic Features Unique to Singing Voices," in Proc. WASPAA 2007, pp. 215-218, Oct. 2007
【非特許文献４】Y. Ohishi, H. Kameoka, K. Kashino, and K. Takeda, “Parameter Estimation Method of F0 Control Model for Singing Voices," in Proc INTERSPEECH 2008, pp. 139-142, Sep. 2008
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上述した非特許文献４の従来技術では、モデルパラメータの推定性能が悪かった。すなわち、観測される基本周波数系列と、パラメータによって再合成される基本周波数系列との誤差が大きくなるという問題がある。これは、基本周波数系列には音高の立ち上がりやビブラートのような様々な動特性が混在しているため、２次系のパラメータを最尤法によって直接推定すると、ある特定の動特性に引っ張られるオーバーフィッティングの問題が発生するためである。これに対し、基本周波数系列をフレーム分割して、フレームごとにモデルパラメータを推定する手法も提案されたが、各動特性を生み出す２次系の影響範囲が系列上で不明確であるため、結局モデルパラメータを適切に推定できないという問題がある。
【０００９】
本発明は、このような事情に鑑みてなされたもので、時系列信号の動特性特徴を抽出する際に、入力となる階段状信号および線形２次系の制御パラメータがいずれも未知の下で、観測される基本周波数系列だけから、モデルパラメータを精度よく推定することができる信号解析装置、信号解析方法及び信号解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
本発明は、観測信号［ｏ］を、入力信号［ｆ］とフィルタ特性を表すインパルス応答信号［ｈ］との畳みこみによって得られる信号生成系の出力信号［ｙ］と残差信号［ε］との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータ［ｕ］と、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータ［｛ａ，ｂ｝または｛ａ０，ａ１，ａ２｝または｛ｗ１，ｗ２，…，ｗＩ｝］と、前記残差信号を表すモデルを構成する残差信号パラメータ［β］とを推定することにより時系列信号の動特性特徴を抽出する信号解析装置であって、前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成部と、前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータ［Θ］とし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離部と、前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの、対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られるＱ関数［式２７］を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新部と、前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離部と前記モデルパラメータ更新部による処理とを再度行わせるパラメータ収束判定部と、前記パラメータ収束判定部により前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力部とを備えることを特徴とする。
【００１１】
本発明は、  前記入力信号はステップ信号［式１７］であり、前記出力信号は多次元ガウス分布に従うものとして確率的にモデル化［式２１］され、前記残差信号は、ガウス性白色雑音として確率的にモデル化［式２２］されることを特徴とする。
【００１２】
本発明は、前記信号生成系のフィルタ特性は、差分法によって導出されるフィルタ［式５］で表され、前記フィルタ特性パラメータは、前記固有周波数の二乗に反比例するパラメータ［ａ＝１／Ω２，  式５］と、前記減衰率に比例し前記固有周波数に反比例するパラメータ［ｂ＝２ζ／Ω，  式５］とであることを特徴とする。
【００１３】
本発明は、前記信号生成系のフィルタ特性は、自己回帰過程に基づいて構成されるフィルタ［式１０］で表され、前記フィルタ特性パラメータは、自己回帰パラメータ［式９］であることを特徴とする。
【００１４】
本発明は  前記信号生成系のフィルタ特性は、複数の２次系フィルタの重み付き線形和によって構成されるフィルタ［式１５］で表され、前記フィルタ特性パラメータは、前記各２次系フィルタの重み［｛ｗ１，ｗ２，…，ｗＩ｝］であることを特徴とする。
【００１５】
本発明は、前記モデルパラメータ更新部は、補助変数から構成される前記目的関数の補助関数を、前記固有周波数の二乗に反比例するパラメータ［ａ］と前記補助関数を前記減衰率に比例し前記固有周波数に反比例するパラメータ［ｂ］とでそれぞれ微分して得られる方程式からなる連立方程式［式３６と式３７］をフィルタ特性パラメータについて解くことにより、フィルタ特性パラメータ［｛ａ，ｂ｝］の値を更新するフィルタ特性パラメータ更新部と、前記補助関数を、前記入力信号パラメータで微分して得られる方程式［式３８］を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、前記補助関数を、前記残差信号パラメータで微分して得られる方程式［式３８］を解くことにより、残差信号パラメータを更新する残差信号パラメータ更新部とから構成されることを特徴とする。
【００１６】
本発明は、前記モデルパラメータ更新部は、前記目的関数を、前記フィルタ特性パラメータに含まれる各自己回帰パラメータ［｛ａ０，ａ１，ａ２｝］でそれぞれ微分して得られる方程式［式４０，  式４１，  式４２］からなる連立方程式をフィルタ特性パラメータについて解くことにより、フィルタ特性パラメータの値を更新するフィルタ特性パラメータ更新部と、前記目的関数を、前記入力信号パラメータで微分して得られる方程式［式４３］を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、前記目的関数を、前記残差信号パラメータで微分して得られる方程式［式４３］を解くことにより、残差信号パラメータを更新する残差信号パラメータ更新部とから構成されることを特徴とする。
【００１７】
本発明は、前記モデルパラメータ更新部は、補助変数から構成される前記目的関数の補助関数を、前記フィルタ特性パラメータである各２次系フィルタの重み［｛ｗ１，ｗ２，…，ｗＩ｝］でそれぞれ微分して得られる方程式からなる非線形連立方程式［式４９または式５７］を、フィルタ特性パラメータについて解くことにより、フィルタ特性パラメータの値を更新するフィルタ特性パラメータ更新部と、前記補助関数を、前記入力信号パラメータで微分して得られる方程式［式５０または式５８］を解くことにより、入力信号パラメータを更新する入力信号パラメータ更新部と、前記補助関数を、前記残差信号パラメータで微分して得られる方程式［式５０または式５８］を解くことにより、残差信号パラメータを更新する残差信号パラメータ更新部とから構成されることを特徴とする。
【００１８】
本発明は、前記信号分離部は、前記観測信号と前記モデルパラメータが与えられた時の、前記出力信号および前記残差信号から構成される完全データの期待値［式（２８）］と、前記完全データの自己相関［式（２９）］とを用いて、前記観測信号を出力信号と残差信号とに分離することを特徴とする。
【００１９】
本発明は、観測信号を、入力信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる信号生成系の出力信号と残差信号との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータと、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータと、前記残差信号を表すモデルを構成する残差信号パラメータとを推定することにより時系列信号の動特性特徴を抽出する信号解析装置における信号解析方法であって、前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成ステップと、前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータとし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離ステップと、前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られる関数を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新ステップと、前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離ステップと前記モデルパラメータ更新ステップによる処理を再度行わせるパラメータ収束判定ステップと、前記パラメータ収束判定ステップにより前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力ステップとを有することを特徴とする。
【００２０】
本発明は、観測信号を、入力信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる信号生成系の出力信号と残差信号との和で表し、前記観測信号から、前記入力信号を表すモデルを構成する入力信号パラメータと、前記フィルタ特性を表すモデルを構成するフィルタ特性パラメータと、前記残差信号を表すモデルを構成する残差信号パラメータとを推定することにより時系列信号の動特性特徴を抽出する信号解析装置上のコンピュータに信号解析を行わせる信号解析プログラムであって、前記観測信号から前記入力信号パラメータと前記フィルタ特性パラメータと前記残差信号パラメータの初期値を生成するパラメータ初期値生成ステップと、前記フィルタ特性パラメータと前記入力信号パラメータと前記残差信号パラメータとの組をモデルパラメータとし、前記モデルパラメータを用いて、前記観測信号を、前記入力信号パラメータと前記フィルタ特性パラメータによって構成される前記信号生成系の出力信号と前記残差信号パラメータによって構成される前記残差信号とに分離する信号分離ステップと、前記観測信号と、前記モデルパラメータと、前記出力信号および前記残差信号の組とが与えられたときの対数尤度関数の条件付き期待値に前記モデルパラメータの事前確率を足し合わせて得られる関数を目的関数として、当該目的関数を前記モデルパラメータに関して最大化するように前記モデルパラメータを更新するモデルパラメータ更新ステップと、前記モデルパラメータが所定の基準を満たしているか否かを判定し、所定の基準を満たしていないと判定された場合に所定の基準を満たすまで、前記信号分離ステップと前記モデルパラメータ更新ステップによる処理を再度行わせるパラメータ収束判定ステップと、前記パラメータ収束判定ステップにより前記モデルパラメータが所定の規準を満たすと判定された場合に、当該モデルパラメータを出力するパラメータ出力ステップとを前記コンピュータに行わせることを特徴とする。
【発明の効果】
【００２１】
本発明によれば、入力となる階段状信号および線形２次系の制御パラメータがいずれも未知の下で、観測される基本周波数系列のみからモデルパラメータを精度よく推定することができるという効果が得られる。
【図面の簡単な説明】
【００２２】
【図１】線形２次系モデルの概念を示す説明図である。
【図２】信号生成系の概念を示す説明図である。
【図３】本発明の第１の実施形態における歌声の基本周波数系列に対する信号解析装置の構成を示すブロック図である。
【図４】歌声の基本周波数系列に対する信号分解モデルを示す説明図である。
【図５】本発明の第５の実施形態における音声信号のメル周波数ケプストラム係数（ＭＦＣＣ）系列に対する信号解析装置の構成を示すブロック図である。
【図６】第４の実施形態におけるモデルパラメータによって生成される信号μと観測基本周波数系列ｏとの比較結果を示す図である。
【図７】第４の実施形態によって推定されたζとΩの歌唱者ごとの平均値を示す図である。
【図８】第４の実施形態によって推定されたｕの歌唱者ごとの頻度分布を示す図である。
【発明を実施するための形態】
【００２３】
以下、図面を参照して、本発明の一実施形態による信号解析装置を説明する。まず、本発明の信号解析装置において、入力となる階段状信号および線形２次系の制御パラメータがいずれも未知の下で、観測される基本周波数系列のみからモデルパラメータを精度よく推定する原理について説明する。
【００２４】
図１は、線形２次系モデルの概念を示す説明図である。図１に示すように、観測信号ｏ（ｔ）を複数の区分に分割し、それぞれの区分における観測信号を、ステップ信号（入力信号）ｆ（ｔ）、フィルタ特性（線形２次系のインパルス応答）ｈ（ｔ）、残差信号（剰余信号）ε（ｔ）の３つの要素に分解するものとする。
【００２５】
ここで、図２を参照して、入力信号ｆ（ｔ）と観測信号ｏ（ｔ）との関係を説明する。図２は、信号生成系の概念を示す説明図である。入力信号ｆ（ｔ）は、目標値からなる信号を意味する。線形の２次系フィルタｈ（ｔ）は、２次系に従う、観測信号の時間的な立ち上がりやオーバーシュートなどの動特性を表現する。線形２次系フィルタの出力信号ｙ（ｔ）は、ｆ（ｔ）とｈ（ｔ）の畳み込み
【数２】
0005318042.tif 000003
によって計算される系の出力信号を表す。残差信号ε（ｔ）は観測信号と出力信号の残差信号を表す。ただし、観測信号は離散的な信号系列であるため、これらの要素に分解するために、式（１）の２次系の伝達関数を離散時間表現する４つの手法のいずれかを用いる。
【００２６】
ここで、線形２次系の伝達関数を離散時間表現する４つの手法について説明する。まず、差分法による離散時間表現法（手法１という）について説明する。式（１）をラプラス逆変換すると、
【数３】
0005318042.tif 000004
となり、ｆ（ｔ）とｙ（ｔ）の関係を表す２階微分方程式が導かれる。ここでｆ（ｔ）とｙ（ｔ）をサンプリング周期Δで離散化して、
【数４】
0005318042.tif 000005
と表し（Ｎは信号の長さを表す）、差分法を利用すると、
【数５】
0005318042.tif 000006
のように各時刻の微分値を近似できる。
【００２７】
この近似を利用すると式（３）は、
【数６】
0005318042.tif 000007
と書ける。ここで、（ａＡ＋ｂＢ＋Ｃ）が２次系の逆フィルタのインパルス応答に相当し、Ａ、Ｂ、ＣはそれぞれＮ×Ｎの行列であり、
【数７】
0005318042.tif 000008
と書ける。
【００２８】
もちろん、差分法は、式（４）以外にも複数の近似方法（中央差分、前進差分、後退差分、Ｓｉｎｃ関数利用など）を用いることができ、結果的にＡ、Ｂ、Ｃそれぞれの行列の構成を変更することに相当する。手法１では、式（５）のａ、ｂがフィルタ特性パラメータとなる。
【００２９】
次に、自己回帰モデルに基づく離散時間表現法（手法２という）について説明する。式（１）のｓ領域からｚ領域への写像を考える。
【数８】
0005318042.tif 000009
と近似すると、ｓ領域とｚ領域の関係は、
【数９】
0005318042.tif 000010
となる。式（１）の逆フィルタの伝達関数はｚ領域で、
【数１０】
0005318042.tif 000011
と書ける。
【００３０】
これは２次の自己回帰モデルと同形であり、
【数１１】
0005318042.tif 000012
とすると、式（８）は、
【数１２】
0005318042.tif 000013
と書ける。Ｕ０＝ＩＮ（ＩＮはＮ×Ｎの単位行列）であり、Ｕ１、Ｕ２はそれぞれ
【数１３】
0005318042.tif 000014
からなるＮ×Ｎの行列とする。手法２では、式（１０）のａ０，ａ１，ａ２がフィルタ特性パラメータとなる。
【００３１】
次に、複数の振動基底の線形和に基づく離散時間表現法（手法３、手法４という）について説明する。式（１）の伝達関数のラプラス逆変換によって得られるインパルス応答は、ζの値によって、以下のように場合分けされる。
【数１４】
0005318042.tif 000015
【００３２】
これらのインパルス応答をサンプリング周期Δに基づいて離散化すると、系の入出力関係は、ｙ＝Φｆのような形式で記述できる。例えば、ζ＝１の場合、Φは下三角行列
【数１５】
0005318042.tif 000016
で書ける。ただし、ｈ（ｔ）は式（１２）のように複数の場合からなるので、行列Φを以下のように構成する。
【数１６】
0005318042.tif 000017
【００３３】
ここでは、予め手動でζ、Ωを決定し、Ｉ個の振動現象を表すインパルス応答｛Φ（１），Φ（２），・・・，Φ（Ｉ）｝を計算する。そして、これらの逆行列Υ（ｉ）：＝（Φ（ｔ））−１（逆フィルタのインパルス応答を表す。これらを以後、振動基底と呼ぶ）の線形重み付き和で、Φ−１を近似する。手法３では、式（１４）のｗ：＝｛ｗ１，ｗ２，・・・，ｗＩ｝をフィルタ特性パラメータとし、これらのパラメータを回帰問題の枠組みで推定するものとする。
【００３４】
一方、手法４では、フィルタ特性パラメータｗ：＝｛ｗ１，ｗ２，・・・，ｗＩ｝がスパースなものとする。これは、Φ−１が、ほんの数個の振動基底のみによって表現されることを意味する。後に説明するが、これはｗの事前確率を想定することで容易に実現されうる。
【００３５】
系のインパルス応答ΦをΦ（ｉ）の線形重み付き和で表現してもよいが、Υ（ｉ）の線形重み付き和で逆フィルタのインパルス応答Φ−１を表現した理由は、後に説明するパラメータ学習アルゴリズムの導出の複雑さを解消するためである。それゆえに、系の入出力関係は、
【数１７】
0005318042.tif 000018
と表現される。
【００３６】
以上に説明した４つの手法（式（５）、式（１０）、式（１５））より、式（１）の伝達関数が、すべてΨｙ＝ｆの形式に変換された。Ψはそれぞれ、
【数１８】
0005318042.tif 000019
と表現される。
【００３７】
次に、ガウス過程に基づく線形２次系の統計的モデリングについて説明する。式（５）、式（１０）、式（１５）によって表現される系の入出力関係を、ガウス過程（文献：C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. MIT Press, Cambridge, Mass, USA, 2006.）に基づいて統計的にモデル化する。
【００３８】
ここで、入力ステップ信号のモデル化について説明する。入力信号ｆはステップ信号を想定する。そのために、常に同じ値をもつベクトル
【数１９】
0005318042.tif 000020
を用意する。ここで、スカラー値ｕは相対音高を表すパラメータ、１ＮはＮ個の１の値が並ぶベクトルとする。このベクトルｕを平均とする多次元ガウス分布から生成される確率変数として、入力信号ｆを確率的にモデル化する。
【数２０】
0005318042.tif 000021
【００３９】
ここで、αは分散を表す超パラメータであり、あらかじめ手動で値を設定する。したがって、請求項で述べた入力信号パラメータはｕのみである。出力ｙは、ガウス分布に従う変数集合ｆの線形結合（ｙ＝Ψ−１ｆ）であるから、ｙ自身もガウス分布に従う。その平均と共分散は、
【数２１】
0005318042.tif 000022
である。よって、ｙが従う確率分布は、
【数２２】
0005318042.tif 000023
と書ける。
【００４０】
このモデルはガウス過程の１つの例となっていることは特筆すべき点である。ガウス過程の重要な点は、ｙのＮ個の要素の同時分布が、平均と共分散といった２次までの統計量で完全に記述される点である。通常のガウス過程では、平均はたいていの場合、零とすることが多く、共分散はカーネル行列の線形和（Multiple Kernel Learning；文献：F. R. Bach, G. R. G. Lanckriet, and M. I. Jordan, “Multiple kernel learning, conic duality, and the smo algorithm," in Proc. ICML 2004, pp. 6-13, July 2004）によって構成されるのが一般的である。この技術により、複数の観測値（ここでは、ｙの要素）がそれぞれ独立同分布に従うと仮定するのではなく、観測値間の相関を考慮できる。このMultiple Kernel Learningを利用したガウス過程による信号モデリングは機械学習の分野で近年注目を集めている。
【００４１】
これに対し、前述したモデルでは、式（２１）の平均と共分散の項に行列Ψ−１を含む。そして、式（１６）のように、Ψがいくつかの基底の線形和で表現されている。したがって、Multiple Kernel Learningとは異なる、特別なガウス過程に基づく信号モデルとここで言及する。
【００４２】
次に、尤度関数と事前確率について説明する。ガウス性白色雑音に従う残差信号
【数２３】
0005318042.tif 000024
を導入し、観測信号
【数２４】
0005318042.tif 000025
は、系の出力信号ｙに残差信号εが加わった信号
【数２５】
0005318042.tif 000026
と仮定する。
【００４３】
ここで、βは残差信号の分散を表す超パラメータとし、これが請求項でいう残差信号パラメータに相当する。ｙとεは互いに独立であると仮定すると、ガウス過程の定義より、観測信号ｏが与えられたときのモデルパラメータθ＝｛Ψ，ｕ，β｝の対数尤度関数は、
【数２６】
0005318042.tif 000027
となる。ここで、μ：＝Ψ−１ｕ，Σ：＝αΨ−１（Ψ−１）Ｔ＋βＩＮとする。
【００４４】
Θの事前確率Ｐ（Θ）は、各要素の独立性Ｐ（Θ）＝Ｐ（Ψ）Ｐ（ｕ）Ｐ（β）を仮定し、ｕとβはそれぞれ一様分布に従うものとする。Ｐ（Ψ）は、前述した手法１〜４のフィルタ特性パラメータの独立性を仮定して、
【数２７】
0005318042.tif 000028
とし、Ｐ（ａ）、Ｐ（ｂ）、Ｐ（ａ０）、Ｐ（ａ１）、Ｐ（ａ２）はそれぞれ一様分布を仮定する。手法３では、同様にＰ（ｗ１），Ｐ（ｗ２），・・・，Ｐ（ｗＩ）はそれぞれ一様分布を仮定する。一方、手法４では、パラメータｗにスパースな制約をもたせるために、一般化正規分布
【数２８】
0005318042.tif 000029
に従うものとする。ただし、ｐ、λは一般化正規分布の形状を規定する定数であり、０＜ｐ＜２のときｐ（ｗ）は優ガウス的となり、スパースネスを測るための尺度となる。
【００４５】
次に、ＥＭ法に基づくパラメータ学習アルゴリズムについて説明する。基本周波数系列ｏが与えられたときに、事後確率Ｐ（Θ｜ｏ）∝Ｐ（ｏ｜Θ）Ｐ（Θ）を最大化するパラメータΘの推定値を決定したい。しかしながら、Θの事後（ＭＡＰ）推定値に関する最適解を解析的に求めることは難しい。その理由は、観測信号ｏが出力信号ｙと残差信号εの和で構成される、尤度関数が、Ψを構成するフィルタ特性パラメータに関して非線形となるなどが挙げられる。ここでは、それぞれの問題に対処するために、次の２つの方策を適用する。
【００４６】
１つ目は、ＥＭ法のＥ−ｓｔｅｐで、ｏをｙとεに分配することである。２つ目は、Ｍ−ｓｔｅｐに補助関数法（文献：H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A New Sparse Representation for Acoustic Signals," in Proc. ICASSP 2009, pp. 3437-3440, April 2009）を適用して、Ｑ関数の補助関数を設計することである。
【００４７】
次に、完全データの定義について説明する。このＭＡＰ推定問題にＥＭ法を適用する際の最初のステップは完全データを定義することである。ここでは、ｙとεを完全データｘと見なして、ＥＭ法を適用する。不完全データと完全データの関係を、
【数２９】
0005318042.tif 000030
と表現する。ｘと現在のパラメータΘ´が与えられたときの、対数尤度関数の条件付き期待値を計算し、さらにｌｏｇＰ（Θ）を加算すると、以下のようなＱ関数を得る。
【数３０】
0005318042.tif 000031
【００４８】
なお、＝の上にｃのついた記号は、定数項を除いて左辺と右辺が等しいことを表し、
【数３１】
0005318042.tif 000032
である。ここで、ｔｒ（・）は行列のトレースを表し、条件付ガウス分布の性質から、
【数３２】
0005318042.tif 000033
と書ける。ここで、式（２８）は、不完全データｏ（観測信号）と現在のパラメータが与えられた時の、完全データｘ（出力信号と残差信号で構成される）の期待値であり、式（２９）は、不完全データｏ（観測信号）と現在のパラメータが与えられた時の、完全データｘ（出力信号と残差信号で構成される）の自己相関である。
【００４９】
ＥＭ法のＥ−ｓｔｅｐでは、直前に更新されたモデルパラメータをΘ´に代入し、Ｅ［ｘ｜ｏ；Θ´］とＥ［ｘｘＴ｜ｏ；Θ´］を計算する。後の計算のため、ｙ、εに対応するようにＥ［ｘ｜ｏ；Θ´］とＥ［ｘｘＴ｜ｏ；Θ´］を
【数３３】
0005318042.tif 000034
のように区分表現する。￣ｘｙと￣ｘεは（￣はそれぞれｘの頭に付く、以下の説明においても同様）はどちらも長さＮのベクトルであり、ＲｙとＲεはどちらもＮ×Ｎの正方行列を表す。
【００５０】
次に、Ｑ関数を前述の手法１〜手法４の各モデルパラメータに関して、最大化するＭ−ｓｔｅｐ更新式について説明する。まず、手法１のＭ−ｓｔｅｐ更新式について説明する。Ψ＝ａＡ＋ｂＢ＋Ｃなので、式（２７）から関連する項を取り出すと、最大化する目的関数は
【数３４】
0005318042.tif 000035
と改めて定義される。ここで、Ａｎ，ｎとは行列Ａのｎ行ｎ列目の要素を表す。パラメータ集合Θに関する最大化問題を解く更新式を前述の補助関数法により導く。式（３１）で与えられる目的関数の補助関数を以下の不等式を用いて導く。
【数３５】
0005318042.tif 000036
ここで、等号成立は、
【数３６】
0005318042.tif 000037
の場合である。補助変数γ１＝｛γａ，１，γｂ，１，γａ，２，γｂ，２，・・・，γａ，Ｎ，γｂ，Ｎ｝を導入すると、補助関数は、
【数３７】
0005318042.tif 000038
とすれば、ｆ１（ａ，ｂ，ｕ，β）≧ｆ１＋（ａ，ｂ，ｕ，β，γ１）が成り立ち、式（３３）が成立するとき、
【数３８】
0005318042.tif 000039
となるため、ｆ１＋（ａ，ｂ，ｕ，β，γ１）は補助関数の定義を満たす。
【００５１】
そこで、ｆ１＋（ａ，ｂ，ｕ，β，γ１）をａ、ｂそれぞれに関して微分して０とおくと、
【数３９】
0005318042.tif 000040
を得る。式（５）から分かるように、ａとｂは共に正の値となる制約の下で、Coordinate descent法を適用して、式（３６）と式（３７）からなる連立方程式を解き、ａとｂを求める。
【００５２】
一方、ｆ１＋（ａ，ｂ，ｕ，β，γ１）をｕ、βそれぞれに関して微分して０とおくと、
【数４０】
0005318042.tif 000041
を得る。これらの式を用いて、ｕとβを更新する。
【００５３】
以上の手法１に関する、ＥＭ法に基づくパラメータ学習アルゴリズムをまとめると、以下のようになる。
初期化： パラメータΘ＝｛ａ，ｂ，u，β｝に初期値を与える。Ｅ−ｓｔｅｐ： Ｅ［ｘ｜ｏ；Θ´］、Ｅ［ｘｘＴ｜ｏ；Θ´］、γ１＝｛γａ，１，γｂ，１，γａ，２，γｂ，２，・・・，γａ，Ｎ，γｂ，Ｎ｝の更新。Ｍ−ｓｔｅｐ：式（３６）、（３７）、（３８）より、モデルパラメータΘ＝｛ａ，ｂ，ｕ，β｝の更新。収束判定：式（３４）の補助関数の値が収束していなければ、Θ´＝ΘとしてＥ−ｓｔｅｐに戻る。
【００５４】
次に、手法２のＭ−ｓｔｅｐ更新式について説明する。Ψ＝ａ０Ｕ０＋ａ１Ｕ１＋ａ２Ｕ２なので、式（２７）から関連する項を取り出すと、最大化する目的関数は、
【数４１】
0005318042.tif 000042
と改めて定義される。ｆ２（ａ０，ａ１，ａ２，ｕ，β）をａ０，ａ１，ａ２に関してそれぞれ微分して０とおくと，
【数４２】
0005318042.tif 000043
を得る。式（９）に示すように、ａ０とａ２は正の値、ａ１は負の値となる制約の下で、Coordinate descent法を適用して、式（４０）、（４１）、（４２）からなる連立方程式を解き、ａ０，ａ１，ａ２を求める。
【００５５】
一方、ｆ２（ａ０，ａ１，ａ２，ｕ，β）をｕ、βそれぞれに関して微分して０とおくと、
【数４３】
0005318042.tif 000044
を得る。これらの式を用いて、ｕとβを更新する。
【００５６】
以上の手法２に関する、ＥＭ法に基づくパラメータ学習アルゴリズムをまとめると、以下のようになる。初期化：パラメータΘ＝｛ａ０，ａ１，ａ２，ｕ，β｝に初期値を与える。Ｅ−ｓｔｅｐ：Ｅ［ｘ｜ｏ；Θ´］、Ｅ［ｘｘＴ｜ｏ；Θ´］の更新。
Ｍ−ｓｔｅｐ：式（４０）、（４１）、（４２）、（４３）より、モデルパラメータΘ＝｛ａ０，ａ１，ａ２，ｕ，β｝の更新。収束判定：式（３９）の目的関数の値が収束していなければ、Θ´＝ΘとしてＥ−ｓｔｅｐに戻る。
【００５７】
次に、手法３のＭ−ｓｔｅｐ更新式について説明する。Ψ＝ｗ１Υ（１）＋ｗ２Υ（２）＋・・・＋ｗＩΥ（Ｉ）なので、式（２７）から関連する項を取り出すと、最大化する目的関数は、
【数４４】
0005318042.tif 000045
と改めて定義される。ここで、Υｎ，ｎ（ｉ）はΥ（ｉ）のｎ行ｎ列目の対角要素を表す。この最大化問題を解く更新式を補助関数法により導く。式（４４）で与えられる目的関数の補助関数を以下の不等式を用いて導く。
【数４５】
0005318042.tif 000046
【００５８】
ここで、補助変数γ３＝｛γ１，１，・・・，γＩ，Ｎ｝を定義し、補助関数を
【数４６】
0005318042.tif 000047
と定義する。このとき、
【数４７】
0005318042.tif 000048
が成り立ち、等号成立は、
【数４８】
0005318042.tif 000049
のときであるため、式（４６）は補助関数の定義を満たす。
【００５９】
式（４６）をｗｉ´に関して微分して０とおくと，
【数４９】
0005318042.tif 000050
を得る。Coordinate descent法を利用して、ｉ´＝１，２，・・・，Ｉに関する非線形連立方程式を解くと、ｗを求めることができる。一方、ｆ３＋（ｗ，ｕ，β，γ３）をｕ、βそれぞれに関して微分して０とおくと、
【数５０】
0005318042.tif 000051
を得る。これらの式を用いて、ｕとβを更新する。
【００６０】
以上の手法３に関する、ＥＭ法に基づくパラメータ学習アルゴリズムをまとめると、以下のようになる。
初期化：パラメータΘ＝｛ｗ，ｕ，β｝に初期値を与える。Ｅ−ｓｔｅｐ：Ｅ［ｘ｜ｏ；Θ´］、Ｅ［ｘｘＴ｜ｏ；Θ´］、γ３の更新。Ｍ−ｓｔｅｐ：式（４９）、（５０）より、モデルパラメータΘ＝｛ｗ，ｕ，β｝の更新。収束判定：式（４６）の補助関数の値が収束していなければ、Θ´＝ΘとしてＥ−ｓｔｅｐに戻る。
【００６１】
次に、手法４のＭ−ｓｔｅｐ更新式について説明する。Ψ＝ｗ１Υ（１）＋ｗ２Υ（２）＋・・・＋ｗＩΥ（Ｉ）なので、式（２７）から関連する項を取り出すと、最大化する目的関数は、
【数５１】
0005318042.tif 000052
と改めて定義される。ここで、Υｎ，ｎ（ｉ）はΥ（ｉ）のｎ行ｎ列目の対角要素を表す。この最大化問題を解く更新式を補助関数法により導く。式（５１）で与えられる目的関数の補助関数を以下の２つの不等式を用いて導く。
【数５２】
0005318042.tif 000053
【００６２】
ここで、補助変数￣ｗ：＝｛￣ｗ１，￣ｗ２，・・・，￣ｗＩ｝、γ４＝｛γ１，１，・・・，γＩ，Ｎ｝を定義し、補助関数を
【数５３】
0005318042.tif 000054
と定義する。このとき、
【数５４】
0005318042.tif 000055
が成り立ち、等号成立は、
【数５５】
0005318042.tif 000056
のときであるため、式（５４）は補助関数の定義を満たす。
式（５４）をｗｉ´に関して微分して０とおくと、
【数５６】
0005318042.tif 000057
を得る。Coordinate descent法を利用して、ｉ´＝１，２，・・・，Ｉに関する非線形連立方程式を解くと、ｗを求めることができる。
【００６３】
一方、ｆ４＋（ｗ，ｕ，β，γ４）をｕ、βそれぞれに関して微分して０とおくと、
【数５７】
0005318042.tif 000058
を得る。これらの式を用いて、ｕとβを更新する。
【００６４】
以上の手法４に関する、ＥＭ法に基づくパラメータ学習アルゴリズムをまとめると、以下のようになる。
初期化：パラメータΘ＝｛ｗ，ｕ，β｝に初期値を与える。Ｅ−ｓｔｅｐ：Ｅ［ｘ｜ｏ；Θ´］、Ｅ［ｘｘＴ｜ｏ；Θ´］、￣ｗ、γ４の更新。Ｍ−ｓｔｅｐ：式（５７）、（５８）から、モデルパラメータΘ＝｛ｗ，ｕ，β｝の更新。収束判定：式（５４）の補助関数の値が収束していなければ、Θ´＝ΘとしてＥ−ｓｔｅｐに戻る。
【００６５】
＜第１の実施形態＞
次に、図３を参照して、本発明の第１の実施形態による信号解析装置の構成を説明する。図３は同実施形態の構成を示すブロック図である。図３に示すように、信号解析装置は、コンピュータ装置によって構成し、基本周波数抽出部１、セグメント分割部２、パラメータ初期値生成部３、信号分離部４、フィルタ特性パラメータ更新部５、入力信号パラメータ更新部６、残差信号パラメータ更新部７、パラメータ収束判定部８およびパラメータ出力部９を備える。フィルタ特性パラメータ更新部５と、入力信号パラメータ更新部６と、残差信号パラメータ更新部７は、モデルパラメータ更新部１０を構成する。第１の実施形態は、前述した手法１を用いて信号解析を行う構成である。
【００６６】
基本周波数抽出部１は、入力される歌声音響信号から観測基本周波数時系列を抽出する。この処理は、周知技術により実現でき、例えば、文献：A de Cheveign´e and H. Kawahara, “YIN, a fundamental frequency estimator for speech and music," Journal of the Acoustical Society of America, vol.111, no.4, pp. 1917-1930, 2002で提案される基本周波数推定法ＹＩＮ  を利用して、歌声音響信号から５ｍｓごとに基本周波数を推定する。
【００６７】
セグメント分割部２は、推定された基本周波数系列をいくつかのセグメントに分割する。図４に示すように、各セグメントは、ある音高から別の音高へ立ち上がる瞬間を始点および終点とする。セグメントへの分割方法は手作業、もしくｋ−ｍｅａｎ法、またはビタビアルゴリズムなどを利用する。分割されたセグメントｏ＝［ｏ１，ｏ２，・・・，ｏＮ］Ｔ（Ｎはセグメントにおける基本周波数系列の長さを表し、セグメントごとに変化する）
ごとにモデルパラメータを推定する。前処理として、セグメントの先頭の基本周波数値ｏ１を、セグメントのすべての基本周波数値から減算し正規化を行う。
【００６８】
パラメータ初期値生成部３は、モデルパラメータΘ１＝｛ａ，ｂ，ｕ，β｝の初期値を決定する。ａ，ｂは、ζ＝１．０，Ω＝０．１のときに、式（５）より計算される、ａ＝１００，ｂ＝２０を初期値とする。ｕは、観測信号ｏの要素の中央値を初期値とする。βは、β＝１００を初期値とする。これらはすべて実験的に決定する。また、αは、α＝２に固定する。
【００６９】
信号分離部４では、ガウス過程の定義およびＥＭアルゴリズムに基づいて導出される出力信号と残差信号の期待値を信号分離基準として、観測信号を出力信号と残差信号に分離する。ここでは、現在のモデルパラメータΘ１´＝｛ａ，ｂ，ｕ，β｝を利用して、式（２８）、（２９）を計算し、式（３０）に基づいて、Ｅ［ｘ｜ｏ；Θ１´］とＥ［ｘｘＴ｜ｏ；Θ１´］を￣ｘｙ，￣ｘε，Ｒｙ，Ｒεに分割する。また、式（３３）に基づいて、補助変数γ１も計算する。
【００７０】
フィルタ特性パラメータ更新部５は、フィルタ特性パラメータであるａ、ｂの値を更新する。式（５）から分かるように、ａとｂは共に正の値となる制約の下で、Coordinate descent法を適用して、式（３６）と式（３７）からなる連立方程式を解き、ａとｂを求める。まず、初期値としてａ＝０、ｂ＝０と設定する。そして、式（３６）をａに関する方程式と見なして、
【数５８】
0005318042.tif 000059
を計算する。
【００７１】
次に、式（３７）をｂに関する方程式と見なして、
【数５９】
0005318042.tif 000060
を計算する。式（５９）、（６０）による更新を交互に繰り返し、ａとｂの値がそれぞれ変化しなくなるまで更新を続ける。
【００７２】
入力信号パラメータ更新部６は、フィルタ特性パラメータ更新部５で更新されたａ、ｂを利用して、式（１６）のΨを計算し、式（３８）に基づいて、入力信号パラメータｕの値を更新する。
【００７３】
残差信号パラメータ更新部７は、式（３８）に基づいて、残差信号パラメータβの値を更新する。
【００７４】
パラメータ収束判定部８は、信号分離部４で計算された￣ｘｙ，Ｒｙ，Ｒεおよび、モデルパラメータ更新部１０によってそれぞれ更新されたモデルパラメータΘ１＝｛ａ，ｂ，ｕ，β｝を利用して、式（３１）の目的関数の値を計算する。更新前のモデルパラメータを用いて計算した式（３１）の目的関数の値と更新後のモデルパラメータを用いて計算した式（３１）の目的関数の値との誤差が、所定の閾値以下であれば、収束したと判定する。収束していればパラメータ出力部９は、モデルパラメータΘ１＝｛ａ，ｂ，ｕ，β｝を出力した後に、次のセグメントにおけるモデルパラメータ推定に移行するため、セグメント分割部２の処理へ移行する。一方、収束しない場合は、Θ１´＝Θ１として信号分離部４の処理に戻る。
【００７５】
なお、収束したか否かを判定する方法としては、目的関数を用いる方法以外に、モデルパラメータ各々の値を更新前と更新後とで比較しても良いし、予め定めた繰り返し回数に到達したか否かで判定を行っても良い。
【００７６】
＜第２の実施形態＞
次に、本発明の第２の実施形態による信号解析装置の構成を説明する。第２の実施形態は、前述した手法２を用いて信号解析を行う構成である。第２の実施形態における信号解析装置の構成は、図３に示す構成と同様であり、基本周波数抽出部１及びセグメント分割部２の処理動作は、第１の実施形態と同じである。第２の実施形態は、その他の構成の処理動作が異なる。
【００７７】
パラメータ初期値生成部３は、モデルパラメータΘ２＝｛ａ０，ａ１，ａ２，ｕ，β｝の初期値を決定する。ａ０，ａ１，ａ２は、ζ＝１．０，Ω＝０．１のときに、式（９）より計算される、ａ０＝１２１，ａ１＝−２２０，ａ２＝１００を初期値とする。ｕは、観測信号ｏの要素の中央値を初期値とする。βは、β＝１００を初期値とする。これらはすべて実験的に決定する。また、αは、α＝２に固定する。
【００７８】
信号分離部４は、ガウス過程の定義およびＥＭアルゴリズムに基づいて導出される出力信号と残差信号の期待値を信号分離基準として、観測信号を出力信号と残差信号に分離する。現在のモデルパラメータΘ２´＝｛ａ０，ａ１，ａ２，ｕ，β｝を利用して、式（２８）、（２９）を計算し、式（３０）に基づいて、Ｅ［ｘ｜ｏ；Θ２´］とＥ［ｘｘＴ｜ｏ；Θ２´］を￣ｘｙ，￣ｘε，Ｒｙ，Ｒεに分割する。
【００７９】
フィルタ特性パラメータ更新部５は、フィルタ特性パラメータであるａ０，ａ１，ａ２の値を更新する。式（９）から分かるように、ａ０とａ２は正の値、ａ１は負の値となる制約の下で、Coordinate descent法を適用して、式（４０）、（４１）、（４２）からなる連立方程式を解き、ａ０，ａ１，ａ２を求める。まず、初期値としてａ０＝０，ａ１＝０，ａ２＝０と設定する。そして、式（４０）をａ０に関する方程式と見なして、
【数６０】
0005318042.tif 000061
を計算する。
【００８０】
次に、式（４１）をａ１に関する方程式と見なして、
【数６１】
0005318042.tif 000062
を計算する。
【００８１】
次に、式（４２）をａ２に関する方程式と見なして、
【数６２】
0005318042.tif 000063
を計算する。そして、式（６１）、（６２）、（６３）による更新を順番に繰り返し、ａ０，ａ１，ａ２の値がそれぞれ変化しなくなるまで更新を続ける。
【００８２】
入力信号パラメータ更新部６は、フィルタ特性パラメータ更新部５で更新されたａ０，ａ１，ａ２を利用して、式（１６）のΨを計算し、式（４３）に基づいて、入力信号パラメータｕの値を更新する。
【００８３】
残差信号パラメータ更新部７は、式（４３）に基づいて、残差信号パラメータβの値を更新する。
【００８４】
パラメータ収束判定部８は、信号分離部４で計算された￣ｘｙ，Ｒｙ，Ｒεおよび、モデルパラメータ更新部１０によって更新されたモデルパラメータΘ２＝｛ａ０，ａ１，ａ２，ｕ，β｝を利用して、式（３９）の目的関数の値を計算する。その値が収束していれば、パラメータ出力部９は、モデルパラメータΘ２＝｛ａ０，ａ１，ａ２，ｕ，β｝を出力し、次のセグメントにおけるモデルパラメータ推定に移行するため、セグメント分割部２の処理動作へ移行する。一方、収束しない場合は、Θ２´＝Θ２として信号分離部４の処理動作に移行する。
【００８５】
＜第３の実施形態＞
次に、本発明の第３の実施形態による信号解析装置の構成を説明する。第３の実施形態は、前述した手法３を用いて信号解析を行う構成である。第３の実施形態における信号解析装置の構成は、図３に示す構成と同様であり、基本周波数抽出部１及びセグメント分割部２の処理動作は、第１の実施形態と同じである。第３の実施形態は、その他の構成の処理動作が異なる。
【００８６】
パラメータ初期値生成部３は、モデルパラメータΘ３＝｛ｗ，ｕ，β｝の初期値を決定する。まず、｛Υ（１），Υ（２），・・・，Υ（Ｉ）｝を作成するために、ζは０から２までの間を０．０２刻みで、Ωは０．０５から０．３までの間を０．００５刻みで変化させる。その結果、Ｉ＝３１００となる。ｗ＝｛ｗ１，ｗ２，・・・ｗＩ｝の初期値はすべて１／Ｉに設定する。ｕは、観測信号ｏの要素の中央値を初期値とする。βは、β＝１００を初期値とする。これらはすべて実験的に決定する。また、αは、α＝２に固定する。
【００８７】
信号分離部４は、ガウス過程の定義およびＥＭアルゴリズムに基づいて導出される出力信号と残差信号の期待値を信号分離基準として、観測信号を出力信号と残差信号に分離する。現在のモデルパラメータΘ３´＝｛ｗ，ｕ，β｝を利用して、式（２８）、（２９）を計算し、式（３０）に基づいて、Ｅ［ｘ｜ｏ；Θ３´］とＥ［ｘｘＴ｜ｏ；Θ３´］を￣ｘｙ，￣ｘε，Ｒｙ，Ｒεに分割する。また、式（４８）に基づいて、補助変数γ３も計算する。
【００８８】
フィルタ特性パラメータ更新部５は、フィルタ特性パラメータであるｗの値を更新する。Coordinate descent法を利用して、式（４９）のｉ´＝１，２，・・・，Ｉに関する非線形連立方程式を解くと、ｗを求めることができる。まず、初期値として｛ｗ１，ｗ２，・・・ｗＩ｝をすべて０に設定する。そして、式（４９）をｗ１´に関する方程式と見なして、
【数６３】
0005318042.tif 000064
と変形する。ｉ´＝１，２，・・・，Ｉに関して、式（６４）による更新を順番に繰り返し、ｗ＝｛ｗ１，ｗ２，・・・ｗＩ｝の値がそれぞれ変化しなくなるまで更新を続ける。
【００８９】
入力信号パラメータ更新部６は、フィルタ特性パラメータ更新部５で更新されたｗを利用して、式（１６）のΨを計算し、式（５０）に基づいて、入力信号パラメータｕの値を更新する。
【００９０】
残差信号パラメータ更新部７は、式（５０）に基づいて、残差信号パラメータβの値を更新する。
【００９１】
パラメータ収束判定部８は、信号分離部４で計算された￣ｘｙ，Ｒｙ，Ｒεおよび、モデルパラメータ更新部１０によって更新されたモデルパラメータΘ３＝｛ｗ，ｕ，β｝を利用して、式（４４）の目的関数の値を計算する。その値が収束していれば、パラメータ出力部９は、モデルパラメータΘ３＝｛ｗ，ｕ，β｝を出力し、次のセグメントにおけるモデルパラメータ推定に移行するため、セグメント分割部２の処理動作へ移行する。一方、収束しない場合は、Θ３´＝Θ３として信号分離部４の処理動作に移行する。
【００９２】
＜第４の実施形態＞
次に、本発明の第４の実施形態による信号解析装置の構成を説明する。第４の実施形態は、前述した手法４を用いて信号解析を行う構成である。第４の実施形態における信号解析装置の構成は、図３に示す構成と同様であり、基本周波数抽出部１及びセグメント分割部２の処理動作は、第１の実施形態と同じである。第４の実施形態は、その他の構成の処理動作が異なる。
【００９３】
パラメータ初期値生成部３は、モデルパラメータをΘ４＝｛ｗ，ｕ，β｝とする以外は第３の実施形態と同じである。
【００９４】
信号分離部４は、ガウス過程の定義およびＥＭアルゴリズムに基づいて導出される出力信号と残差信号の期待値を信号分離基準として、観測信号を出力信号と残差信号に分離する。現在のモデルパラメータΘ４´＝｛ｗ，ｕ，β｝を利用して、式（２８）、（２９）を計算し、式（３０）に基づいて、Ｅ［ｘ｜ｏ；Θ４´］とＥ［ｘｘＴ｜ｏ；Θ４´］を￣ｘｙ，￣ｘε，Ｒｙ，Ｒεに分割する。また、式（５６）に基づいて、補助変数γ４も計算する。
【００９５】
フィルタ特性パラメータ更新部５は、フィルタ特性パラメータであるｗの値を更新する。Coordinate descent法を利用して、式（５７）のｉ´＝１，２，・・・，Ｉに関する非線形連立方程式を解くと、ｗを求めることができる。まず、初期値として｛ｗ１，ｗ２，・・・ｗＩ｝をすべて０に設定する。そして、式（５７）をｗ１´に関する方程式と見なして、
【数６４】
0005318042.tif 000065
と変形する。ｉ´＝１，２，・・・，Ｉに関して、式（６５）による更新を順番に繰り返し、ｗ＝｛ｗ１，ｗ２，・・・ｗＩ｝の値がそれぞれ変化しなくなるまで更新を続ける。
【００９６】
入力信号パラメータ更新部６は、フィルタ特性パラメータ更新部５で更新されたｗを利用して、式（１６）のΨを計算し、式（５８）に基づいて、入力信号パラメータｕの値を更新する。
【００９７】
残差信号パラメータ更新部７は、式（５８）に基づいて、残差信号パラメータβの値を更新する。
【００９８】
パラメータ収束判定部８は、前記信号分離部で計算され￣ｘｙ，Ｒｙ，Ｒεおよび、モデルパラメータ更新部１０によって更新されたモデルパラメータΘ４＝｛ｗ，ｕ，β｝を利用して、式（５１）の目的関数の値を計算する。その値が収束していれば、パラメータ出力部９は、モデルパラメータΘ４＝｛ｗ，ｕ，β｝を出力し、次のセグメントにおけるモデルパラメータ推定に移行するため、セグメント分割部２の処理動作へ移行する。一方、収束しない場合は、Θ４´＝Θ４として信号分離部４の処理動作に移行する。
【００９９】
＜第５の実施形態＞
次に、図５を参照して、本発明の第５の実施形態による信号解析装置の構成を説明する。第４の実施形態は、前述した手法１〜４を使用して、通常の音声信号（話声、歌声を含む）の信号解析を行う構成である。第５の実施形態における信号解析装置の構成は、図３に示す基本周波数抽出部１に代えてメルケプストラム係数抽出部１１を備えた点が、第１〜４の実施形態と異なる。第５の実施形態における信号解析装置は、音声信号から抽出されるメル周波数ケプストラム係数（ＭＦＣＣ）の時系列を、歌声の基本周波数系列と同様に、階段状の入力信号にフィルタが畳み込まれたものと見なし、観測ＭＦＣＣ信号から、音素列を表す入力信号とフィルタ特性の分離特徴抽出を行う。
【０１００】
メルケプストラム係数抽出部１１は、音声信号を周波数分析し、メルケプストラム係数（通常は１２次元程度のベクトル）の時系列を抽出する。第５の実施形態における信号解析装置は、ＭＦＣＣベクトルの各要素の時系列信号ごとに解析処理を行うものである。
【０１０１】
セグメント分割部２、パラメータ初期値生成部３、信号分離部４、フィルタ特性パラメータ更新部５、入力信号パラメータ更新部６、残差信号パラメータ更新部７、パラメータ収束判定部８およびパラメータ出力部９の処理動作は、前述した第１〜第４の実施形態のいずれかに基づくものとする。
【０１０２】
＜実験結果＞
次に、本発明の効果および作用を示すため、本発明の実施形態による信号解析装置を用いた実験結果を以下に説明する。ここでは、上記第１の実施形態、第２の実施形態、第４の実施形態および非特許文献４の従来法を実装して、基本周波数系列の入力ステップ信号（音高目標値系列）とフィルタ特性を表すインパルス応答信号（歌唱動的変動成分）への分解性能を評価する。
【０１０３】
１つ目の評価実験では、本発明が局所最小化問題を解決できているかを確認する。まず、パラメータｕをランダムに１００個決定し、式（１８）に基づいて、１００個の入力ステップ信号を人工的に作成する。ここでは、Ｎ＝３００、α＝２とした。同様に、ζとΩをそれぞれ１００個ランダムに決定し、式（１２）に基づいて、１００個のインパルス応答信号を人工的に作成する。これらの人工的な信号と式（２２）に基づいて、１００個の観測信号を人工的に作成する。ここでは、β＝１００とした。本実験では、これらの観測信号ごとにモデルパラメータを推定する。
【０１０４】
評価尺度として、観測信号ごとに、
（１）推定されたモデルパラメータｕと、観測信号の作成のためにランダムに決定されたｕとの二乗誤差
（２）推定されたモデルパラメータによって構成される系のインパルス応答信号と、観測信号の作成のためにランダムに決定されたζとΩに基づくインパルス応答信号の二乗平均平方根誤差（RootMean Square Error,RMSE）
を計算する。どちらも誤差が小さくなれば、観測信号の入力ステップ信号とインパルス応答信号への分解性能が高い（局所最小化問題を回避できている）と言える。表１はそれぞれの誤差の平均値を示し、最も誤差が小さくなったのは、第４の実施形態よる信号解析装置であった。
【０１０５】
【表１】
0005318042.tif 000066
【０１０６】
２つ目の評価実験では、歌声音響信号から抽出される実際の基本周波数系列を利用してパラメータ学習アルゴリズムの収束性能を評価する。歌声データベースとして、クラシックの声楽家、ポップス歌手、音楽的な訓練を受けていない素人（それぞれ男女１名ずつの計６名）の歌声からなるデータベースを利用した。歌唱者は、伴奏なしで歌唱曲を歌唱した。歌唱曲は、「きらきら星」、「喜びの歌（Ｂｅｅｔｈｏｖｅｎの交響曲第９番第４楽章の歌の部分を岩佐東一郎氏によって作詞されたもの）」である。日本語歌詞による歌唱（２パターンの歌詞）、ハミングによる歌唱を収録した。基本周波数は、前述したＹＩＮを利用して、５ｍｓごとに推定される。なお、Ｈｚで表される周波数ｏＨｚを、次のようにｃｅｎｔで表される対数スケールの周波数ｏｃｅｎｔに変換する。
【数６５】
0005318042.tif 000067
【０１０７】
この変換により、半音は１００ｃｅｎｔに相当する。次に、図４に示すように、推定された基本周波数系列を手動でセグメントに分割する。その結果、セグメントの総数は１３２３となった。セグメントごとにモデルパラメータを推定する。評価尺度として、セグメントごとに、観測信号ｏと、推定されたモデルパラメータによって再合成される信号μ（式（２３）を参照）との二乗平均平方根誤差（ＲＭＳＥ）を計算した。そのＲＭＳＥの平均値を表１の右部分に示す。この評価実験においても第４の実施形態よる信号解析装置が最も小さい誤差となった。
【０１０８】
図６は、観測信号と、第４の実施形態のモデルパラメータによって再合成される信号μを示す。観測信号を見ると、各セグメントには、音高の立ち上がりに関する動特性、音高が安定するときに振動するビブラートのような動特性が複雑に重ね合わさっていることがわかる。本発明は、ＥＭ法を利用して、それぞれを分離する手段を持つため、従来法に比べてパラメータの推定性能が向上した。特に、第４の実施形態の信号解析装置は、様々な振動基底の線形和で２次系のフィルタ特性を表現し、さらにその重みパラメータにスパースな制約を持たせているため、他の実施形態に比べて、観測信号へのオーバーフィッティングの問題を解消し、誤差が最も小さくなっている。
【０１０９】
図７は、第４の実施形態よる信号解析装置によって推定されたζとΩの歌唱者ごとの平均値を示したものである。各セグメントで推定されるｗの要素の最も大きい値ｗｉに着目し、それに対応するζとΩのすべてのセグメントにわたる平均値を歌唱者ごとに計算した。ζが小さな値であることは、その振動現象がオーバーシュートのような減衰振動であることを意味する。一方で、Ωが小さな値であることは、音高の立ち上がり時間が長いことを意味する。したがって、素人歌唱者は歌唱技術が乏しいため、ζとΩの値が他の歌唱者に比べて、ともに小さい値になった。図８は、第４の実施形態よる信号解析装置によって推定されたｕの歌唱者ごとの頻度分布を示す。声楽家やポップス歌手に関しては、半音（１００ｃｅｎｔ）の整数倍の位置に分布のピークが観測される。一方で、素人歌唱者に関しては、そのピークが不鮮明となっている。これも素人歌唱者は歌唱技術が乏しいため、正しい音階（正確な半音単位）で歌唱することが難しいことを意味する。
【０１１０】
このように、観測信号は、入力ステップ信号とフィルタ特性を表すインパルス応答信号との畳みこみによって得られる系の出力信号と、残差信号との和で構成されるものとし、基本周波数系列の場合は、音が立ち上がる（音高が変化する）時点を始点終点として、それによって分割されるセグメントごとに信号の生成過程とする。また、観測信号には、様々な動特性が混在するため、ガウス過程（ベイジアンアプローチ）の観点から信号の生成過程をモデル化し、観測信号の、出力信号と残差信号への分離過程を適用する。そして、分離された出力信号に基づいて、入力信号パラメータやフィルタ特性パラメータを推定する。一方で、分離された残差信号に基づいて、残差信号パラメータを推定する。この推定の流れを何度も繰り返すことにより、最終的にモデルパラメータを精度よく推定することができる。
【０１１１】
以上説明したように、入力信号をステップ信号と仮定し、また、系のフィルタ特性が線形２次系に従うことを仮定して、観測される時系列信号から、信号生成系の入力信号およびフィルタ特性（インパルス応答信号）、残差信号を推定するようにしたため、歌声の基本周波数（Ｆ０）系列から、歌唱者が歌おうとする音高目標値系列とビブラートやオーバーシュートのような歌唱動的変動成分の特徴抽出を行うことができる。
【０１１２】
なお、図３、図５における各処理部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより信号解析処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータシステム」は、ホームページ提供環境（あるいは表示環境）を備えたＷＷＷシステムも含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（ＲＡＭ）のように、一定時間プログラムを保持しているものも含むものとする。
【０１１３】
また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。また、上記プログラムは、前述した機能の一部を実現するためのものであってもよい。さらに、前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるもの、いわゆる差分ファイル（差分プログラム）であってもよい。
【産業上の利用可能性】
【０１１４】
時系列信号の動特性特徴を抽出することによって、信号解析を行うことが不可欠な用途に適用できる。
【符号の説明】
【０１１５】
１・・・基本周波数抽出部、２・・・セグメント分割部、３・・・パラメータ初期値生成部、４・・・信号分離部、５・・・フィルタ特性パラメータ更新部、６・・・入力信号パラメータ更新部、７・・・残差信号パラメータ更新部、８・・・パラメータ収束判定部、９・・・パラメータ出力部、１０・・・モデルパラメータ更新部、１１・・・メルケプストラム係数抽出部
【図１】
0005318042.tif 000068
【図２】
0005318042.tif 000069
【図３】
0005318042.tif 000070
【図４】
0005318042.tif 000071
【図５】
0005318042.tif 000072
【図６】
0005318042.tif 000073
【図７】
0005318042.tif 000074
【図８】
0005318042.tif 000075