(57)【要約】【課題】正弦波重畳モデルの最尤パラメータを求めるための新しい推定アルゴリズムを導出する。【解決手段】入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する。まず、入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として設定する。次に、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求める。そして、前記補助変数ｍの更新を行うステップと、前記パラメータΘの更新を行うステップと、を交互に実行する。これは、前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで繰り返す。このように補助関数を用いて、順に更新していくことによって、正弦波重畳モデルのパラメータΘが求められる。【選択図】図１
2008304718.tif 000002
【特許請求の範囲】
【請求項１】
入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する方法において、
入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として設定するステップと、
前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップと、
前記補助変数ｍの更新を行うステップと、
前記パラメータΘの更新を行うステップと、
前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行うステップ及び前記パラメータΘの更新を行うステップ、を繰り返し実行する繰り返しステップと、
を含むことを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項２】
請求項１記載の正弦波重畳モデルパラメータ推定方法において、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップは、更に、以下の要件１及び要件２を満たす補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めることを特徴とする正弦波重畳モデルパラメータ推定方法。
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
【請求項３】
請求項１記載の正弦波重畳モデルパラメータ推定方法において、
正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数１】
2008304718.tif 000199
とし、それ以外のパラメータを要素としたベクトルを
【数２】
2008304718.tif 000200
とするとき、補助関数Φ＋（Θ、ｍ）が、
【数３】
2008304718.tif 000201
のような形態である正弦波重畳モデルパラメータ推定方法。
【請求項４】
請求項２記載の正弦波重畳モデルパラメータ推定方法において、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めることを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項５】
請求項１記載の正弦波重畳モデルパラメータ推定方法において、
前記目的関数を設定するステップは、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分を目的関数Φ（Θ）として設定することを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項６】
請求項５記載の正弦波重畳モデルパラメータ推定方法において、
前記前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、
前記パラメータΘの更新を行うステップは、前記ｄも更新することを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項７】
請求項１記載の正弦波重畳モデルパラメータ推定方法において、
前記補助関数を設定するステップは、ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成し、その第２の補助関数を、前記目的関数Φ（Θ）に対する補助関数として設定し、
前記補助変数の更新を行うステップは、前記補助変数ｍとともに前記補助変数αも更新することを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項８】
請求項１記載の正弦波重畳モデルパラメータ推定方法において、
前記繰り返しステップにおいて、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、パラメータの推定を終了することを特徴とする正弦波重畳モデルパラメータ推定方法。
【請求項９】
複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離方法において、
前記入力音声信号を、所定フレームに区分けするステップと、
前記区分けしたフレーム毎に、前記請求項１〜８のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定するステップと、
前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得るステップと、
を含むことを特徴とする音源分離方法。
【請求項１０】
入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として用い、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを用い、前記入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する装置において、
前記補助変数ｍの更新を行う手段と、
前記パラメータΘの更新を行う手段と、
前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行う手段、及び、前記パラメータΘの更新を行う手段、を順番に動作させる手段と、
を含み、前記順番に動作させる手段は、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、動作を終了し、その時点でのパラメータを出力することを特徴とする正弦波重畳モデルパラメータ推定装置。
【請求項１１】
請求項１０記載の正弦波重畳モデルパラメータ推定装置において、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、以下の要件１及び要件２を満たすことを特徴とする正弦波重畳モデルパラメータ推定装置。
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
【請求項１２】
請求項１０記載の正弦波重畳モデルパラメータ推定装置において、
正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数４】
2008304718.tif 000202
とし、それ以外のパラメータを要素としたベクトルを
【数５】
2008304718.tif 000203
とするとき、前記補助関数Φ＋（Θ、ｍ）が、
【数６】
2008304718.tif 000204
のような形態である正弦波重畳モデルパラメータ推定装置。
【請求項１３】
請求項１１記載の正弦波重畳モデルパラメータ推定装置において、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍであることを特徴とする正弦波重畳モデルパラメータ推定装置。
【請求項１４】
請求項１０記載の正弦波重畳モデルパラメータ推定装置において、
前記目的関数Φ（Θ）は、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分であることを特徴とする正弦波重畳モデルパラメータ推定装置。
【請求項１５】
請求項１４記載の正弦波重畳モデルパラメータ推定装置において、
前記前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、
前記パラメータΘの更新を行う手段は、前記ｄも更新することを特徴とする正弦波重畳モデルパラメータ推定装置。
【請求項１６】
請求項１０記載の正弦波重畳モデルパラメータ推定装置において、
ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成することによって、その第２の補助関数を前記補助関数として用い、
前記補助変数の更新を行う手段は、前記補助変数ｍとともに前記補助変数αも更新することを特徴とする正弦波重畳モデルパラメータ推定装置。
【請求項１７】
複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離装置において、
前記入力音声信号を、所定フレームに区分けする手段と、
前記区分けしたフレーム毎に、前記請求項１〜８のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定する手段と、
前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得る手段と、
を含むことを特徴とする音源分離装置。
【請求項１８】
コンピュータを、入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として用い、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを用い、前記入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する装置として動作させるプログラムにおいて、
前記補助変数ｍの更新を行う手順と、
前記パラメータΘの更新を行う手順と、
前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行う手順、及び、前記パラメータΘの更新を行う手順、を順番に動作させる手順と、
を、コンピュータに実行させ、前記順番に動作させる手順において、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、動作を終了し、その時点でのパラメータを出力することを特徴とするプログラム。
【請求項１９】
請求項１８記載のプログラムにおいて、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、以下の要件１及び要件２を満たすことを特徴とするプログラム。
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
【請求項２０】
請求項１８記載のプログラムにおいて、
正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数７】
2008304718.tif 000205
とし、それ以外のパラメータを要素としたベクトルを
【数８】
2008304718.tif 000206
とするとき、前記補助関数Φ＋（Θ、ｍ）が、
【数９】
2008304718.tif 000207
のような形態であるプログラム。
【請求項２１】
請求項１９記載のプログラムにおいて、
前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍであることを特徴とするプログラム。
【請求項２２】
請求項１８記載のプログラムにおいて、
前記目的関数Φ（Θ）は、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分であることを特徴とするプログラム。
【請求項２３】
請求項２２記載のプログラムにおいて、
前記前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、
前記パラメータΘの更新を行う手順においては、前記ｄも更新することを特徴とするプログラム。
【請求項２４】
請求項１８記載のプログラムにおいて、
ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成することによって、その第２の補助関数を、前記コンピュータに前記補助関数として用いさせ、
前記補助変数の更新を行う手順は、前記補助変数ｍとともに前記補助変数αも更新することを特徴とするプログラム。
【請求項２５】
コンピュータを、複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離装置として動作させるプログラムにおいて、
前記入力音声信号を、所定フレームに区分けする手順と、
前記区分けしたフレーム毎に、前記請求項１〜８のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定する手順と、
前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得る手順と、
を、前記コンピュータに実行させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、正弦波重畳モデルのパラメータ推定及びその応用に関する。
【背景技術】
【０００２】
正弦波重畳モデル（下記非特許文献１）は、ＭｃＡｕｌａｙらが分析合成系に適用し、高品質な合成音声を生成できることを示して以来、分析合成系に留まらず、テキスト音声合成、音声加工、符号化などの用途として広く用いられてきた。特に、高品質な合成音声を生成できることは、正弦波重畳モデルが音声や楽音などの音響信号を極めて良く表現するモデルであることの根拠となっているため、音源分離への応用に対する期待も大きい。一方、音声の自己相関関数を少数の正弦波の重ね合わせで表現しようとする複合正弦波モデル（下記非特許文献２１）もまた良好な音声符号化方式として知られている。応用場面に依らずこの枠組（正弦波重畳モデルを用いた信号分析）において共通するのは、いかにして正弦波重畳モデルのパラメータを精度良く推定するかが最も重要な問題となっている点であり、その推定精度があらゆる応用における性能に直結する。
【０００３】
さて、ＭｃＡｕｌａｙらが採用した正弦波重畳モデルは、周波数と振幅がいずれも定常であることを仮定したＫ個の周波数成分を重畳した解析信号（式（１））
【数１０】
2008304718.tif 000003
である。ただし、
【数１１】
2008304718.tif 000004
は、それぞれｋ番目の正弦波成分の周波数、複素振幅を表す。短時間分析区間
【数１２】
2008304718.tif 000005
における対象信号の解析信号を
【数１３】
2008304718.tif 000006
とするとき、
【数１４】
2008304718.tif 000007
が、下記式（２）
【数１５】
2008304718.tif 000008
と表され、
【数１６】
2008304718.tif 000009
がＧａｕｓｓ性白色雑音
【数１７】
2008304718.tif 000010
であるという仮定のもとで
【数１８】
2008304718.tif 000011
の最尤パラメータを求めることが具体的な問題設定である。この場合、
【数１９】
2008304718.tif 000012
であるため、
【数２０】
2008304718.tif 000013
の最尤解は、Ｌ２ノルム誤差を最小にする解（式（３））
【数２１】
2008304718.tif 000014
に相当する。ただし、ここでは
【数２２】
2008304718.tif 000015
を
【数２３】
2008304718.tif 000016
において１、それ以外では０をとる矩形窓としている。
【０００４】
上述した式（１）（数１０）が示すように、正弦波重畳モデルは、
【数２４】
2008304718.tif 000017
に関しては線形に依存するが、
【数２５】
2008304718.tif 000018
に関しては非線形に依存している。したがって、
【数２６】
2008304718.tif 000019
が固定の場合には
【数２７】
2008304718.tif 000020
の最尤解が解析的に得られることはただちに分かるが、この
【数２８】
2008304718.tif 000021
が固定であったとしても
【数２９】
2008304718.tif 000022
の最尤解は解析的に得られない。
【０００５】
この点が、正弦波重畳モデルのパラメータ最適化における難しさの本質である。そのため、統計的信号処理の分野では
【数３０】
2008304718.tif 000023
の最尤解を得るための方法がこれまで長年に渡り数多く研究されている。そのような研究は、例えば、下記非特許文献１〜非特許文献１４に見られる。
【０００６】
さて、最初に述べたＭｃＡｕｌａｙらの手法では
【数３１】
2008304718.tif 000024
の推定値を得るのに、対象信号の離散パワースペクトル密度（ピリオドグラム）が最大となるピーク成分の周波数と振幅と位相を決定し、そのピーク成分を信号から減算する操作をＫ回繰り返すという簡便な方法が用いられた（下記非特許文献１参照）。
【０００７】
単一正弦波のピリオドグラムの最大値を与える周波数が最尤推定量であり、かつその推定量が不遍推定量であるとするＲｉｆｅらの指摘（下記非特許文献１７〜下記非特許文献１９）が上記の周波数推定法が式（３）（数２１）の近似解であることの正当性に対する１つの根拠であると言える。
【０００８】
しかし、
（ａ）離散的なピリオドグラムのピークが真の連続ピリオドグラムの最大値とは必ずしも対応しないこと、
（ｂ）周波数成分が複数ある場合には周波数成分間の干渉のため上記理論が成立しなくなること、
（ｃ）複数の周波数成分が近接する場合にはエネルギー拡散のため各々のピークが正しく検出できないことが起こり得ること、
などを考慮するとこの単純な周波数推定法と比較してより一層高精度な推定法の開発が望まれることは、世の中の技術進歩の流れを考慮すれば当然のこととも考えられよう。
【０００９】
そのような中で、ピーク成分のパラメータを直接推定値として扱うのではなく、ピーク近傍の複数点を結ぶ補間曲線から極大点を求めてパラメータ推定の高精度化を試みたもの（非特許文献４〜非特許文献９）がその簡便さから近年特によく用いられ、数多く提案されている。
【００１０】
しかしながら、これらピーク近傍の複数点を結ぶ補間曲線から極大点を求めてパラメータ推定を行う手法（例えば、非特許文献４〜非特許文献９）は、上記の（ｂ）と（ｃ）の問題は依然解決されておらず、ＭｃＡｕｌａｙらの方法と同様、上記式（３）（数２１）の近似解を与えるに過ぎない。そのため、これらの研究では、その近似精度を向上するための零詰め量や窓関数の設計方法などが、議論の対象の中心となっている（例えば、非特許文献８、非特許文献９）。
【００１１】
その一方で、上記式（３）（数２１）の解を最急降下法やＮｅｗｔｏｎ法などの勾配法、あるいはＧｉｂｂｓサンプラーや、マルコフ連鎖モンテカルロ（ＭＣＭＣ）法などの確率的サンプリング手法のような非線形最適化法により数値的に探索する方法も提案されている（例えば、非特許文献３、非特許文献１０〜１４）。
【００１２】
ところで、ＭｃＡｕｌａｙらのモデルがＫ個の純音信号を混合したものであるのに対して、Ｋ個の調波信号（Ｎ個の調波成分で構成され、ｎ次調波成分の周波数が基本周波数
【数３２】
2008304718.tif 000025
のｎ倍であるような信号）を重畳した解析信号（式（４））
【数３３】
2008304718.tif 000026
も同様に考えることができる。このモデルは、対象の混合信号が調波信号のみから構成されているような場合の１ｃｈ音源分離におけるモデルとしてしばしば用いられる（非特許文献１１〜非特許文献１５参照）。上記式（４）（数３３）においてＮ＝１としたとき、式（１）と同じモデルを想定したことになるため、ＭｃＡｕｌａｙらのモデルは式（４）の特殊なケースと見ることができる。
【００１３】
しかしこの式（４）のモデルでは、各正弦波の周波数がそれぞれ自由度をもたず、
【数３４】
2008304718.tif 000027
のような拘束関係をもつため、
【数３５】
2008304718.tif 000028
の最尤解を得るのは上記のＭｃＡｕｌａｙらのモデルの場合より更に難しくなる。例えば、上記のうちの前者のアプローチのように、研究の中にはピーク抽出後の処理によって  
【数３６】
2008304718.tif 000029
を推定する手法を提案する研究もなされている。しかし、この手法では、抽出されたピークがどの調波信号の何番目の調波成分に対応するかを決定するための場当たり的な閾値判定に頼らざるを得なくなり、得られた
【数３７】
2008304718.tif 000030
の解の最尤性に関する議論が煩雑化する。
【００１４】
そのため、このモデルを用いた従来の音源分離法では、後者のアプローチ（勾配法やサンプリング手法）が採用されることが多い（下記非特許文献１１〜下記非特許文献１３参照）。しかしながら、このような数値計算法には、常に局所最適解の問題がつきまとう。
【００１５】
例えば、勾配法の場合には、収束に至るまでの反復計算を無数通りの初期点から行わない限り、また、確率的サンプリング手法の場合には無数回の反復試行を行わない限り、式（３）の大域的最適解を得られるという保証はない。そのため、少ない計算量でいかに解探索を行えるかが問題となる（計算量が少なく済むということはその分だけ多くの異なる初期値条件から解探索が行える）が、勾配法やサンプリング法のようないわば力任せな数値計算しか適用されていないのが現在の状況である。
【００１６】
先行特許文献 
下記特許文献１には、正弦波重畳モデルパラメータを、接続区間の間で漸次変化するように構成し、２つの原音性の周波数特性や基本周波数が異なっていても不連続感を近くさせない技術が開示されている。
【００１７】
また、下記特許文献２には、正弦波重畳合成法において、合成音声がブザー音的なるのを防止し得る技術が開示されている。
【００１８】
また、下記特許文献３には、独立成分分析を利用した周波数領域ブラインド音源分離法において、雑音を抑制する技術が開示されている。
【００１９】
【特許文献１】特開２００４−１０２１１８号公報
【特許文献２】特開２００２−０４１０７６号公報（特許第３５９２６１７号）
【特許文献３】特開２００３−０９９０９３号公報（特許第３８３１２２０号）
【非特許文献１】R. J. McAulay et al., IEEE Trans. Acoust., Speech, Signal Process., 34(4), 744, 1986.
【非特許文献２】D. C. Rife et al., Bell System Technical Journal, 55(9),1389, 1976.
【非特許文献３】T. J. Abatzoglou, IEEE Trans. Acoust., Speech, Signal Process., 33(1), 77, 1985.
【非特許文献４】J. O. Smith III et al., In Proc. ICMC'87, 290, 1987.
【非特許文献５】Ph. Depalle et al., In Proc. IEEE WASPAA'97, 1997.
【非特許文献６】Y. V. Zakharov et al., IEEE Electronic Letters, 35(19),1608, 1999.
【非特許文献７】S. S. Abeysekera, In Proc. IEEE Workshop on Statistical Signal Processing (SSP), 417, 2001.
【非特許文献８】J. C. Brown et al., J. Acoust. Soc. Am., 94(2), 662, 1993.
【非特許文献９】安部ら，信学技報，EA2004-61, 7, 2004.
【非特許文献１０】D. Starer et al., IEEE Trans. Acoust. Speech, Signal Process., 40(6), 1528, 1992.
【非特許文献１１】D. Chazan et al., In Proc. ICASSP'93, 2, 728, 1993.
【非特許文献１２】P. Jinachitra, In Proc. ICASSP'03, 6, 609, 2003.
【非特許文献１３】S. Godsill et al., In Proc. ICASSP'02, 2, 1769, 2002.
【非特許文献１４】T. Tolonen, In Proc. AES 106th Convention, 1999.
【非特許文献１５】T. Virtanen et al., In Proc. ICASSP'00, 2, 765, 2000.
【非特許文献１６】D. J. Thomson, In Proc. IEEE, 70(9), 1055, 1982.
【非特許文献１７】D. C. Rife et al., IEEE Trans. Info. Theory, 20(5), 591,1974.
【非特許文献１８】B. G. Quinn, IEEE Trans. Signal Process., 42, 1264, 1994.
【非特許文献１９】B. G. Quinn, IEEE Trans. Signal Process., 45(3), 814, 1997.
【非特許文献２０】M. Feder et al., IEEE Trans., Acoust., Speech, Signal Process., 36(4), 477, 1988.
【非特許文献２１】嵯峨山ら，通信学会誌，  J64-A(2), 105, 1981.
【発明の開示】
【発明が解決しようとする課題】
【００２０】
以上述べたように、正弦波重畳モデルは音声や楽音を極めて良く表現するモデルであるにもかかわらず、そのパラメータ推定法に関しては議論する余地が極めて大きい。このような背景のもと、正弦波重畳モデルの最尤パラメータを求めるための新しい推定アルゴリズムを導出することが本発明の目的である。
【課題を解決するための手段】
【００２１】
（１）本発明は、上記課題を解決するために、入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する方法において、前記入力信号と前記正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として設定するステップと、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップと、前記補助変数ｍの更新を行うステップと、前記パラメータΘの更新を行うステップと、前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行うステップ及び前記パラメータΘの更新を行うステップ、を繰り返し実行する繰り返しステップと、を含むことを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２２】
（２）また、本発明は、上記（１）記載の正弦波重畳モデルパラメータ推定方法において、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップは、更に、以下の要件１及び要件２を満たす補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めることを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２３】
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
（３）上記（１）１記載の正弦波重畳モデルパラメータ推定方法において、正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数３８】
2008304718.tif 000031
とし、それ以外のパラメータを要素としたベクトルを
【数３９】
2008304718.tif 000032
とするとき、補助関数Φ＋（Θ、ｍ）が、
【数４０】
2008304718.tif 000033
のような形態である正弦波重畳モデルパラメータ推定方法である。
【００２４】
（４）また、本発明は、上記（２）記載の正弦波重畳モデルパラメータ推定方法において、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めるステップは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを求めることを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２５】
（５）また、本発明は、上記（１）記載の正弦波重畳モデルパラメータ推定方法において、前記目的関数を設定するステップは、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分を目的関数Φ（Θ）として設定することを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２６】
（６）また、本発明は、上記（５）記載の正弦波重畳モデルパラメータ推定方法において、前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、前記パラメータΘの更新を行うステップは、前記ｄも更新することを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２７】
（７）また、本発明は、上記（１）記載の正弦波重畳モデルパラメータ推定方法において、前記補助関数を設定するステップは、ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成し、その第２の補助関数を、前記目的関数Φ（Θ）に対する補助関数として設定し、前記補助変数の更新を行うステップは、前記補助変数ｍとともに前記補助変数αも更新することを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２８】
（８）また、本発明は、上記（１）記載の正弦波重畳モデルパラメータ推定方法において、前記繰り返しステップにおいて、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、パラメータの推定を終了することを特徴とする正弦波重畳モデルパラメータ推定方法である。
【００２９】
（９）また、本発明は、上記課題を解決するために、複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離方法において、前記入力音声信号を、所定フレームに区分けするステップと、前記区分けしたフレーム毎に、前記（１）〜（８）のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定するステップと、前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得るステップと、を含むことを特徴とする音源分離方法。
【００３０】
（１０）また、本発明は、上記課題を解決するために、入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として用い、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを用い、前記入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する装置において、前記補助変数ｍの更新を行う手段と、前記パラメータΘの更新を行う手段と、前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行う手段、及び、前記パラメータΘの更新を行う手段、を順番に動作させる手段と、を含み、前記順番に動作させる手段は、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、動作を終了し、その時点でのパラメータを出力することを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３１】
（１１）また、本発明は、上記（１０）記載の正弦波重畳モデルパラメータ推定装置において、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、以下の要件１及び要件２を満たすことを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３２】
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
【００３３】
（１２）上記（１０）記載の正弦波重畳モデルパラメータ推定装置において、正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数４１】
2008304718.tif 000034
とし、それ以外のパラメータを要素としたベクトルを
【数４２】
2008304718.tif 000035
とするとき、前記補助関数Φ＋（Θ、ｍ）が、
【数４３】
2008304718.tif 000036
のような形態である正弦波重畳モデルパラメータ推定装置である。
（１３）本発明は、上記課題を解決するために、上記（１１）記載の正弦波重畳モデルパラメータ推定装置において、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍであることを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３４】
（１４）また、本発明は、上記（１０）記載の正弦波重畳モデルパラメータ推定装置において、前記目的関数Φ（Θ）は、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分であることを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３５】
（１５）また、本発明は、上記（１４）記載の正弦波重畳モデルパラメータ推定装置において、前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、前記パラメータΘの更新を行う手段は、前記ｄも更新することを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３６】
（１６）また、本発明は、上記（１０）記載の正弦波重畳モデルパラメータ推定装置において、ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成することによって、その第２の補助関数を前記補助関数として用い、前記補助変数の更新を行う手段は、前記補助変数ｍとともに前記補助変数αも更新することを特徴とする正弦波重畳モデルパラメータ推定装置である。
【００３７】
（１７）また、本発明は、上記課題を解決するために、複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離装置において、前記入力音声信号を、所定フレームに区分けする手段と、前記区分けしたフレーム毎に、前記（１）〜（８）のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定する手段と、前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得る手段と、を含むことを特徴とする音源分離装置である。
【００３８】
（１８）本発明は、上記課題を解決するために、コンピュータを、入力信号と正弦波重畳モデルとの２乗誤差の積分を目的関数Φ（Θ）として用い、前記目的関数Φ（Θ）に対して、目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）である補助関数Φ＋（Θ、ｍ）及びその補助変数ｍを用い、前記入力信号を正弦波重畳モデルで表現する場合のパラメータを推定する装置として動作させるプログラムにおいて、前記補助変数ｍの更新を行う手順と、前記パラメータΘの更新を行う手順と、前記補助変数又は前記パラメータΘの更新量が所定値未満になるまで、前記補助変数ｍの更新を行う手順、及び、前記パラメータΘの更新を行う手順、を順番に動作させる手順と、を、コンピュータに実行させ、前記順番に動作させる手順において、前記補助変数ｍ又は前記パラメータΘの更新量が所定値未満になったと判断された場合は、動作を終了し、その時点でのパラメータを出力することを特徴とするプログラムである。
【００３９】
（１９）本発明は、上記（１８）記載のプログラムにおいて、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、以下の要件１及び要件２を満たすことを特徴とするプログラムである。
【００４０】
（要件１）上記目的関数Φ（Θ）＝＜補助関数Φ＋（Θ、ｍ）という要件の等号が成立する際の前記補助変数ｍが解析的に得られること、
（要件２）前記補助関数Φ＋（Θ、ｍ）を最小化する前記パラメータΘが解析的に得られること。
（２０）上記（１８）記載のプログラムにおいて、正弦波重畳モデルの中の周波数又は基本周波数パラメータを要素にしたベクトルを
【数４４】
2008304718.tif 000037
とし、それ以外のパラメータを要素としたベクトルを
【数４５】
2008304718.tif 000038
とするとき、前記補助関数Φ＋（Θ、ｍ）が、
【数４６】
2008304718.tif 000039
のような形態であるプログラムである。
【００４１】
（２１）また、本発明は、上記（１９）記載のプログラムにおいて、前記補助関数Φ＋（Θ、ｍ）及びその補助変数ｍは、前記目的関数Φ（Θ）のパラメータのベクトルの２次式で表される補助関数Φ＋（Θ、ｍ）及びその補助変数ｍであることを特徴とするプログラムである。
【００４２】
（２２）本発明は、上記（１８）記載のプログラムにおいて、前記目的関数Φ（Θ）は、入力信号と正弦波重畳モデルとの２乗誤差に、ガウス窓関数を乗じた値の積分であることを特徴とするプログラムである。
【００４３】
（２３）また、本発明は、上記（２２）記載のプログラムにおいて、前記ガウス窓関数の窓の広がりを表す定数であるｄを、変数として扱い、前記パラメータΘの更新を行う手順においては、前記ｄも更新することを特徴とするプログラムである。
【００４４】
（２４）本発明は、上記（１８）記載のプログラムにおいて、ノルム誤差に基づき、前記目的関数Φ（Θ）から、第１の補助関数Φ＋（Θ、ｍ）を作成し、更に、凹関数の性質に基づき前記第１の補助関数Φ＋（Θ、ｍ）に対して更に第２の補助関数Φ＋〜（Θ、ｍ、α）及び補助変数αを作成することによって、その第２の補助関数を、前記コンピュータに前記補助関数として用いさせ、前記補助変数の更新を行う手順は、前記補助変数ｍとともに前記補助変数αも更新することを特徴とするプログラムである。
【００４５】
（２５）本発明は、上記課題を解決するために、コンピュータを、複数の音源からの信号を含む単独の入力音声信号中のそれぞれの音源を、正弦波重畳モデルのパラメータ推定の手法によってそれぞれ認識し、分離する音源分離装置として動作させるプログラムにおいて、前記入力音声信号を、所定フレームに区分けする手順と、前記区分けしたフレーム毎に、前記（１）〜（８）のいずれかに記載の正弦波重畳モデルパラメータ推定方法を用いて、そのフレームに含まれる所定のピッチの音声信号のパラメータをそれぞれの音源毎に推定する手順と、前記音源毎に前記フレーム毎に推定したパラメータをまとめて、各音源毎に音声信号を再構築することによって、各音源毎の信号をそれぞれ得る手順と、を、前記コンピュータに実行させることを特徴とするプログラムである。
【発明の効果】
【００４６】
本発明によれば、正弦波重畳モデルパラメータ推定を、補助関数を用いた繰り返し処理で実行できるので、パラメータをより正確に推定することができた。
【００４７】
また、その推定によって、複数の音源からの音を含む音声信号中から、音源の分離を行えば、１チャネルの音声信号からも音源分離することができる。
【発明を実施するための最良の形態】
【００４８】
以下、図面を用いて、発明の実施の形態を説明する。
【００４９】
第１  正弦波重畳モデルと目的関数
まず、本実施の形態で採用する正弦波重畳モデルと、目的関数とを説明する。
【００５０】
解析対象である混合信号を次のような正弦波重畳モデル
【数４７】
2008304718.tif 000040
でモデル化する。この中で未知パラメータはＫ個の調波信号の基本周波数
【数４８】
2008304718.tif 000041
ｎ次調波成分の複素振幅
【数４９】
2008304718.tif 000042
（すなわち、
【数５０】
2008304718.tif 000043
は振幅、
【数５１】
2008304718.tif 000044
は初期位相を表す。）であり、これらをまとめて下記式（６）のように表記する。
【００５１】
【数５２】
2008304718.tif 000045
ここで、考える問題は、観測混合信号ｙ（ｔ）のｔ＝０における最尤の瞬時特徴量
【数５３】
2008304718.tif 000046
を求めることであり、これは下記（７）式
【数５４】
2008304718.tif 000047
を最小化する
【数５５】
2008304718.tif 000048
に相当する。さて、窓関数
【数５６】
2008304718.tif 000049
は、先に用いた矩形窓の代わりに、下記式（８）の時刻ｔ＝０を中心としたＧａｕｓｓ窓（ガウス窓）
【数５７】
2008304718.tif 000050
とすると、
【数５８】
2008304718.tif 000051
は周波数領域ではＰａｒｓｅｖａｌの等式により、下記式（９）
【数５９】
2008304718.tif 000052
のように書くことができる。ただし、
【数６０】
2008304718.tif 000053
であり、ｄは窓の拡がりを表す定数であって、このｄは分析条件によって決まる定数である。また、
【数６１】
2008304718.tif 000054
は
【数６２】
2008304718.tif 000055
のＦｏｕｒｉｅｒ変換であり、すなわち、Ｇａｂｏｒ変換によるｔ＝０における時間周波数分解成分を意味する。さて、次章より、上記式（９）を最小化するパラメータ最適化アルゴリズムを導出していく。
【００５２】
なお、正弦波重畳モデルのフーリエ変換式
【数６３】
2008304718.tif 000056
において、ｄをｋ、ｎごとに自由度を持つパラメータ
【数６４】
2008304718.tif 000057
と見なして、
【数６５】
2008304718.tif 000058
と置き、
【数６６】
2008304718.tif 000059
を
【数６７】
2008304718.tif 000060
や
【数６８】
2008304718.tif 000061
と同様に推定対象としても良い。
【００５３】
第２．パラメータ推定アルゴリズム
２．１  補助関数を用いた目的関数の非増加更新法
本実施の形態で提案するパラメータ推定アルゴリズムは、ＥＭアルゴリズムからヒントを得て導かれた、補助関数を用いた目的関数の非増加更新法の考え方を基礎とする。なお、本アルゴリズムは、最適なパラメータが得られる可能性が高く、パラメータ最適化アルゴリズムと呼ぶことも好ましい。
【００５４】
次節以降の定式化に向けた準備として、まず補助関数の定義を行ったのちに、補助関数に関する補題を示す。
【００５５】
定義１．  パラメータ
【数６９】
2008304718.tif 000062
に関して最小化したい目的関数を
【数７０】
2008304718.tif 000063
とすると、下記式（１０）
【数７１】
2008304718.tif 000064
を満たすような
【数７２】
2008304718.tif 000065
を目的関数
【数７３】
2008304718.tif 000066
の補助関数と定義し、また、そのときの
【数７４】
2008304718.tif 000067
を補助変数と定義する。
【００５６】
補題１．  目的関数を
【数７５】
2008304718.tif 000068
とし、この目的関数の補助関数を
【数７６】
2008304718.tif 000069
とし、補助変数を
【数７７】
2008304718.tif 000070
とする。すると、この補助関数の補助変数に関する最小化（下記式（１１）と、パラメータ
【数７８】
2008304718.tif 000071
に関する最小化（下記式１２）
【数７９】
2008304718.tif 000072
を行うと、目的関数
【数８０】
2008304718.tif 000073
は単調減少する。この目的関数
【数８１】
2008304718.tif 000074
が下に有界ならば、上記の一連の更新の反復計算は収束性が保証される。
【００５７】
証明：  任意のパラメータを
【数８２】
2008304718.tif 000075
とし、上記式（１１）と式（１２）（数７９）により、
【数８３】
2008304718.tif 000076
が単調減少することを以下示す。まず、上記式（１１）により、下記式（１３）
【数８４】
2008304718.tif 000077
であり、上記式（１２）の更新により明らかに、下記式（１４）
【数８５】
2008304718.tif 000078
が成立する。補助関数の定義式（１０）に基づけば、下記式（１５）
【数８６】
2008304718.tif 000079
なので、
【数８７】
2008304718.tif 000080
である。以上より、各更新による目的関数の非増加が示された。（証明終わり）
この目的関数更新法をある最適化問題に応用するにあたり、鍵となるのは、以下の点である。
【００５８】
第一に、上記不等式（１０）の等号が成立するときの補助変数が解析的に求められるか否かという点、
第二に、補助関数を最小化するモデルパラメータが解析的に求められるかどうかという点、
これらの点を肯定的に満たすことができれば、最適化問題への応用の可能性が大である。換言すれば、そのようになるようにうまく補助関数を設計することが重要である。このような要件が、本実施の形態で提案する補助関数において特徴的なことであり、従来知られていなかったことである。ところで、この第２章の冒頭で述べたＥＭアルゴリズムは、この反復最適化法の原理の１つの例に相当する。
【００５９】
２．２  Ｆｅｄｅｒ  の補題（非特許文献２０参照）
上で述べた式（９）（数５９）の補助関数は、例えば、Ｆｅｄｅｒらが示唆した次のような補題によって作ることができる。
【００６０】
補題２．  
【数８８】
2008304718.tif 000081
においてある複素数値関数ｍｉ（ｘ）が下記式（１７）
【数８９】
2008304718.tif 000082
を満たすとき、
【数９０】
2008304718.tif 000083
が成り立つ。特にこの式は、下記式（１８）
【数９１】
2008304718.tif 000084
が成立するき、等号が成立する。ただし、
【数９２】
2008304718.tif 000085
は
【数９３】
2008304718.tif 000086
を満たす任意の係数である。
【００６１】
証明：  Ｌａｇｒａｎｇｅ  未定乗数法を用いた変分法により、上記（数９０）の不等式の右辺のｍｉ（ｘ）に関する最小値が左辺と等しくなることを示す。また、そのときのｍｉ（ｘ）を求める。
【００６２】
そこでまず、上記不等式（数９０）の右辺に条件（１７）のためのＬａｇｒａｎｇｅ  未定乗数項を付加した汎関数である下記式（１９）
【数９４】
2008304718.tif 000087
の極値を求めることとする。
【００６３】
【数９５】
2008304718.tif 000088
の
【数９６】
2008304718.tif 000089
に関する変分は、下記式（２０）
【数９７】
2008304718.tif 000090
と記述されるが、この式が恒等的に０であるためには、
【数９８】
2008304718.tif 000091
である必要がある。よって、下記式（２１）
【数９９】
2008304718.tif 000092
を０とおくと、下記式（２２）
【数１００】
2008304718.tif 000093
を得る。また、条件（１７）より、下記式（２３）
【数１０１】
2008304718.tif 000094
が得られる。これを式（２２）に代入すると、最終的に上記式（１８）を得る。これを上記不等式の右辺に代入すると、左辺と等しくなることを確認することができる。次に、この極値が最小解であるかどうかは、
【数１０２】
2008304718.tif 000095
の
【数１０３】
2008304718.tif 000096
に関するＨｅｓｓｅ行列
【数１０４】
2008304718.tif 000097
が正定値であることから示すことができる。（証明終わり）
２．３  補助関数の設計
さて、上記補題２より、
【数１０５】
2008304718.tif 000098
とし、
【数１０６】
2008304718.tif 000099
と置くと、上記式（９）に関して次のような不等式（２４）
【数１０７】
2008304718.tif 000100
が成立し、等号は
【数１０８】
2008304718.tif 000101
が
【数１０９】
2008304718.tif 000102
のとき成立する。ここで上記不等式（２４）の右辺を
【数１１０】
2008304718.tif 000103
と置く。定義１より、
【数１１１】
2008304718.tif 000104
は目的関数
【数１１２】
2008304718.tif 000105
の補助関数であり、
【数１１３】
2008304718.tif 000106
は補助変数である。そして式（２５）は、補題１の式（１１）（補助変数の更新式）に相当することは明らかである。
【００６４】
ところで、この不等式が意味するところは、混合信号と正弦波重畳モデルの複素スペクトル領域のＬ２ノルム誤差が、混合信号を分配フィルタ
【数１１４】
2008304718.tif 000107
によって任意に分解した各成分
【数１１５】
2008304718.tif 000108
と各正弦波成分とのＬ２ノルム誤差の総和をとったものの最小値となっているということである。そして、式（２５）は、混合信号と正弦波重畳モデルとの間にある誤差を配分率
【数１１６】
2008304718.tif 000109
に応じて｛ｋ，ｎ｝番目の正弦波成分に付加したものを
【数１１７】
2008304718.tif 000110
とすることが、
【数１１８】
2008304718.tif 000111
の最適な分解の仕方であることを表している。以後、表記の簡略化のため、
【数１１９】
2008304718.tif 000112
とする。補題１より、次に、
【数１２０】
2008304718.tif 000113
を
【数１２１】
2008304718.tif 000114
に関して最小化することを検討する。
【００６５】
【数１２２】
2008304718.tif 000115
は式（２６）
【数１２３】
2008304718.tif 000116
と書けるが、上式の第二項の積分はＧａｕｓｓ積分により解析的に計算可能であることが理解されよう。すなわち、下記式（２７）
【数１２４】
2008304718.tif 000117
であり、この項は
【数１２５】
2008304718.tif 000118
に依存しない。以上のことから、式（２６）は次のように書くことができる。
【００６６】
【数１２６】
2008304718.tif 000119
ここで、式（２８）の補助関数は依然として
【数１２７】
2008304718.tif 000120
に関して非線形であり、このままではまだ
【数１２８】
2008304718.tif 000121
の更新式に関しては解析的に求められない。そこで、
【数１２９】
2008304718.tif 000122
の更新式に関する検討の前に、
【数１３０】
2008304718.tif 000123
を固定した状態での振幅
【数１３１】
2008304718.tif 000124
と位相
【数１３２】
2008304718.tif 000125
の更新式に関しては解析的に得られることをまず確認しておこう。
【００６７】
まず、
【数１３３】
2008304718.tif 000126
を解くと、下記式（２９）
【数１３４】
2008304718.tif 000127
を得る。振幅と位相は、それぞれ上で求まった
【数１３５】
2008304718.tif 000128
を用いて
【数１３６】
2008304718.tif 000129
と書ける。なお、２階偏導関数が正であることから、この極値が最小解であることが明らかである。次節では、指数関数の性質に着目し、それに基づき、
【数１３７】
2008304718.tif 000130
の更新式を解析的に得るための他の補助関数が作れることを示す。
【００６８】
２．４  凹関数の性質に基づく補助関数の設計
上記式（２８）において
【数１３８】
2008304718.tif 000131
の更新式が解析的に得られない原因となっているのは、非線形項
【数１３９】
2008304718.tif 000132
である。ここで、本願発明者らが着目したのは、−ｅ−ｘのような型の関数が微分可能な凹
関数である点と、微分可能な凹関数に対して次のような定理が成り立つということである。
【００６９】
補題３．  実関数ｆ（ｘ）が微分可能な凹関数のとき、任意の点
【数１４０】
2008304718.tif 000133
において次の不等式（３０）が成り立つ。
【００７０】
【数１４１】
2008304718.tif 000134
証明：凹関数の定義より、任意の異なる２点
【数１４２】
2008304718.tif 000135
及び、任意の実数
【数１４３】
2008304718.tif 000136
に対して、下記（３１）式
【数１４４】
2008304718.tif 000137
が成り立つとき、ｆ（ｘ）を凹関数と言う。上式は、下記（３２）式
【数１４５】
2008304718.tif 000138
と書き換えることができるが、
【数１４６】
2008304718.tif 000139
のもとでは、ｆ（ｘ）は微分可能であるので、下記（３３）式
【数１４７】
2008304718.tif 000140
が成り立つ。したがって、この（３３）式を式（３２）に代入すると、下記（３４）式
【数１４８】
2008304718.tif 000141
を得る。（証明終わり）
さて、−ｅ−ｘは微分可能な凹関数であるから、補題３より、任意の点
【数１４９】
2008304718.tif 000142
において下記式（３５）
【数１５０】
2008304718.tif 000143
が成り立つ。したがって、式（３５）のｘを
【数１５１】
2008304718.tif 000144
とおき、
【数１５２】
2008304718.tif 000145
を実関数
【数１５３】
2008304718.tif 000146
と対応させると、
【数１５４】
2008304718.tif 000147
において下記
【数１５５】
2008304718.tif 000148
のような不等式（３６）が立てられる。この式（３６）  より、
【数１５６】
2008304718.tif 000149
のように、式（２８）の上限を得る。この不等式右辺を
【数１５７】
2008304718.tif 000150
と置くと、下記式（３７）
【数１５８】
2008304718.tif 000151
であるので、定義１より、これもまた
【数１５９】
2008304718.tif 000152
の補助関数として扱うことができる。また、
【数１６０】
2008304718.tif 000153
及び
【数１６１】
2008304718.tif 000154
が、この場合、いずれも補助変数となる。そして（数１５６）の等号成立は、ｍが式（２５）の場合、及び、下記式（３８）
【数１６２】
2008304718.tif 000155
の場合である。
【００７１】
本実施の形態において特徴的なことは、以上説明したような補助関数を導いたことである。
【００７２】
特に、本実施の形態では、ガウス窓を窓関数として利用したが、このガウス窓関数は指数関数で表されている。そして、この指数関数が凹関数である点に着目し、凹関数について成り立つ上記不等式を利用して補助関数を作成した点が、本実施の形態において特徴的な事項であり、従来行われていなかった事項である。本実施の形態では、このようにして補助関数を導いたことによって、
【数１６３】
2008304718.tif 000156
を最小化する
【数１６４】
2008304718.tif 000157
を解析的に求めることができたのである。この点が、従来の技術と顕著に相違する点である。
【００７３】
そのため、このガウス窓関数を好適に適用することができる信号であれば、本実施の形態の正弦波重畳モデルパラメータの推定手法を適用することが好ましい。
【００７４】
ここでは実際にピッチ周波数
【数１６５】
2008304718.tif 000158
の更新式が解析的に得られることを確認的に示しておく。まず、
【数１６６】
2008304718.tif 000159
を解くと、
【数１６７】
2008304718.tif 000160
であるならば、下記式（３９）
【数１６８】
2008304718.tif 000161
が得られる。
【００７５】
さらに、ｄを自由パラメータ
【数１６９】
2008304718.tif 000162
と置いた場合、この補助関数を最小化する
【数１７０】
2008304718.tif 000163
は、やはり、
【数１７１】
2008304718.tif 000164
のように、解析的に得ることができる。
【００７６】
ｄは、操作者が分析条件によって決めるものであるので、
【数１７２】
2008304718.tif 000165
の真値は既知である（
【数１７３】
2008304718.tif 000166
を求めること自体には意味がない）が
【数１７４】
2008304718.tif 000167
を小さい値として初期値設定して
【数１７５】
2008304718.tif 000168
や
【数１７６】
2008304718.tif 000169
とともに更新していくことで
【数１７７】
2008304718.tif 000170
や
【数１７８】
2008304718.tif 000171
が局所解に陥るのを回避する効果がもたらされる。
【００７７】
このｄは、もともとは、周期信号モデルにかかるガウス窓関数の広がりを表すパラメータである。従って、ｄを十分小さくとるということは、周期信号のｔ＝０近辺だけを抜き出し、裾を押さえるという意味がある。
【００７８】
周期信号とは文字通り、周期的な関数であるが、ｄを小さくとることによって、観測信号は周期的でも、モデルの方の裾の周期性を小さく抑えることができる。
【００７９】
この状態では、
【数１７９】
2008304718.tif 000172
と
【数１８０】
2008304718.tif 000173
とのラフな解が見つけられることができると考えられる。そして、ｄを更新していくということは、徐々に周期信号の裾の方まで正確にフィットさせていくことに相当するので、最終的に
【数１８１】
2008304718.tif 000174
と
【数１８２】
2008304718.tif 000175
とがより正確に求まると考えられる。
【００８０】
なお、ｄを変数とする手法の他にも、
【数１８３】
2008304718.tif 000176
と
【数１８４】
2008304718.tif 000177
の反復推定中に、ｄをある一定のスピードで小さい値から徐々に大きくしていき最終的に真値にする、というやり方も効果が期待できる。
【００８１】
２．５  本実施の形態におけるパラメータ推定アルゴリズムのまとめ
以上述べたように、本実施の形態の正弦波重畳モデルのパラメータ推定アルゴリズムは次のようにまとめられる。また、この反復計算における各ステップでの目的関数値の変遷の様子が図１に示されている。
【００８２】
この図１に示すように、ｍの更新を行い、また、αの更新を行って補助関数を小さくすれば、そのときのパラメータの値は、目的関数も小さくさせる。
【００８３】
さて、本実施の形態におけるパラメータ推定アルゴリズムは、以下のようなステップから構成される。
【００８４】
Ｓｔｅｐ０：
【数１８５】
2008304718.tif 000178
の初期設定。
【００８５】
Ｓｔｅｐ１：上述した式（２５）によって
【数１８６】
2008304718.tif 000179
を更新する。
【００８６】
Ｓｔｅｐ２：上述した式（３８）によって
【数１８７】
2008304718.tif 000180
を更新する。
【００８７】
Ｓｔｅｐ３：上述した式（２９）によって
【数１８８】
2008304718.tif 000181
を更新する。
【００８８】
Ｓｔｅｐ４：上述した式（３９）によって
【数１８９】
2008304718.tif 000182
を更新する。その後、上記Ｓｔｅｐ１に戻る。
【００８９】
２．６  補助関数の要件
さて、これまで、本実施の形態で提案する補助関数の要件を述べてきた。
【００９０】
（１）目的関数＝＜補助関数
（２）上記式の等号成立するときの補助変数が解析的に得られること
（３）補助関数を最小化する（求めたい）パラメータが解析的に得られること
これらの要件を満たす補助関数をこれまで説明してきたが、特に、上で説明した補助関数の特徴の中で、上記３要件のうち（３）を満たしている要因は、補助関数が推定したいパラメータのベクトルの２次式で表されている点である。このような形で、且つ、上記３要件のうち（１）と（２）を満たす補助関数であれば、上で説明した補助関数以外の関数でも、上記３要件を満たすことができよう。
【００９１】
２．７  補助関数の形式
なお、請求の範囲では、補助関数の表現形式として、
【数１９０】
2008304718.tif 000183
という表現形式を示している。本実施の形態で提案する補助関数はこの式の好適な一態様の一つに相当することを示そう。
すなわち、
【数１９１】
2008304718.tif 000184
とおき、
【数１９２】
2008304718.tif 000185
とおき、
【数１９３】
2008304718.tif 000186
とおき、
【数１９４】
2008304718.tif 000187
とおけば、本実施の形態の補助関数が、上記の表現開式になることが理解されよう。
３．評価実施例
３．１  シミュレーションによる収束性の確認
、上述した正弦波重畳モデルパラメータ推定方法を音源分離に応用した場合の具体例について結果を紹介する。この実施例によって、本実施の形態で提案した方法と、従来の勾配法によるパラメータ推定の初期値依存性や収束速度の比較結果を説明する。
【００９２】
使用した従来の勾配法
ここで比較対象の勾配法を用いたパラメータ推定法（以後、単に勾配法と呼ぶ。）は、Ｊｉｎａｃｈｉｔｒａ  の手法（非特許文献１２）を参考にしたもので、式（１１）により
【数１９５】
2008304718.tif 000188
を更新するステップ、式（２９）によって
【数１９６】
2008304718.tif 000189
を更新するステップ、式（２８）を
【数１９７】
2008304718.tif 000190
に関する最急降下更新により減少させるステップ、の各ステップによって構成される。
【００９３】
この従来手法と本実施の形態で提案する手法を比べることによって、本特許で提案したピッチ周波数推定方式が局所解回避能力などの面で有効かどうかを示すことができる。この比較実験では、パラメータの真値があらかじめ分かっている信号（合成信号）を分析対象とした。具体的には、２つの周期信号（ピッチ周波数が２７０Ｈｚ、２００Ｈｚ）は１０個の調波成分からなり、各成分の振幅と位相は乱数生成により決め、これらを加算して混合信号を作成した。ｎ次調波成分の振幅と位相の乱数値の範囲はそれぞれ
【数１９８】
2008304718.tif 000191
及び
【数１９９】
2008304718.tif 000192
とした。また、正弦波重畳モデルにおいてＫ＝２，Ｎ＝１０とした。この合成信号（１６ｋＨｚ  サンプリング）  からＧａｂｏｒ変換により短時間複素スペクトル
【数２００】
2008304718.tif 000193
を得た。ただし、Ｇａｂｏｒ変換における拡散パラメータｄは０．０６７とした。
【００９４】
様々な初期値条件のもとで、ピッチ周波数
【数２０１】
2008304718.tif 000194
が、本実施の形態で提案する方法及び従来の勾配法によって各ステップで更新されていく様子をそれぞれ図２、図３に示す（ピッチ周波数以外のパラメータの更新の様子については省略する）。図２及び図３のグラフは、横軸は、繰り返しの数であり、縦軸は、想定しているピッチである。繰り返し回数が０の場合は、それぞれのモデルは初期値のパラメータを備えている。初期の段階は、種々のピッチ（周波数）を初期値と考えている。
【００９５】
これらグラフは、ピッチの
【数２０２】
2008304718.tif 000195
の更新値の変遷を表している。なお、振幅
【数２０３】
2008304718.tif 000196
の初期値は０とした。図２、図３を比べてみると、勾配法では、
【数２０４】
2008304718.tif 000197
の初期値が真値（２７０Ｈｚ、２７０Ｈｚ）にある程度近い値でないと、真値以外の停留点に陥ってしまうことが多いが、本実施の形態で提案する方法では、広範囲の初期値から真値に素早く収束していく様子が確認できる。このシミュレーション実験の結果は、勾配法を原理として用いる従来から知られている多くの従来手法よりも、本実施の形態で提案する方法が局所解回避能力と収束速度の面で優れていることを示す１つの例証と言えよう。
【００９６】
図２のグラフ（本実施の形態）は、初期値がどの値でも、繰り返しを３０回−４０回程度繰り返せば、概ね、想定していた値、すなわち２７０Ｈｚと、２００Ｈｚに収束していくことが理解されよう。図２（１）は、２７０Ｈｚを推定した場合の各モデルのパラメータ（ピッチ周波数）の変化の様子を示すグラフである。図２（２）は、２００Ｈｚと推定した場合の各モデルのパラメータ（ピッチ周波数）の変化の様子を示すグラフである。このように、本実施の形態で提案する手法にでは、様々な初期値から、速やかに想定しているピッチへ収束している。
【００９７】
図３のグラフ（従来の勾配法）では、初期値の値によっては、収束しなかったり、誤ったピッチ周波数に収束してしまう現象が見られた。すなわち２７０Ｈｚや、２００Ｈｚに近い初期値から始めた場合は収束していっているが、離れた周波数からシミュレーションを行った場合は、倍音のピッチに収束してしまったりする現象が見られた。図３（１）は、２７０Ｈｚを推定した場合の各モデルのパラメータ（ピッチ周波数）の変化の様子を示すグラフである。図３（２）は、２００Ｈｚと推定した場合の各モデルのパラメータ（ピッチ周波数）の変化の様子を示すグラフである。
【００９８】
３．２  混合音声の１ｃｈブラインド音源分離
次に、本実施の形態で提案する正弦波重畳モデルのパラメータ推定法を、１ｃｈブラインド音源分離方法を実現した。そして、その結果から、正弦波重畳モデルパラメータ推定法並びに１ｃｈブラインド音源分離方法の基本性能を確認した。
【００９９】
ここで対象とする混合信号は、ＡＴＲ研究用日本語音声データベースのセットＢ中の音声データを用い、男性同士の音声信号を加算したもの、女性同士の音声信号を加算したもの、男性と女性の音声波形を加算したもの、の３種類に混合信号の分離を試みた。
【０１００】
すべての音声データは１６ｋＨｚでリサンプリングし、フレーム間隔１０ｍｓのＧａｂｏｒ変換により周波数解析した。前節同様、Ｇａｂｏｒ変換における拡散パラメータｄは０．０６７とした。
【０１０１】
なお、フレームは一般にオーバーラップしており、フレーム間隔が１０ｍｓで、フレームの長さが３０ｍｓという場合もある。この場合は、時刻０ｍｓ−３０ｍｓが第１フレーム、時刻１０ｍｓ−４０ｍｓが第２フレーム、時刻２０ｍｓ−５０ｍｓが第３フレームの
ように設定される。
【０１０２】
このように、音声信号は１チャネルのみで、１０ｍｓｅｃのフレームに分けられており、各フレーム毎に、正弦波重畳モデルパラメータの推定を行う。
【０１０３】
また、正弦波重畳モデルの各調波信号の調波成分数Ｎは３０とした。すなわち、第３０次調波信号まで取り扱うことにした。
【０１０４】
ここで採用したアルゴリズムの内容は次の通りである。
【０１０５】
まず、正弦波重畳モデルの初期調波信号数Ｋを１０から開始した。これは、１フレーム中に１０種類のピッチの信号が含まれることを念頭に置いたものである。しかし、パラメータの反復推定の過程で、複数の調波信号モデルのピッチ周波数パラメータが
（１）ある一定値より近くなった場合、
（２）ほぼ整数倍関係になった場合、
には低い方の調波信号だけを残し、残りは除外していくというアルゴリズムを採用した。その結果、一般的にはＫの数は小さくなっていく。そして、収束した後、エネルギー全パワーの大きい２つの調波信号だけを残し（つまりＫ＝２）、再びパラメータ推定を行う。このようなアルゴリズムを採用したのは、あらかじめ２種の信号の混合信号であることを前提としたシミュレーションだからである。
【０１０６】
このようにして、最終的に得られた２つの調波信号を分離推定信号とする。また、
【数２０５】
2008304718.tif 000198
の初期値は、まず観測混合信号の複素スペクトルの実部又は虚部において極大／極小点を与える周波数をすべて見つけ、その中でパワーの大きいものから順に選んだ１０個の周波数とした。上記は、各短時間窓（フレーム）で分離信号を推定する手続きであるが、各々のフレームで分離した信号がどの音源に対応したものであるかは本来は定かではない。  今回の実験では、この対応づけができている状況での仮の音源分離性能を調べるために、分離信号がどの音源に対応するかを混合する前の各信号との近さを調べて判定することにした。
【０１０７】
以上の条件のもとで、図４に示す混合音声を実際に分離した結果の例を図５に示す。
【０１０８】
図４・図５は、横軸が時間を表し、縦軸が音声信号の振幅を表すグラフである。まず、図４（１）は、女性話者の音声であり、図４（２）は男性話者の音声である。図４（３）はこれらの混合信号である。図５は、図４（３）の混合信号を分離した結果であり、図５（１）は、分離した女性話者の音声信号、図５（２）は分離した男性話者の音声信号である。
【０１０９】
男性話者Ａと女性話者Ｂの混合音声（男性話者Ａから見てＳＮ比が−０．３ｄＢ）の分離後の双方のＳＮ比は、７．２ｄＢ，６．４ｄＢ（７．５ｄＢ，６．１ｄＢの改善）、女性話者Ａと女性話者Ｂの混合音声（女性話者Ａから見てＳＮ比が１．５ｄＢ）の分離後の双方のＳＮ  比は６．０ｄＢ，  ４．８ｄＢ（４．５ｄＢ，  ６．３ｄＢ  の改善）、男性話者Ａと男性話者Ｂの混合音声（男性話者Ａから見てＳＮ比が−０．３ｄＢ）の分離後の双方のＳＮ  比は４．８ｄＢ，４．３ｄＢ（５．１ｄＢ，４．０ｄＢの改善）であった。
【０１１０】
本実施の形態で提案する方法では、話者間のピッチ周波数の違いが、音源分離の手がかりとなっているため、同性同士の混合音声の分離精度が異性同士のそれよりいくらか下回ってしまうのは概ね予想通りの結果と言えよう。
【０１１１】
本実施の形態で提案した手法は、フレーム毎に独立にパラメータ推定をしているため、場合によっては、分離信号の位相変化が不連続になったり、振幅が急激に変化する可能性がある。そこで、変化が急峻とならないように、なめらかに変化させたパラメータを用いて原音源を推定することも好適である。
【０１１２】
第３．まとめ
本特許では、正弦波重畳モデルのパラメータ推定問題の核である周波数推定ないし基本周波数推定の難しさの本質が、正弦波重畳モデルが周波数パラメータに関して非線形である点にあることに注目した。そして、本願発明者らは、ＥＭアルゴリズムからヒントを得て、補助関数を用いた新しい反復推定アルゴリズムを導いた。
【０１１３】
更に、この反復推定アルゴリズムを用いて、フレーム毎に異なる音源からの音声のパラメータを推定し、音源毎の音声を得る（音源分離）作業を行い、本アルゴリムの有効性を確認した。
【０１１４】
特に、従来の音源分離では、複数チャネルの音声信号を用いて、各音源の方向を推定してその推定方向から各音源の音声信号を分離・推定する場合が多かったが、本実施の形態にかかる手法では、各時間フレーム毎に正弦波重畳モデルのパラメータ推定を行って、ピッチの相違から各音源の信号を分離している。したがって、１チャネルの音声信号のみを用いつつ、実際に有用な音源分離を実現することができた。
【０１１５】
第４．  その他
（１）上記実施の形態では、正弦波重畳モデルパラメータの推定方法について述べたが、同方法はコンピュータによって実行することが好ましい。したがって、上記処理をコンピュータに実行させるプログラムもまた本発明の１種である。同様にこのようなプログラムを所定の記録媒体に格納しておけば、その記録媒体からコンピュータにプログラムを実行させ、上記方法を実行させることができる。また、この記録媒体からプログラムをコンピュータにインストールすれば、上記方法を実行する装置が得られる。この装置（すなわち正弦波重畳モデルパラメータ推定装置）もまた本発明の１種である。そのような装置においては、上記方法の各ステップ・各処理を実行する手段が、プログラムで実現されている。
【０１１６】
（２）上記実施の形態では、正弦波重畳モデルパラメータの推定方法について述べ、音声信号に関して、１ｃｈの信号に対する音源分離のアプリケーション例を説明した。
【０１１７】
しかし、上記正弦波重畳モデルパラメータの推定方法は、正弦波重畳モデルが適用できれば、どのような「波」にも利用できることは言うまでもない。例えば、地震波の解析、超音波検査装置や超音波診断装置、各種ＣＴ装置等、の信号処理にも応用することが可能である。すなわち、地震波、振動解析、超音波探傷、電気通信における波形解析・診断、医療診断解析等の分野における信号処理に応用可能である。
【図面の簡単な説明】
【０１１８】
【図１】本実施の形態の正弦波重畳モデルのパラメータ推定アルゴリズムの反復計算における各ステップでの目的関数値の変遷の様子を示す説明図である。
【図２】ピッチ周波数が、本実施の形態で提案する方法によって各ステップで更新されていく様子を示す説明図である。
【図３】ピッチ周波数が、従来の勾配法によって各ステップで更新されていく様子を示す説明図である。
【図４】混合音声の波形の例を示すグラフである。
【図５】図４の混合音声を音源分離した結果を示すグラフである。
【符号の説明】
【０１１９】
μ１、μ２  音声信号に含まれる所定の信号のピッチ
【図１】
2008304718.tif 000208
【図２】
2008304718.tif 000209
【図３】
2008304718.tif 000210
【図４】
2008304718.tif 000211
【図５】
2008304718.tif 000212