(57)【要約】【課題】観測信号の基本周波数が低い時系列信号のモデルパラメータの低周波成分の推定精度を確保できる擬似周期信号のパラメータを推定する信号処理装置を提供する。【解決手段】波数成分出力部１０は、時系列信号をウェーブレット変換して観測時間周波数成分を出力し、初期パラメータ生成部２０が演算初期値および誤差配分比重定数を設定する。補助変数更新部３０は、観測時間周波数成分を各擬似周期信号モデルの各調波成分に分配する補助変数を更新する。周波数集合更新部４０は、各擬似周期信号モデルの調波成分複素振幅パラメータと観測時間周波数成分の分配成分とが複素平面上で鋭角をなす全ての周波数値の周波数値集合を出力する。擬似周期信号モデルパラメータは、擬似周期信号モデルパラメータ更新部５０によって更新され、収束判断条件を満たすと収束判定部６０に判定されると、擬似周期信号モデルパラメータ出力部７０によって出力される。【選択図】図１
2010039215.tif 000002
【特許請求の範囲】
【請求項１】
入力される時系列信号を擬似周期信号または擬似周期信号の重ね合わせとみなし、前記擬似周期信号の特徴を示すパラメータを推定する信号処理装置であって、
前記入力される時系列信号から、マザーウェーブレットの周波数特性が対数正規分布型の関数で定義されるウェーブレット変換によって変換された観測時間周波数成分を出力する観測時間周波数成分出力部と、
前記観測時間周波数成分出力部によって出力された前記観測時間周波数成分を入力として離散時刻ごとの対数基本周波数パラメータと調波成分複素振幅パラメータとを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する初期パラメータ生成部と、
前記観測時間周波数成分出力部によって出力された前記観測時間周波数成分および、設定された前記擬似周期信号モデルパラメータと前記誤差配分比重定数を入力として、前記観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタと各分配成分に重みづけするために用いられるスペクトル重みとを含む補助変数を更新する補助変数更新部と、
前記観測時間周波数成分と、前記擬似周期信号モデルパラメータと前記補助変数更新部によって出力された前記補助変数を入力として、離散時間ごとに前記擬似周期信号モデルの各調波成分の前記調波成分複素振幅パラメータと前記調波成分複素振幅パラメータに対応する前記観測時間周波数成分の分配成分とが複素平面上で鋭角をなす全ての周波数値を含む周波数値集合を出力する周波数集合更新部と、
前記観測時間周波数成分と前記補助変数と前記擬似周期信号モデルパラメータと前記周波数値集合と前記誤差配分比重定数を入力として擬似周期信号モデルパラメータを設定し出力する擬似周期信号モデルパラメータ更新部と、
前記擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には前記補助変数更新部に補助変数の更新を行わせる収束判定部と、
前記収束判定部において所定の規準を満たすと判定された前記擬似周期信号モデルパラメータを出力する擬似周期信号モデルパラメータ出力部と、を備え、
前記補助変数更新部と前記周波数集合更新部と前記擬似周期信号モデルパラメータ更新部は、
初回の演算処理で参照する擬似周期信号モデルパラメータは、前記初期パラメータ生成部によって設定された擬似周期信号モデルパラメータを参照し、繰り返される演算処理で参照する擬似周期信号モデルパラメータは、前記擬似周期信号モデルパラメータ更新部によって設定された擬似周期信号モデルパラメータを参照する
ことを特徴とする信号処理装置。
【請求項２】
前記収束判定部は、
前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、予め定められる目的関数によって定められる目的関数値の変化率が予め定められた所定値以下となる
ことを特徴とする請求項１に記載の信号処理装置。
【請求項３】
前記収束判定部は、
前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、前記擬似周期信号モデルパラメータが反復して演算されて更新される前記擬似周期信号モデルパラメータの変化率が予め定められた所定値以下となる
ことを特徴とする請求項１に記載の信号処理装置。
【請求項４】
前記収束判定部は、
前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、前記擬似周期信号モデルパラメータの反復演算回数が予め定められた所定の回数を満たす
ことを特徴とする請求項１に記載の信号処理装置。
【請求項５】
入力される時系列信号を擬似周期信号または擬似周期信号の重ね合わせとみなし、前記擬似周期信号の特徴を示すパラメータを推定する信号処理方法であって、
前記入力される時系列信号から、マザーウェーブレットの周波数特性が対数正規分布型の関数で定義されるウェーブレット変換によって変換された観測時間周波数成分を出力する観測時間周波数成分出力過程と、
前記観測時間周波数成分出力過程によって出力された前記観測時間周波数成分を入力として離散時刻ごとの対数基本周波数パラメータと調波成分複素振幅パラメータとを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する初期パラメータ生成過程と、
前記観測時間周波数成分出力過程によって出力された前記観測時間周波数成分および、設定された前記擬似周期信号モデルパラメータと前記誤差配分比重定数を入力として、前記観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタと各分配成分に重みづけするために用いられるスペクトル重みとを含む補助変数を更新する補助変数更新過程と、
前記観測時間周波数成分と、前記擬似周期信号モデルパラメータと前記補助変数更新過程によって出力された前記補助変数を入力として、離散時間ごとに前記擬似周期信号モデルの各調波成分の前記調波成分複素振幅パラメータとそれに対応する前記観測時間周波数成分の分配成分とが複素平面上で鋭角をなす全ての周波数値を含む周波数値集合を出力する周波数集合過程と、
前記観測時間周波数成分と前記補助変数と前記擬似周期信号モデルパラメータと前記周波数値集合と前記誤差配分比重定数を入力として擬似周期信号モデルパラメータを設定し出力する擬似周期信号モデルパラメータ更新過程と、
前記擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には前記補助変数更新過程を行わせる収束判定過程と、
前記収束判定過程において所定の規準を満たすと判定された前記擬似周期信号モデルパラメータを出力する擬似周期信号モデルパラメータ出力過程と、を備え、
前記補助変数更新過程と前記周波数集合更新過程と前記擬似周期信号モデルパラメータ更新過程は、
初回の演算処理で参照する擬似周期信号モデルパラメータは、前記初期パラメータ生成過程によって設定された擬似周期信号モデルパラメータを参照し、繰り返される演算処理で参照する擬似周期信号モデルパラメータは、前記擬似周期信号モデルパラメータ更新過程によって設定された擬似周期信号モデルパラメータを参照する
ことを特徴とする信号処理方法。
【請求項６】
コンピュータを請求項１から請求項４のいずれかに記載の信号処理装置として機能させる信号処理プログラム。
【請求項７】
請求項６に記載の信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、主に音声信号における時系列信号に含まれる擬似周期信号のパラメータを推定する信号処理装置、信号処理方法、信号処理プログラムおよび記録媒体に関する。
【背景技術】
【０００２】
人間のコミュニケーション媒体としての音、例えば音声などの有意な情報の大部分は、擬似周期信号で表現できることが多い（例えば、非特許文献１参照）。そのため、信号処理装置には、実世界音響信号から有意な擬似周期信号を分離抽出したり、強調したりする信号処理技術の確立が重要であると考えられる。こうした動機のもと、これまでに擬似周期信号重畳モデルを仮定し、時間領域モデル化誤差のｌ２ノルムをＧａｂｏｒ変換（Ｇａｕｓｓ（ガウス）窓による短時間Ｆｏｕｒｉｅｒ（フーリエ）変換）領域で等価表現し、フレーム毎に擬似周期信号パラメータを最尤推定する安定かつ効率的なアルゴリズムが開発されている（例えば、非特許文献２参照）。
【非特許文献１】Ｒ．Ｊ．ＭｃＡｕｌａｙ  ｅｔ  ａｌ．，ＩＥＥＥ  Ｔｒａｎｓ．  Ａｃｏｕｓｔ．，Ｓｐｅｅｃｈ，Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓ．，３４（４），ｐ．７４４−７５４，１９８６．
【非特許文献２】亀岡弘和、他２名、「正弦波重畳モデルのパラメータ最適化アルゴリズムの導出」、電子情報通信学会技術研究報告、Ｖｏｌ．１０６、ＥＡ２００６−９７、ｐ．４９−５４、２００６．
【発明の開示】
【発明が解決しようとする課題】
【０００３】
ところで、擬似周期信号とは、瞬時瞬時においては周期性があるが、その振幅や周期が時々刻々穏やかに変化する変調信号のことをいう。このような性質の信号は式（１）に示されるようにモデル化することができる。
【０００４】
【数１】
2010039215.tif 000003
【０００５】
式（１）において、ｋは擬似周期信号のインデックス、ｎは擬似周期信号中の調波成分のインデックスを示し、Ａｋ，ｎ（ｔ）は、ｋ番目の擬似周期信号のｎ番目の調波成分の瞬時複素振幅を示し、ｊは虚数単位（ｊ２＝−１）を示す。また、φｋ（ｔ）は、基本波成分（１番目の調波成分）の瞬時位相を示し、ｎ番目の調波成分の瞬時位相がそのｎ倍の速度で進むような拘束が設けられているのがこのモデルの特徴になる。瞬時位相の時間微分は瞬時周波数に相当するため、調波成分の瞬時周波数は常に基本波成分の瞬時周波数（基本周波数）の整数倍になる。したがって、以上のモデルは全体としては必ずとも周期性はないが、瞬時瞬時においては調波性を有することになる。
観測信号をｙ（ｔ）で示し、これをＫ個の擬似周期信号ｓｋ（ｔ）を重ね合わせたモデルによって最も良く近似させる問題を考える。このモデルをｓ（ｔ）で示し、擬似周期信号ｓｋ（ｔ）との関係を式（２）に示す。
【０００６】
【数２】
2010039215.tif 000004
【０００７】
また、観測信号ｙ（ｔ）とモデルｓ（ｔ）の関係を、式（３）に示す。
【０００８】
【数３】
2010039215.tif 000005
【０００９】
式（３）において、ｎ（ｔ）はＧａｕｓｓ（ガウス）性白色雑音とすると、式（４）の特性を有するものと仮定する。
【００１０】
【数４】
2010039215.tif 000006
【００１１】
式（４）で示されるようにｎ（ｔ）の分布は、正規分布Ν（０，ν２）にしたがうこととすると、モデルパラメータ（以降、正弦波パラメータという）は、式（５）に示される。
【００１２】
【数５】
2010039215.tif 000007
【００１３】
式（５）において、μｋ（ｔ）は、瞬時基本周波数を示し、この式に示される「≡」の符号の意味は、左辺に示される式は、右辺に示される式で定義されることを示す。なお、右辺の文字「φｋ（ｔ）」の上に「・」を付けた文字は、「φｋ（ｔ）」の１回微分を示す。
この式（５）における、最尤解は式（６）を解くことで得ることができる。
【００１４】
【数６】
2010039215.tif 000008
【００１５】
非特許文献２では、モデルｓ（ｔ）の短時間Ｆｏｕｒｉｅｒ（フーリエ）変換（以下、ＳＴＦＴ（Short Time Fourier Transformation）という）モデルを用いて、式（６）に示した目的関数をＧａｂｏｒ（ガボール）変換（ＳＴＦＴ）領域で等価表現し、これを規準として擬似周期信号のＳＴＦＴモデルのパラメータをフレームごとに推定するアルゴリズムが導かれている。
【００１６】
しかしながら、ＳＴＦＴ変換は定帯域フィルタバンク分析に相当するため、非特許文献２で示される技術では、観測信号の各周波数成分に対し、高周波および低周波の正弦波パラメータを同じ長さの短時間窓内でフィッティングしようとしていたことになる。その結果、低周波成分を解析するのに十分な長さを確保することができず、必然的に低周波成分の推定性能が相対的に劣悪なものとならざるを得ないという問題がある。
また、ＳＴＦＴ変換では観測信号の各周波数成分の周波数分解能が低くなることになるという問題がある。この問題は、時間周波数解析における、周波数分解能と時間分解能の積は一定であること（時間周波数計測の不確定性）によるもので、低周波の正弦波パラメータを高周波のそれと同程度の精度で推定するためには、周波数に反比例した時間解析区間が必要になるという問題がある。
【００１７】
本発明は、上記問題を解決すべくなされたもので、その目的は、観測信号の基本周波数が低い時系列信号においても、モデルパラメータの低周波成分の推定精度を確保できる擬似周期信号のパラメータを推定する信号処理装置、信号処理方法、信号処理プログラムおよび記録媒体を提供することにある。
【課題を解決するための手段】
【００１８】
本発明は、入力される時系列信号を擬似周期信号または擬似周期信号の重ね合わせとみなし、前記擬似周期信号の特徴を示すパラメータを推定する信号処理装置であって、前記入力される時系列信号から、マザーウェーブレットの周波数特性が対数正規分布型の関数で定義されるウェーブレット変換によって変換された観測時間周波数成分を出力する観測時間周波数成分出力部と、前記観測時間周波数成分出力部によって出力された前記観測時間周波数成分を入力として離散時刻ごとの対数基本周波数パラメータと調波成分複素振幅パラメータとを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する初期パラメータ生成部と、前記観測時間周波数成分出力部によって出力された前記観測時間周波数成分および、設定された前記擬似周期信号モデルパラメータと前記誤差配分比重定数を入力として、前記観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタと各分配成分に重みづけするために用いられるスペクトル重みとを含む補助変数を更新する補助変数更新部と、前記観測時間周波数成分と、前記擬似周期信号モデルパラメータと前記補助変数更新部によって出力された前記補助変数を入力として、離散時間ごとに前記擬似周期信号モデルの各調波成分の前記調波成分複素振幅パラメータとそれに対応する前記観測時間周波数成分の分配成分が複素平面上で鋭角をなす全ての周波数値を含む周波数値集合を出力する周波数集合更新部と、前記観測時間周波数成分と前記補助変数と前記擬似周期信号モデルパラメータと前記周波数値集合と前記誤差配分比重定数を入力として擬似周期信号モデルパラメータを設定し出力する擬似周期信号モデルパラメータ更新部と、前記擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には前記補助変数更新部に補助変数の更新を行わせる収束判定部と、前記収束判定部において所定の規準を満たすと判定された前記擬似周期信号モデルパラメータを出力する擬似周期信号モデルパラメータ出力部と、を備え、前記補助変数更新部と前記周波数集合更新部と前記擬似周期信号モデルパラメータ更新部は、初回の演算処理で参照する擬似周期信号モデルパラメータは、前記初期パラメータ生成部によって設定された擬似周期信号モデルパラメータを参照し、繰り返される演算処理で参照する擬似周期信号モデルパラメータは、前記擬似周期信号モデルパラメータ更新部によって設定された擬似周期信号モデルパラメータを参照することを特徴とする信号処理装置である。
【００１９】
また、本発明は、上記記載の発明において、前記収束判定部は、前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、予め定められる目的関数によって定められる目的関数値の変化率が予め定められた所定値以下となることを特徴とする。
【００２０】
また、本発明は、上記記載の発明において、前記収束判定部は、前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、前記擬似周期信号モデルパラメータが反復して演算されて更新される前記擬似周期信号モデルパラメータの変化率が予め定められた所定値以下となることを特徴とする。
【００２１】
また、本発明は、上記記載の発明において、前記収束判定部は、前記擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、前記擬似周期信号モデルパラメータの反復演算回数が予め定められた所定の回数を満たすことを特徴とする。
【００２２】
また、本発明は、入力される時系列信号を擬似周期信号または擬似周期信号の重ね合わせとみなし、前記擬似周期信号の特徴を示すパラメータを推定する信号処理方法であって、前記入力される時系列信号から、マザーウェーブレットの周波数特性が対数正規分布型の関数で定義されるウェーブレット変換によって変換された観測時間周波数成分を出力する観測時間周波数成分出力過程と、前記観測時間周波数成分出力過程によって出力された前記観測時間周波数成分を入力として離散時刻ごとの対数基本周波数パラメータと調波成分複素振幅パラメータとを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する初期パラメータ生成過程と、前記観測時間周波数成分出力過程によって出力された前記観測時間周波数成分および、設定された前記擬似周期信号モデルパラメータと前記誤差配分比重定数を入力として、前記観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタと各分配成分に重みづけするために用いられるスペクトル重みとを含む補助変数を更新する補助変数更新過程と、前記観測時間周波数成分と、前記擬似周期信号モデルパラメータと前記補助変数更新過程によって出力された前記補助変数を入力として、離散時間ごとに前記擬似周期信号モデルの各調波成分の前記調波成分複素振幅パラメータとそれに対応する前記観測時間周波数成分の分配成分が複素平面上で鋭角をなす全ての周波数値を含む周波数値集合を出力する周波数集合過程と、前記観測時間周波数成分と前記補助変数と前記擬似周期信号モデルパラメータと前記周波数値集合と前記誤差配分比重定数を入力として擬似周期信号モデルパラメータを設定し出力する擬似周期信号モデルパラメータ更新過程と、前記擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には前記補助変数更新過程を行わせる収束判定過程と、前記収束判定過程において所定の規準を満たすと判定された前記擬似周期信号モデルパラメータを出力する擬似周期信号モデルパラメータ出力過程と、を備え、前記補助変数更新過程と前記周波数集合更新過程と前記擬似周期信号モデルパラメータ更新過程は、初回の演算処理で参照する擬似周期信号モデルパラメータは、前記初期パラメータ生成過程によって設定された擬似周期信号モデルパラメータを参照し、繰り返される演算処理で参照する擬似周期信号モデルパラメータは、前記擬似周期信号モデルパラメータ更新過程によって設定された擬似周期信号モデルパラメータを参照することを特徴とする信号処理方法である。
【００２３】
また、本発明は、コンピュータを請求項１から請求項４のいずれかに記載の信号処理装置として機能させる信号処理プログラムである。
【００２４】
また、本発明は、請求項６に記載の信号処理プログラムを記録したコンピュータ読み取り可能な記録媒体である。
【発明の効果】
【００２５】
本発明は、時系列信号から、ウェーブレット変換によって変換された観測時間周波数成分から離散時刻ごとの擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する。また、観測時間周波数成分を各擬似周期信号モデルの各調波成分分配する分配フィルタと各分配成分のスペクトル重みとを更新して、擬似周期信号モデルの調波成分複素振幅パラメータと対応する観測時間周波数成分が複素平面上で鋭角をなす全ての周波数値を抽出し、擬似周期信号モデルパラメータを設定し、所定の規準を満たすまで演算を繰り返す。
これにより、基本周波数が低い時系列信号においても、モデルパラメータの低周波成分の推定精度を確保できる擬似周期信号のパラメータを推定することができる。
【００２６】
また、本発明は、上記記載の発明において、擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、予め定められる目的関数によって定められる目的関数値の変化率が予め定められた所定値以下となることとする。
これにより、繰り返し行われる擬似周期信号モデルパラメータの演算処理で、擬似周期信号モデルパラメータが推定値近傍に収束して、擬似周期信号モデルパラメータの変化率が予め定められた所定値以下になった段階で推定処理を打ち切ることにより、演算時間を短縮することができ、効率よく擬似周期信号のパラメータを推定することができる。
【００２７】
また、本発明は、上記記載の発明において、擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、擬似周期信号モデルパラメータが反復して演算されて更新される擬似周期信号モデルパラメータの変化率が予め定められた所定値以下となることとする。
これにより、繰り返し行われる擬似周期信号モデルパラメータの演算処理で、擬似周期信号モデルパラメータが推定値近傍に収束して、擬似周期信号モデルパラメータの変化率が予め定められた所定値以下になった段階で推定処理を打ち切ることにより、演算時間を短縮することができ、効率よく擬似周期信号のパラメータを推定することができる。
【００２８】
また、本発明は、上記記載の発明において、擬似周期信号モデルパラメータの演算結果の判定において必要とされる所定の規準を、擬似周期信号モデルパラメータの反復演算回数が予め定められた所定の回数を満たすこととする。
これにより、繰り返し行われる擬似周期信号モデルパラメータの演算処理の反復演算回数を予め定められた所定の回数で打ち切ることで演算時間を所定の時間内とすることができ、効率よく擬似周期信号のパラメータを推定することができる。
【発明を実施するための最良の形態】
【００２９】
以下、本発明の信号処理装置１００について図面を参照して説明する。
信号処理装置１００は、入力される音響信号にしたがって、その入力された音響信号を時系列信号ｙ（ｔ）として処理することにより時系列信号ｙ（ｔ）に含まれる擬似周期信号のパラメータを推定し、出力する。
（第１実施形態）
図１は、第１実施形態による信号処理装置１００を示す概略ブロック図である。
信号処理装置１００は、観測時間周波数成分出力部１０、初期パラメータ生成部２０、補助変数更新部３０、周波数値集合更新部４０、擬似周期信号モデルパラメータ更新部５０、収束判定部６０、擬似周期信号モデルパラメータ出力部７０を備える。
また、信号処理装置１００は、図示されないＣＰＵ（Central Processing Unit）とＲＡＭ（Random Access Memory）とＲＯＭ（Read Only Memory）とＥＥＰＲＯＭ（Electrically Erasable and Programmable Read Only Memory）とＩ／Ｏ（Input/Output）などが備えられており、演算処理を実行するＣＰＵがＲＡＭを一時的なワークエリアや情報記憶領域として使用したりプログラム領域として使用したりしながら、ＲＯＭに書き込まれた基本プログラムを適宜実行し、Ｉ／Ｏを介して接続されている外部機器および内部機器などを制御している。また、信号処理装置１００は、内部の記憶領域として備えるＲＯＭならびにＥＥＰＲＯＭを処理プログラム、各種テーブルとしての情報、処理結果など情報の記憶領域として利用され、各機能部の処理により記録され、記憶された情報が各機能部の処理により参照される。
【００３０】
信号処理装置１００において観測時間周波数成分出力部１０は、時間ｔにおいて入力される音響信号を時系列信号ｙ（ｔ）として処理することにより、その時系列信号ｙ（ｔ）に基づいて周波数ωにおいてウェーブレット変換によって変換された観測時間周波数成分Ｗｙ（ｘ，ｔ）を出力する。対数周波数ｘは、周波数ωをｌｏｇで演算した変数である。
この処理を式（７）に示す。
【００３１】
【数７】
2010039215.tif 000009
【００３２】
式（７）において、Ｆ−１［・］はＦｏｕｒｉｅｒ（フーリエ）逆変換を示す。Ｙ（ω）は、式（８）で定義される。
【００３３】
【数８】
2010039215.tif 000010
【００３４】
式（８）において、Ｆ［・］はＦｏｕｒｉｅｒ（フーリエ）変換を示す。この式に示される「≡」の符号の意味は、左辺に示される式は、右辺に示される式で定義されることを示し、以下同様とする。
また、Ψ＊（ω）は、式（９）に示される対数正規分布型の関数として定義される。
【００３５】
【数９】
2010039215.tif 000011
【００３６】
式（９）において、「Ψ＊（ω）」は「Ψ（ω）」との複素共役を示す。また、σは任意に定められる定数とする。
観測時間周波数成分出力部１０は、時系列信号ｙ（ｔ）を入力信号として、マザーウェーブレットの周波数特性が対数正規分布型の関数Ψ＊（ω）によって定義されたウェーブレット変換によって変換された観測時間周波数成分Ｗｙ（ｘ，ｔ）を出力する。
【００３７】
初期パラメータ生成部２０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）を入力として離散時刻ごとの対数基本周波数パラメータと調波成分複素振幅パラメータとを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数を設定する。
初期パラメータ生成部２０は、Ｋ個の擬似周波信号とＮ個の調波成分によって示される、対数基本周波数パラメータ（Ωｋ（ｔ））１≦ｋ≦Ｋと調波成分複素振幅パラメータ（Ｂｋ，ｎ（ｔ））１≦ｋ≦Ｋ，１≦ｎ≦Ｎの初期値を生成する。これらの初期値は、乱数により生成することとしてもよいが、例えば、対数基本周波数パラメータ（Ωｋ（ｔ））の初期値は、時間ｔごとに、観測時間周波数成分Ｗｙ（ｘ，ｔ）が対数周波数ｘ方向に極大あるいは極小になる点をすべて検出し、その中で｜Ｗｙ（ｘ，ｔ）｜で示される観測時間周波数成分Ｗｙ（ｘ，ｔ）の絶対値の大きいものから順に選択したＫ個の周波数と置くのが好適である。
これを（Ωｋ（ｔ））ｉｎｉｔで示すと、調波成分複素振幅パラメータＢｋ，ｎ（ｔ）の初期値は、Ｗｙ（（Ωｋ（ｔ））ｉｎｉｔ＋ｌｏｇｎ，ｔ）とするのが好適である。
また、誤差配分比重定数λｋ，ｎ（ｔ）は、∀ｔ（以下、「すべてのｔにおいて」を示す），Σｋ，ｎλｋ，ｎ（ｔ）＝１を満たし、かつ非負である任意の定数であるので、例えば、式（１０）と置くことができる。
【００３８】
【数１０】
2010039215.tif 000012
【００３９】
初期パラメータ生成部２０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）を入力として離散時刻ごとの対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）とを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数λｋ，ｎ（ｔ）を設定する。
【００４０】
補助変数更新部３０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）および、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータと誤差配分比重定数を入力として、観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタと各分配成分に重みづけするために用いられるスペクトル重みとを含む補助変数を更新する。
補助変数更新部３０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）および、擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、誤差配分比重定数をλｋ，ｎ（ｔ）を用いて分配フィルタｍｋ，ｎ（ｘ，ｔ）の更新値を式（１１）に示される式で演算し更新する。
【００４１】
【数１１】
2010039215.tif 000013
【００４２】
また、対数基本周波数パラメータΩｋ（ｔ）を入力信号に用いて、スペクトル重みαｋ，ｎ（ｘ，ｔ）の更新値を式（１２）に示される式で演算し更新する。
【００４３】
【数１２】
2010039215.tif 000014
【００４４】
補助変数更新部３０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）および、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、誤差配分比重定数をλｋ，ｎ（ｔ）を入力として、観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタｍｋ，ｎ（ｘ，ｔ）と各分配成分に重みづけするために用いられるスペクトル重みαｋ，ｎ（ｘ，ｔ）とを含む補助変数を更新する。
【００４５】
周波数値集合更新部４０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータと補助変数更新部３０によって出力された補助変数を入力として、離散時間ごとに擬似周期信号モデルの各調波成分の調波成分複素振幅パラメータとそれに対応する観測時間周波数成分の分配成分とが複素平面上で鋭角をなす全ての周波数値を含む周波数値集合を出力する。
周波数値集合更新部４０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、擬似周期信号モデルパラメータの要素である調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と補助変数更新部３によって出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）を用いて周波数値集合
【００４６】
【数１３】
2010039215.tif 000015
【００４７】
（以下、Ｓ（ｔ）ｋ，ｎと示す）を式（１３）により定める。
【００４８】
【数１４】
2010039215.tif 000016
【００４９】
式（１３）が意味するところは、全体集合「Ｒ」に属する任意の対数周波数ｘにおいて、
【００５０】
【数１５】
2010039215.tif 000017
【００５１】
となる、すなわち、調波成分複素振幅パラメータＢｋ，ｎ（ｔ）とそれに対応する分配成分
【００５２】
【数１６】
2010039215.tif 000018
【００５３】
とが、複素平面上で鋭角をなす全ての対数周波数ｘの集合Ｓ（ｔ）ｋ，ｎとする、ということである。
式（１３）における「Ｒｅ［ａ］」は、複素数ａの実部を示す。
周波数値集合更新部４０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータの要素である調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と補助変数更新部３０によって出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）とを入力とする。
また、周波数値集合更新部４０は、その入力によって、離散時間ごとに擬似周期信号モデルの各調波成分の調波成分複素振幅パラメータＢｋ，ｎ（ｔ）とそれに対応する観測時間周波数成分の分配成分ｍｋ，ｎ（ｘ，ｔ）Ｗｙ（ｘ，ｔ）とが複素平面上で鋭角をなす全ての周波数値を含む周波数値集合Ｓ（ｔ）ｋ，ｎを出力する。
【００５４】
擬似周期信号モデルパラメータ更新部５０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータと補助変数、および、周波数値集合と誤差配分比重定数を入力として擬似周期信号モデルパラメータを更新し出力する。
擬似周期信号モデルパラメータ更新部５０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、補助変数更新部３０で出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）とスペクトル重みαｋ，ｎ（ｘ，ｔ）と、周波数値集合Ｓ（ｔ）ｋ，ｎと誤差配分比重定数λｋ，ｎ（ｔ）を入力として擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）を式（１４）により設定する。
【００５５】
【数１７】
2010039215.tif 000019
【００５６】
また、調波成分複素振幅パラメータＢｋ，ｎ（ｔ）を式（１５）により設定する。
【００５７】
【数１８】
2010039215.tif 000020
【００５８】
擬似周期信号モデルパラメータ更新部５０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、初期パラメータ生成部２０または擬似周期信号モデルパラメータ更新部５０によって出力された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、補助変数更新部３０で出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）とスペクトル重みαｋ，ｎ（ｘ，ｔ）と、周波数値集合Ｓ（ｔ）ｋ，ｎと誤差配分比重定数λｋ，ｎ（ｔ）を入力とする。
また、擬似周期信号モデルパラメータ更新部５０は、擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）を出力する。
【００５９】
収束判定部６０は、擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には補助変数更新部３０からの処理を繰り返し行わせる。
収束判定部６０における所定の規準を満たすか否かの判定基準は、予め定められる目的関数によって定められる目的関数値の変化率が予め定められた所定値以下となることとすることを基準とする。
目的関数を定める式を式（１６）に示す。
【００６０】
【数１９】
2010039215.tif 000021
【００６１】
擬似周期信号モデルパラメータ出力部７０は、収束判定部６０において所定の規準を満たすと判定された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）を出力する。
【００６２】
（目的関数）
続いて、ウェーブレット変換領域における等価な目的関数について説明する。
非特許文献２に示された目的関数は、ウェーブレット変換領域における等価な目的関数として等価表現することができる。まず、説明に先立ち、ここで扱うウェーブレット変換の定義について示す。アドミッシブル条件を満たす中心周波数のマザーウェーブレットをψ（ｕ）とし、時刻ｔに局在する周期ｅ−ｘのウェーブレット基底関数ψｘ，ｔ（ｕ）式（１７）に示す。
【００６３】
【数２０】
2010039215.tif 000022
【００６４】
式（１７）においてｘは対数周波数、ｕは時間を示す。式（１７）を用いることで信号ｙ（ｕ）のウェーブレット変換は、式（１８）、式（１９）として定義することができる。
【００６５】
【数２１】
2010039215.tif 000023
【００６６】
【数２２】
2010039215.tif 000024
【００６７】
式（１８）、式（１９）において、「＜ａ，ｂ＞」は、ａとｂの内積を意味する。ただし、Ｙ（ω）は式（２０）によって定義され、ｙ（ｕ）のフーリエ変換によって定められる。
【００６８】
【数２３】
2010039215.tif 000025
【００６９】
また、Ψｘ，ｔ（ω）は式（２１）によって定義され、ψｘ，ｔ（ｕ）のフーリエ変換によって定められる。
【００７０】
【数２４】
2010039215.tif 000026
【００７１】
また、Ψ（ω）は式（２２）によって定義され、ψ（ｕ）のフーリエ変換によって定められる。
【００７２】
【数２５】
2010039215.tif 000027
【００７３】
式（１９）に示される関係は、Ｐａｒｓｅｖａｌ（パーセバル）の定理によるものである。式（１９）と式（２１）から、式（２３）を導くことができる。
【００７４】
【数２６】
2010039215.tif 000028
【００７５】
式（２３）において、Ψ（ω）は任意であるので、式（２４）として示すことができる。
【００７６】
【数２７】
2010039215.tif 000029
【００７７】
マザーウェーブレットの周波数特性を上記のように定義すると、次に示す定理１が成り立つ。
【００７８】
＜定理１＞
Ｗｙ（ｘ，ｔ）、Ｗｙ（ｘ，ｔ）をそれぞれｙ（ｔ）、ｓ（ｔ）のウェーブレット変換とすると、式（２５）が成り立つ。
【００７９】
【数２８】
2010039215.tif 000030
【００８０】
式（２５）において、「｜・｜」で示される記号は、絶対値を示し、以下、同様とする。
証明：式（２５）に示された等式は、次のように示すことができる。
【００８１】
【数２９】
2010039215.tif 000031
【００８２】
ただし、Ｙ（ω）、Ｓ（ω）は、それぞれｙ（ｔ）、ｓ（ｔ）のフーリエ変換である（証明終）。
【００８３】
上記の定理１から言えることは、ｓ（ｔ）のパラメータ最適化が、ウェーブレット変換領域における誤差のｌ２ノルム規準をもとに行うことができるということである。
【００８４】
続いて、擬似周期信号のウェーブレット変換モデルについて説明する。
ここでは、連続ウェーブレット変換ＷＳ(ｘ，ｔ)を具体的に導出する。定義より、ｓ(ｕ)の連続ウェーブレット変換ＷＳ(ｘ，ｔ)は式（２６）に示すことができる。
【００８５】
【数３０】
2010039215.tif 000032
【００８６】
式（２６）において、ウェーブレット変換では一般に、
【００８７】
【数３１】
2010039215.tif 000033
【００８８】
（以下、ψ＊ｘ，ｔ（ｕ)と記載する）は時刻ｔに局在すること、瞬時位相φｋ(ｕ)と調波成分の瞬時振幅Ａｋ，ｎ(ｕ)は緩やかに変化することを前提として扱うことができることから、時刻ｕ＝ｔ周辺で、瞬時位相φｋ(ｕ)を式（２７）で示し、０次で近似することとする。
【００８９】
【数３２】
2010039215.tif 000034
【００９０】
式（２７）において、「〜と＿」の組み合わせによって示される記号は、右辺を近似したものを左辺として定義することを示す。
また、時刻ｕ＝ｔ周辺で、瞬時振幅Ａｋ，ｎ(ｕ)を式（２８）で示し、１次で近似することとする。
【００９１】
【数３３】
2010039215.tif 000035
【００９２】
式（２７）と式（２８）において０次および１次近似を用いても計算結果には大きな影響を与えないことが見込まれる。瞬時位相の導関数は瞬時周波数なので、
【００９３】
【数３４】
2010039215.tif 000036
【００９４】
（以下、「・φ」と示す）は瞬時基本周波数を示し、これをμｋ（ｕ）と置くことにする。なお、「・φ」および「・・φ」は、文字「φ」の上に「・」「・・」を付けた文字を示し、「・φ」および「・・φ」は、それぞれ「φ」の１回微分、２回微分を示すこととする。
式（２７）、式（２８）より式（２６）は、式（２９）に示されるように近似することができる。
【００９５】
【数３５】
2010039215.tif 000037
【００９６】
また、一般化Ｐａｒｓｅｖａｌの等式より、式（３０）に示される関係がある。
【００９７】
【数３６】
2010039215.tif 000038
【００９８】
すなわち、式（２９）は式（３１）に示すように整理することができる。
【００９９】
【数３７】
2010039215.tif 000039
【０１００】
式（２４）を式（３１）に代入すると、式（３２）が得られる。
【０１０１】
【数３８】
2010039215.tif 000040
【０１０２】
ただし、式（３２）において、Ωｋ（ｔ）は式（３３）と定義する。
【０１０３】
【数３９】
2010039215.tif 000041
【０１０４】
また、式（３２）において、Ｂｋ，ｎ（ｔ）は式（３４）と定義する。
【０１０５】
【数４０】
2010039215.tif 000042
【０１０６】
以上より、式（２７）、式（２８）の仮定の下では式（３５）を最小化することと、式（６）を達成することは等価である。
【０１０７】
【数４１】
2010039215.tif 000043
【０１０８】
ここで、Ｗｙ(ｘ，ｔ)は離散時間観測データであることを想定し、式（３６）を最小化するΘ≡｛Ωｋ（ｔ），｛Ｂｋ，ｎ（ｔ）｝１≦ｎ≦Ｎ｝１≦ｋ≦Ｋを各時刻で個別に求めることを考える。
【０１０９】
【数４２】
2010039215.tif 000044
【０１１０】
従って、時刻ｔのインデックスは、上記に示した式においてすべて省かれている。ここで、以上のようにして導出された式（３７）を、擬似周期信号のウェーブレット変換モデルと呼ぶこととする。
【０１１１】
【数４３】
2010039215.tif 000045
【０１１２】
（最適パラメータ推定アルゴリズム）
続いて、最適パラメータ推定アルゴリズムについて説明する。
最適パラメータ推定アルゴリズムは、以下に示す補助関数を用いた目的関数値の降下法（以後、補助関数法）の考え方を基礎とする。まず補助関数の定義を行い、この定義に基づき補助関数法の収束性を保証する定理を示す。
【０１１３】
＜定義１＞
パラメータΘ＝（Θ１，・・・,ΘI）に関して最小化したい目的関数をΦ（Θ）とすると、式（３８）を満たすとき、Φ＋（Θ，ｍ）を目的関数Φ（Θ）の補助関数、ｍ＝（ｍ１，・・・，ｍＪ）を補助変数と定義する。
【０１１４】
【数４４】
2010039215.tif 000046
【０１１５】
＜定理２(補助関数法)＞
目的関数をΦ(Θ)、目的関数Φ(Θ)の補助関数をΦ+(Θ，ｍ)、補助変数をｍ＝（ｍ１，・・・，ｍＪ）とすると、補助関数の、補助変数に関する最小化を行い、式（３９）に示す。
【０１１６】
【数４５】
2010039215.tif 000047
【０１１７】
式（３９）において、左辺に示した「ｍ」の上に「＾」を付けた文字で示される変数（以下、「＾ｍ」と記載する）は、補助変数Θを固定にして補助変数ｍに関して最小となる補助関数Φ+(Θ，ｍ)の最小値を示す。
また、パラメータΘ＝Θ１，・・・，ΘI に関する最小化を行い式（４０）に示す。
【０１１８】
【数４６】
2010039215.tif 000048
【０１１９】
式（４０）において、左辺に示した「Θ」の上に「＾」を付けた文字で示される変数（以下、「＾Θ」と記載する）は、式（３９）によって得られた補助変数＾ｍを固定にして補助変数Θに関して最小となる補助関数Φ+(Θ，＾ｍ)の最小値を示す。
それぞれ、最小化を行うと目的関数Φ(Θ)は単調減少する。Φ(Θ)が下に有界ならば、上記の一連の更新の反復計算によりΘは停留点に収束する。
【０１２０】
証明：任意のパラメータをΘ（０）とし、式（３９）、式（４０）によりΦ(Θ)が減少することを示す。まず式（３９）により、式（４１）と示すことができる。
【０１２１】
【数４７】
2010039215.tif 000049
【０１２２】
また、式（４０）に示される更新により明らかに式（４２）の関係を示すことができる。
【０１２３】
【数４８】
2010039215.tif 000050
【０１２４】
補助関数の定義を示す式(３８)より、式（４３）の関係を示すことができる。
【０１２５】
【数４９】
2010039215.tif 000051
【０１２６】
ゆえに、式（４４）に示す関係がある。
【０１２７】
【数５０】
2010039215.tif 000052
【０１２８】
以上より、各更新による目的関数の非増加が示された。（証明終）。
【０１２９】
この目的関数更新方法をある最適化問題に応用するにあたり、不等式（３８）の統合が成立するときの補助変数が解析的に求められるかどうか、補助関数を最小化するモデルパラメータが解析的に求められるかが鍵であり、そのようになるまでうまく補助関数を設計することが重要である。
【０１３０】
式（３６）の補助関数は次の定理により導かれる。
【０１３１】
＜定理３＞
Φ(Θ)が式（３６）で与えられるとき、次式（４５）の不等式が成り立つ。
【０１３２】
【数５１】
2010039215.tif 000053
【０１３３】
式（４５）において、右辺の第４項の積分範囲を示す式（４６）は、集合Ｒから集合Ｓを除いた差集合を示すものとする。
【０１３４】
【数５２】
2010039215.tif 000054
【０１３５】
ただし、式（４５）において、式（４７）とする。
【０１３６】
【数５３】
2010039215.tif 000055
【０１３７】
また、式（４７）において、λｋ，ｎはΣｋ，ｎλｋ，ｎ＝１、λ∈（０，１）を満たす任意の定数とする。
式（４７）の不等式の等号は、式（４８）、式（４９）に示される条件のときに成立する。
【０１３８】
【数５４】
2010039215.tif 000056
【０１３９】
【数５５】
2010039215.tif 000057
【０１４０】
証明：定理３は、次に示す２つの補題により直ちに証明できる。
【０１４１】
＜補題１＞
ｚ∈（−∞，∞）においてある複素数値関数ωｉ（ｚ）が式（５０）を満たすとき、式（５１）の関係が成り立つ。
【０１４２】
【数５６】
2010039215.tif 000058
【０１４３】
【数５７】
2010039215.tif 000059
【０１４４】
式（５１）において、式（５２）の関係で示されるとき等号が成立する。
【０１４５】
【数５８】
2010039215.tif 000060
【０１４６】
ただし、式（５２）において、βｉはβｉ∈（０，１）、Σｉβｉ＝１を満たす任意の定数である。
【０１４７】
証明：
Ｌａｇｒａｎｇｅ未定乗数法を用いた変分法により、不等式の右辺のωｉ（ｚ）に関する最小値が左辺と等しくなることを示す。また、そのときのωｉ（ｚ）を求める。
そこで、まず、不等式の右辺に式（５０）で示される条件のためのＬａｇｒａｎｇｅ未定乗数項を付加した式（５３）で示される汎関数の極値を求めることとする。
【０１４８】
【数５９】
2010039215.tif 000061
【０１４９】
式（５３）において、Ｊ［ω］の
【０１５０】
【数６０】
2010039215.tif 000062
【０１５１】
に関する変分は、式（５４）として示される。
【０１５２】
【数６１】
2010039215.tif 000063
【０１５３】
式（５４）が恒等的に０であるためには、式（５５）で示される関係が必要である。
【０１５４】
【数６２】
2010039215.tif 000064
【０１５５】
よって、式（５６）を０と置くと、式（５７）を得る。
【０１５６】
【数６３】
2010039215.tif 000065
【０１５７】
【数６４】
2010039215.tif 000066
【０１５８】
式（５０）に示される条件により、式（５８）である。
【０１５９】
【数６５】
2010039215.tif 000067
【０１６０】
式（５８）を式（５７）に代入すると、最終的に式（５２）を得る。これを式（５１）の右辺に代入すると、左辺と等しくなる。次に、この極値が極小解であるかどうかは、Ｊ［ω］のω（ｚ）に関して式（５９）に示されるＨｅｓｓｅ（ヘッセ）行列が正定値であることから示される。
【０１６１】
【数６６】
2010039215.tif 000068
【０１６２】
（証明終）。
【０１６３】
＜補題２＞
式（６０）に示す条件において、式（６１）に示す不等式が成り立ち、等号はｚ＝αにおいて成立する。
【０１６４】
【数６７】
2010039215.tif 000069
【０１６５】
【数６８】
2010039215.tif 000070
【０１６６】
証明：実関数ｆ（ｘ）を微分可能な凹関数とする。凹関数の定義より、任意の異なる２点
【０１６７】
【数６９】
2010039215.tif 000071
【０１６８】
は、任意の実数γ∈（０，１）に対して、式（６２）が成り立つとき、ｆ（ｘ）は凹関数という。
【０１６９】
【数７０】
2010039215.tif 000072
【０１７０】
式（６２）は、式（６３）に書き換えることができる。
【０１７１】
【数７１】
2010039215.tif 000073
【０１７２】
また、式（６３）は、微分可能であるので、実数γが０に漸近（γ→０）するもとでは、式（６４）に示すことができる。
【０１７３】
【数７２】
2010039215.tif 000074
【０１７４】
式（６４）を式（６３）に代入すると、式（６５）を得る。
【０１７５】
【数７３】
2010039215.tif 000075
【０１７６】
−ｅ−ｚは、ｚの凹関数であるので、ｆ（ｚ）を式（６６）で定義して、式（６５）に代入すると、式（６７）を得る。
【０１７７】
【数７４】
2010039215.tif 000076
【０１７８】
【数７５】
2010039215.tif 000077
【０１７９】
また、式（６７）において明らかにｚ＝ｈにおいて等号は成立する（証明終）。
【０１８０】
続いて、定理２における式（４０）を導出する。
式（４５）により得られた補助関数Φ＋（Θ，ｍ，α）を最小化するΩｋは厳密には求めることができない。しかし、もし、この補助関数において、Ωｋの非多項式項である式（６８）の目的関数に対する寄与が、式（６９）に示される式に比べて著しく小さくなるという傾向があるならば、補助関数Φ＋（Θ，ｍ，α）のうち式（６８）に示した項を無視することにより近似された式を式（７０）に示す。
【０１８１】
【数７６】
2010039215.tif 000078
【０１８２】
【数７７】
2010039215.tif 000079
【０１８３】
【数７８】
2010039215.tif 000080
【０１８４】
式（７０）において、
【０１８５】
【数７９】
2010039215.tif 000081
【０１８６】
は、補助関数Φ＋（Θ，ｍ，α）の近似式であり「〜Φ＋（Θ，ｍ，α）」と示す。
式（７０）に示される〜Φ＋（Θ，ｍ，α）を最小化できるΩｋは、補助関数Φ＋（Θ，ｍ，α）を最小化するΩｋの良い近似解になる。そこで、式（７０）をΩｋにより偏微分して式（７１）に示す。
【０１８７】
【数８０】
2010039215.tif 000082
【０１８８】
式（７１）の値を、０と置くと、式（７２）に示すようにΩｋについて整理できる。
【０１８９】
【数８１】
2010039215.tif 000083
【０１９０】
また、補助関数Φ＋（Θ，ｍ，α）の
【０１９１】
【数８２】
2010039215.tif 000084
【０１９２】
に関する偏微分を０と置いて解くと、式（７３）を得る。
【０１９３】
【数８３】
2010039215.tif 000085
【０１９４】
さて、αＫ，ｎ（ｘ）を式（４９）のとおりに更新した直後にＢＫ，ｎを更新する場合を考えると、式（７３）は、式（７４）となることに注意する。
【０１９５】
【数８４】
2010039215.tif 000086
【０１９６】
すなわち、ＢＫ，ｎの更新値は、式（７５）に示す式の重みつき平均に他ならない。
【０１９７】
【数８５】
2010039215.tif 000087
【０１９８】
したがって、ｍｋ，ｎ（ｘ）Ｗｙ（ｘ）が優勢である（絶対値が大きい）ｘの大部分において、式（７４）にしたがって更新されたＢＫ，ｎとｍｋ，ｎ（ｘ）Ｗｙ（ｘ）は複素平面上で鋭角をなしており、式（７６）は正になっていることが予想される。仮に式（７６）が負になる対数周波数ｘがあったとしても、その対数周波数ｘでのｍｋ，ｎ（ｘ）Ｗｙ（ｘ）は、優勢ではない。
【０１９９】
【数８６】
2010039215.tif 000088
【０２００】
以上のようにして、Ωｋの更新式を導く上で補助関数に関して立てた仮定（補助関数に対して式（６８）の寄与が小さいこと）は、実際の傾向としてある程度の妥当性をもって成立すると見ることができる。
【０２０１】
図を参照し、擬似周期信号重畳モデルのパラメータ最適化アルゴリズムを用いた信号処理の手順について説明する。
図２は、擬似周期信号重畳モデルのパラメータ最適化アルゴリズムを用いた信号処理を示すフローチャートである。
このフローチャートで示される信号処理は、以下に説明する複数の信号処理過程に従って処理される。この信号処理は、信号処理装置１００におけるＣＰＵで実行されるプログラムによって処理される。
【０２０２】
観測時間周波数成分出力過程では、観測時間周波数成分出力部１０は、時系列信号ｙ（ｔ）を入力信号として、マザーウェーブレットの周波数特性が対数正規分布型の関数Ψ＊（ω）によって定義されたウェーブレット変換によって変換された観測時間周波数成分Ｗｙ（ｘ，ｔ）を出力する（ステップＳａ１）。
【０２０３】
初期パラメータ生成過程では、初期パラメータ生成部２０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）を入力として離散時刻ごとの対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）とを含む擬似周期信号モデルパラメータの初期値および誤差配分比重定数λｋ，ｎ（ｔ）を設定する（ステップＳａ２）。
【０２０４】
補助変数更新過程では、補助変数更新部３０は、観測時間周波数成分出力部１０によって出力された観測時間周波数成分Ｗｙ（ｘ，ｔ）および、ステップＳａ２またはステップＳａ５で出力された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、誤差配分比重定数をλｋ，ｎ（ｔ）を入力として、観測時間周波数成分を各擬似周期信号モデルの各調波成分に対応する成分に分配するために用いられる分配フィルタｍｋ，ｎ（ｘ，ｔ）と各分配成分に重みづけするために用いられるスペクトル重みαｋ，ｎ（ｘ，ｔ）とを含む補助変数を更新する（ステップＳａ３）。
【０２０５】
周波数値集合更新過程では、周波数値集合更新部４０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、ステップＳａ２またはステップＳａ５で出力された擬似周期信号モデルパラメータの要素である調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と補助変数更新部３０によって出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）を入力とされる。また、周波数値集合更新部４０は、その入力によって、離散時間ごとに擬似周期信号モデルの各調波成分の調波成分複素振幅パラメータＢｋ，ｎ（ｔ）とそれに対応する観測時間周波数成分の分配成分ｍｋ，ｎ（ｘ，ｔ）Ｗｙ（ｘ，ｔ）とが複素平面上で鋭角をなす全ての周波数値を含む周波数値集合Ｓ（ｔ）ｋ，ｎを出力する（ステップＳａ４）。
【０２０６】
擬似周期信号モデルパラメータ更新過程では、擬似周期信号モデルパラメータ更新部５０は、観測時間周波数成分Ｗｙ（ｘ，ｔ）と、ステップＳａ２または前回のステップＳａ５で出力された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）と、補助変数更新部３０で出力された補助変数の要素である分配フィルタｍｋ，ｎ（ｘ，ｔ）とスペクトル重みαｋ，ｎ（ｘ，ｔ）と、周波数値集合Ｓ（ｔ）ｋ，ｎと誤差配分比重定数λｋ，ｎ（ｔ）を入力とする。また、擬似周期信号モデルパラメータ更新部５０は、擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）を出力する（ステップＳａ５）。
【０２０７】
収束判定過程では、収束判定部６０は、擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には補助変数更新部３０からの処理を繰り返し行わせる。収束判定部６０における所定の規準を満たすか否かの判定基準は、例えば、式（１６）に示されるように、予め定められる目的関数によって定められる目的関数値の変化率が予め定められた所定値以下となることとすることを基準とする。目的関数を定める式は、前述の式（１６）に示す（ステップＳａ６）。
【０２０８】
擬似周期信号モデルパラメータ出力過程では、擬似周期信号モデルパラメータ出力部７０は、収束判定部６０において所定の規準を満たすと判定された擬似周期信号モデルパラメータの要素である対数基本周波数パラメータΩｋ（ｔ）と調波成分複素振幅パラメータＢｋ，ｎ（ｔ）を出力する（ステップＳａ７）。
なお、上記のステップＳａ３、ステップＳａ４，ステップＳａ５における繰り返して演算される際に参照される擬似周期信号モデルパラメータは、初回の演算では、ステップＳａ２の処理で生成された初期値を参照し、ステップＳａ６の判定によって繰り返されることになる２回目以降の演算では、ステップＳａ５によって生成され更新された擬似周期信号モデルパラメータを参照することとする。
【０２０９】
図３を参照し本発明の擬似周期信号重畳モデルのパラメータ最適化アルゴリズムを用いた信号処理（以下、本発明の信号処理という。）を用いた推定性能について示す。
ここでは、基本周波数が非常に近接した２つの周期信号の合成音（１６ｋＨｚ)を分析対象として、擬似周期信号モデルパラメータの推定処理を行う。基本周波数が非常に近接した２つの周期信号は、一方のピッチ（１００Ｈｚ，４００Ｈｚ，８００Ｈｚ，)に対して、もう一方のピッチが約２（１／１２）倍である合成音（１６ｋＨｚ)とする。
さらに、本発明の信号処理による手法と「非特許文献２」における従来方式による手法での性能を比較して示す。
「非特許文献２」における従来方式による手法は、短時間Ｆｏｕｒｉｅｒ（フーリエ）変換（以下、ＳＴＦＴ（Short Time Fourier Transformation）という）モデルを用いて、式（６）に示した目的関数をＧａｂｏｒ（ガボール）変換（ＳＴＦＴ）領域で等価表現し、これを規準として擬似周期信号のＳＴＦＴモデルのパラメータをフレームごとに推定するアルゴリズムを用いたものであり、観測信号の各周波数成分に対し、高周波および低周波の正弦波パラメータを同じ長さの短時間窓内でフィッティングさせる手法を示す。
【０２１０】
信号処理の対象とする正弦波重畳モデルにおいてＫ＝２，Ｎ＝１０とし、本発明の信号処理による手法は６０Ｈｚ（ヘルツ）〜８ｋＨｚ（キロヘルツ）の分析帯域に対して、５１２チャンネルのフィルタバンクを設定し各フィルタバンクでの分析を行った。
また、参考例として比較対象とした「非特許文献２」における従来方式による手法では、Ｇａｂｏｒ変換の実効フレーム長さは１０２４点とした。
【０２１１】
それぞれの周波数組合せに対して、一回の推定毎に、入力信号の１０個の調波成分の振幅と位相を乱数で発生させ合成信号を作りなおし、初期値も様々に変化させて、それぞれ５０回ずつ推定を行った。各推定における、アルゴリズムの反復演算回数は１００回とした。また、周波数に依存しない推定性能に基づいて比較するため、推定される周波数の推定精度を、正解とされる基準周波数を基準にしたセント値を用いて示すこととする。例えば、周波数ωのセント値は、基準周波数ω０に対して、１２００ｌｏｇ１０（ω／ω０）で算出される値である。基準周波数に対して推定された周波数の精度が、５０セント以内に収まるものを、正しく推定が行われたものとして正解とみなし、５０回ずつの推定における正解率を算出した。
【０２１２】
本発明の信号処理による手法によって、設定された基本周波数が（８００Ｈｚ，８４７Ｈｚ)の分析では、初期値依存性が少なく安定的かつ効率的な推定性能が実現できることが示される。また、設定された基本周波数が（１００Ｈｚ，１０５Ｈｚ)と（４００Ｈｚ，４２３Ｈｚ)の分析でも、周波数に依存することなく、ほぼ一定の性能を得ることができることが示される。
参考に比較を行った「非特許文献２」における従来手法では、設定された基本周波数が（８００Ｈｚ，８４７Ｈｚ)の分析では本発明の信号処理による手法との差異が少ないものの、周波数が低く設定された基本周波数が（１００Ｈｚ，１０５Ｈｚ)と（４００Ｈｚ，４２３Ｈｚ)の分析では、性能が著しく低下することが示される。
本発明は、上記比較からもわかるように従来手法による性能限界を改善して、観測信号の基本周波数が低い時系列信号においても、モデルパラメータの低周波成分の推定精度を確保できる擬似周期信号のパラメータを推定することができ、周波数に依存しない処理性能を実現できる信号処理装置１００を提供することができる。
【０２１３】
（第２実施形態）
第２実施形態における信号処理装置２００による信号処理において異なる収束判定処理について示す。
ここで示す収束判定処理は、繰り返して行われる擬似周期信号モデルパラメータの反復演算において擬似周期信号モデルパラメータが推定値近傍の所定の範囲内に収束するまで処理することとする。
図４は、信号処理装置２００を示すブロック図である。
信号処理装置２００のうち、収束判定処理部６２以外の構成は第１実施形態に係る信号処理装置１００と同じであるため、同じ符号を付し、第１実施形態における説明を参照し、収束判定処理部６２について説明する。
信号処理装置２００における収束判定部６２は、擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かを判定し、満たさない場合には補助変数更新部３０からの処理を繰り返し行わせる。収束判定部６２における所定の規準を満たすか否かの判定基準は、擬似周期信号モデルパラメータが反復して演算されて更新される擬似周期信号モデルパラメータの変化率が予め定められた所定値以下となることを基準とする。
これにより、繰り返し行われる擬似周期信号モデルパラメータの演算処理で、擬似周期信号モデルパラメータが推定値近傍に収束して、擬似周期信号モデルパラメータの変化率が予め定められた所定値以下になった段階で推定処理を打ち切ることにより、演算時間を短縮することができ、効率よく擬似周期信号のパラメータを推定する信号処理装置２００を実現することができる。
【０２１４】
（第３実施形態）
第３実施形態における信号処理装置３００による信号処理にいて、異なる収束判定処理について示す。
ここで示す収束判定処理は、繰り返して行われる擬似周期信号モデルパラメータの反復演算回数が予め定められた所定の回数まで処理することとする。
図５は、信号処理装置３００を示すブロック図である。
信号処理装置３００のうち、収束判定処理部６３以外の構成は第１実施形態に係る信号処理装置１００と同じであるため、同じ符号を付し、第１実施形態における説明を参照し、収束判定処理部６３について説明する。
擬似周期信号モデルパラメータの演算回数が演算回数記憶部に記録され、その初期値を１として設定される。また、演算回数の限度値記憶部に予め定められた回数を示す値が記憶される。
信号処理装置３００における収束判定部６３は、演算回数記憶部と限度値記憶部を参照し、参照して得られた演算回数と限度値とを比較する。また、演算回数の比較を行うごとに演算回数記憶部に記録されている値に１加えて演算回数記憶部の値を更新する。
【０２１５】
信号処理装置３００における収束判定部６３は、擬似周期信号モデルパラメータの演算結果の判定において所定の規準を満たすか否かが判定され、満たさない場合には補助変数更新部３０からの処理を繰り返し行わせる。収束判定部６３における所定の規準を満たすか否かの判定基準は、繰り返して行われる擬似周期信号モデルパラメータの反復演算回数が予め定められた所定の回数（限度値）を満たすことを判定基準とする。
これにより、繰り返し行われる擬似周期信号モデルパラメータの演算処理の反復演算回数を予め定められた所定の回数で打ち切ることで演算時間を所定の時間内とすることができ、効率よく擬似周期信号のパラメータを推定する信号処理装置３００を実現することができる。
【０２１６】
本発明の信号処理装置は、本発明の実施形態における信号処理装置１００、信号処理装置２００、信号処理装置３００に対応する。本発明の観測時間周波数成分出力部は、本発明の実施形態における観測時間周波数成分出力部１０に対応する。本発明の初期パラメータ生成部は、本発明の実施形態における初期パラメータ生成部２０に対応する。本発明の補助変数更新部は、本発明の実施形態における補助変数更新部３０に対応する。本発明の周波数集合更新部は、本発明の実施形態における周波数集合更新部４０に対応する。本発明の擬似周期信号モデルパラメータ更新部は、本発明の実施形態における擬似周期信号モデルパラメータ更新部５０に対応する。本発明の収束判定部は、本発明の実施形態における収束判定部６０、収束判定部６２、または収束判定部６３に対応する。本発明の擬似周期信号モデルパラメータ出力部は、本発明の実施形態における擬似周期信号モデルパラメータ出力部７０に対応する。
【０２１７】
本発明の入力される時系列信号を擬似周期信号とみなすステップは、本発明の実施形態における観測時間周波数成分出力過程に対応する。本発明の擬似周期信号の特徴を示すパラメータを推定するステップは、本発明の実施形態における初期パラメータ生成過程、補助変数更新過程、周波数値集合更新過程、擬似周期信号モデルパラメータ更新過程、収束判定過程および擬似周期信号モデルパラメータ出力過程に対応する。
【０２１８】
なお、本発明は、上記の各実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で変更可能である。本発明の信号処理装置における、収束判定部では、必要に応じて第１実施形態から第３実施形態に示した判定処理を複数あわせて備えることができる。
【図面の簡単な説明】
【０２１９】
【図１】本発明の第１実施形態における信号処理装置を示すブロック図である。
【図２】第１実施形態における信号処理手順を示すフローチャートである。
【図３】第１実施形態における信号処理性能を示す図である。
【図４】本発明の第２実施形態における信号処理装置を示すブロック図である。
【図５】本発明の第３実施形態における信号処理装置を示すブロック図である。
【符号の説明】
【０２２０】
１００  信号処理装置
１０  観測時間周波数成分出力部
２０  初期パラメータ生成部
３０  補助変数更新部
４０  周波数値集合更新部
５０  擬似周期信号モデルパラメータ更新部
６０  収束判定部
７０  擬似周期信号モデルパラメータ出力部
【図１】
2010039215.tif 000089
【図２】
2010039215.tif 000090
【図３】
2010039215.tif 000091
【図４】
2010039215.tif 000092
【図５】
2010039215.tif 000093