(57)【要約】【課題】振幅包絡が変化せず、因果性と完全再構成性を充たすモード分解ができるモード分解装置を提供する。【解決手段】モード分解装置は、楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、楽音信号を各モードに分解する。【選択図】図１
2019159018.tif 000002
【特許請求の範囲】
【請求項１】
楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、前記楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、前記楽音信号を各モードに分解するモード分解装置。
【請求項２】
請求項１に記載のモード分解装置であって、
楽音信号sまたはその近似音から極を取り出す極取り出し部と、
iを前記楽音信号sの各モードのインデックスとし、前記極の情報に基づいて定まるi番目の共振周波数fiの共振フィルタの振幅応答|hi|に基づいて、i番目の第１の重みwi(1)を推定する動作を各モードについて繰り返し実行する第１重み推定部と、
jをj≠iを満たす楽音信号sのモードのインデックスとし、j番目の共振周波数fjの共振フィルタの振幅応答|hj|に基づいて、i番目の第２の重みwi(2)を推定する動作を各モードについて繰り返し実行する第２重み推定部と、
残差rに対する重みである残差重みwrを推定する残差重み推定部と、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積であるi番目の第３の重みwiと、i番目の前記モードxiの周波数成分F(xi)と、前記楽音信号sの周波数成分F(s)と、パラメータλと、前記残差重みwrと、前記残差rの周波数成分F(r)に関する最小二乗問題
【数３０】
2019159018.tif 000033
の解を、完全再構成性と因果性を充たすことを制約として求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解部
を含むモード分解装置。
【請求項３】
請求項２に記載のモード分解装置であって、
前記第１重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、μをパラメータとし、i番目の第１の重みwi(1)[k]を、
【数３１】
2019159018.tif 000034
と推定する
モード分解装置。
【請求項４】
請求項２または３に記載のモード分解装置であって、
前記第２重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、i番目の第２の重みwi(2)[k]を
【数３２】
2019159018.tif 000035
と推定する
モード分解装置。
【請求項５】
請求項２から４の何れかに記載のモード分解装置であって、
前記残差重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、残差重みwr[k]を
【数３３】
2019159018.tif 000036
と推定する
モード分解装置。
【請求項６】
モード分解装置が実行するモード分解方法であって、
楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、前記楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、前記楽音信号を各モードに分解するモード分解方法。
【請求項７】
請求項６に記載のモード分解方法であって、
楽音信号sまたはその近似音から極を取り出す極取り出しステップと、
iを前記楽音信号sの各モードのインデックスとし、前記極の情報に基づいて定まるi番目の共振周波数fiの共振フィルタの振幅応答|hi|に基づいて、i番目の第１の重みwi(1)を推定する動作を各モードについて繰り返し実行する第１重み推定ステップと、
jをj≠iを満たす楽音信号sのモードのインデックスとし、j番目の共振周波数fjの共振フィルタの振幅応答|hj|に基づいて、i番目の第２の重みwi(2)を推定する動作を各モードについて繰り返し実行する第２重み推定ステップと、
残差rに対する重みである残差重みwrを推定する残差重み推定ステップと、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積であるi番目の第３の重みwiと、i番目の前記モードxiの周波数成分F(xi)と、前記楽音信号sの周波数成分F(s)と、パラメータλと、前記残差重みwrと、前記残差rの周波数成分F(r)に関する最小二乗問題
【数３４】
2019159018.tif 000037
の解を、完全再構成性と因果性を充たすことを制約として求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解ステップ
を含むモード分解方法。
【請求項８】
請求項７に記載のモード分解方法であって、
前記第１重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、μをパラメータとし、i番目の第１の重みwi(1)[k]を、
【数３５】
2019159018.tif 000038
と推定する
モード分解方法。
【請求項９】
請求項７または８に記載のモード分解方法であって、
前記第２重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、i番目の第２の重みwi(2)[k]を
【数３６】
2019159018.tif 000039
と推定する
モード分解方法。
【請求項１０】
請求項７から９の何れかに記載のモード分解方法であって、
前記残差重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、残差重みwr[k]を
【数３７】
2019159018.tif 000040
と推定する
モード分解方法。
【請求項１１】
コンピュータを請求項１から５の何れかに記載のモード分解装置として機能させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は楽音信号を複数のモードに分解するモード分解装置、モード分解方法、プログラムに関する。
【背景技術】
【０００２】
楽器音のモードの減衰特性は楽器の種類や音の高さ、強さによって異なり、音色に大きな影響を与える。そこで、楽器音に対してモードの減衰特性の解析が盛んに行われている（非特許文献１）。特に打楽器や撥弦楽器の各モードの減衰特性において、発音時の挙動は重要である。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】T. Cheng, S. Dixon and M. Mauch, "Modelling the decay of piano sounds," Int. Conf. Acoust., Speech Signal Process. (ICASSP), pp. 594{598, Apr. 2015.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかし、発音時に非調波成分が存在するため、各モードの発音時の挙動の解析は容易ではない。楽器音sをN個の中心周波数fi(i=1,...,N)の狭帯域なモードmiと非調波成分等の残差rに分解することを考える。このようなモード分解では一般にバンドパスフィルタが用いられる。しかし、フィルタの特性により各モードの振幅包絡が変化する他に、位相遅延やプリリンギングが発生する。因果的なフィルタを用いると、各モードに位相遅延が発生し、モードmi及び残差rの足し合わせと元の楽器音sが一致する完全再構成性が満たされない。一方、位相遅延のないフィルタを用いるとプリリンギングが発生し、因果性が満たされない。バンドパスフィルタ以外の手法として、短時間フーリエ変換(STFT, Short-Time Fourier Transform)や経験的モード分解(EMD, Empirical Mode Decomposition)を用いる手法も提案されているが、下表にまとめる通り、上記の課題をすべて解決するモード分解にはなっていない。
【０００５】
【表１】
2019159018.tif 000003
【０００６】
そこで本発明では、振幅包絡が変化せず、因果性と完全再構成性を充たすモード分解ができるモード分解装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明のモード分解装置は、楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、楽音信号を各モードに分解する。
【発明の効果】
【０００８】
本発明のモード分解装置によれば、振幅包絡が変化せず、因果性と完全再構成性を充たすモード分解ができる。
【図面の簡単な説明】
【０００９】
【図１】実施例１のモード分解装置の構成を示すブロック図。
【図２】実施例１のモード分解装置の動作を示すフローチャート。
【図３】設計した共振フィルタ|h5|の振幅応答（上段左）、第１重みw5(1)（上段右）、第２重みw5(2)と残差重みwr（下段左）、第３重みw5（下段右）を示す図。
【図４】フィルタと提案手法によるモード分解の結果を示す図。
【図５】実施例２のモード分解装置の構成を示すブロック図。
【図６】実施例２のモード分解装置の動作を示すフローチャート。
【図７】実施例３のモード分解装置の構成を示すブロック図。
【図８】実施例３のモード分解装置の動作を示すフローチャート。
【発明を実施するための形態】
【００１０】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【実施例１】
【００１１】
＜理論＞
以下、実施例１のモード分解装置１の基礎となる理論を説明する。因果的な線形フィルタを用いると、位相遅延が起こる。一方。このフィルタと同じ振幅応答|hi|を持つ位相遅延のないフィルタを考えると、フィルタ処理は
【００１２】
【数１】
2019159018.tif 000004
【００１３】
と表せるが、プリリンギングが発生する。ただし、F，F-1，○はそれぞれ離散フーリエ変換、逆離散フーリエ変換、アダマール積を表す。そこで、式（１）で示したフィルタ処理を最小二乗法を用いて再定式化することで、フィルタ処理を拡張する。フィルタの振幅応答|hi|に逆特性|hi|-1が存在する時、式（１）を変形すると
【００１４】
【数２】
2019159018.tif 000005
【００１５】
となる。式（２）を最小二乗問題に置き換えると
【００１６】
【数３】
2019159018.tif 000006
【００１７】
と表せる。この問題の解は式（１）と一致するため、フィルタ処理は最小二乗問題の解として解釈できる。このように、最小二乗問題として定式化することで、制約や他のコスト関数の導入が可能である。
【００１８】
最小二乗問題（式（３））は各モードごとのフィルタ処理に対応し、完全再構成性、因果性は考慮されていない。そこで、これらを制約として、フィルタ処理に対応する最小二乗問題を解くことで、上記の性質を満たしたモード分解を実現する。すなわち、以下の制約付き最小二乗問題
【００１９】
【数４】
2019159018.tif 000007
【００２０】
を解くことで、楽器音sをN個のモードxi及び残差rへ分解する。ただし、nは時間サンプルを、τは楽器の発音時刻に対応するサンプルを表す。また、λは残差の大きさを調整するパラメータである。この問題は凸なので、適当な凸最適化アルゴリズムで解くことができる。式（４）では、全てのモード及び残差を同時に扱い、下段左式の制約を導入することで、完全再構成性を考慮することができる。また、下段右式の制約は各モードxiの因果性を表す。これにより、得られた各モードxi及び残差rが上記の性質を同時に満たすモード分解が可能である。式（４）上段のコスト関数における周波数領域での重みwiを|hi|-1の拡張として考えると、コスト関数の最小化はフィルタ処理の一般化に対応する。
【００２１】
式（４）における重みwiは各モードxiが所望の周波数特性を持つように任意に作成することができる。楽器音のモード分解を考えるため、各モードxiが狭帯域かつ他のモードの成分の混入が少ない信号にする重みwiを作成する。ARモデルを用いて楽器音sを近似し、極の情報から共振周波数fiの共振フィルタを設計した。共振フィルタの振幅応答を|hi|とし、楽器音sのスペクトルのピークでモードmiに対応するピークの大きさをgiとする。これらの情報を利用し、各モードxiを中心周波数fiの狭帯域な信号にする第１の重みwi(1)と、他のモードの混入を防ぐ第２の重みwi(2)を作成する。これら2つの重みのアダマール積wi(1)○wi(2)を式（４）の重みwi（第３の重み）とした（詳細は実施例１において後述）。
【００２２】
＜理論に基づく本発明の装置の一般的表現＞
すなわち、以下の実施例に記載するモード分解装置１は、楽音信号sを各モードm1,...,mN（あるいはx1,...,xN）に分解するフィルタの振幅応答|hi|の逆特性|hi|-1（第３の重みwiとして拡張）と各モードm1,...,mN（あるいはx1,...,xN）の周波数成分F(mi)（あるいはF(xi)）のアダマール積と、楽音信号sの周波数成分F(s)に関する最小二乗問題であって、完全再構成性（（４）下段左式）と因果性（（４）下段右式）を充たすことを制約とした制約付き最小二乗問題を解くことにより、楽音信号を各モードに分解する装置であるといえる。
【００２３】
＜モード分解装置１＞
以下、図１を参照して実施例１のモード分解装置１の構成を説明する。同図に示すように、本実施例のモード分解装置１は、極取り出し部１１と、第１重み推定部１２と、第２重み推定部１３と、残差重み推定部１４と、モード分解部１５を含む構成である。
【００２４】
以下、図２を参照して各部の動作について説明する。まず、極取り出し部１１は、楽音信号sまたはその近似音から極を取り出す（Ｓ１１）。
【００２５】
第１重み推定部１２は、iを楽音信号sの各モードのインデックスとし、極の情報に基づいて定まるi番目の共振周波数fiの共振フィルタの振幅応答|hi|に基づいて、i番目の第１の重みwi(1)を推定する動作を各モードについて繰り返し実行する（Ｓ１２）。具体的には、第１重み推定部１２は、kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、μを上限を表すパラメータとし、i番目の第１の重みwi(1)[k]を、
【００２６】
【数５】
2019159018.tif 000008
【００２７】
と推定する（Ｓ１２）。
【００２８】
第２重み推定部１３は、jをj≠iを満たす楽音信号sのモードのインデックスとし、j番目の共振周波数fjの共振フィルタの振幅応答|hj|に基づいて、i番目の第２の重みwi(2)を推定する動作を各モードについて繰り返し実行する（Ｓ１３）。具体的には、第２重み推定部１３は、kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、i番目の第２の重みwi(2)[k]を
【００２９】
【数６】
2019159018.tif 000009
【００３０】
と推定する（Ｓ１３）。
【００３１】
残差重み推定部１４は、残差rに対する重みである残差重みwrを推定する（Ｓ１４）。具体的には、残差重み推定部１４は、kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、残差重みwr[k]を
【００３２】
【数７】
2019159018.tif 000010
【００３３】
と推定する（Ｓ１４）。
【００３４】
モード分解部１５は、i番目の第１の重みwi(1)とi番目の第２の重みwi(2)のアダマール積であるi番目の第３の重みwiと、i番目のモードxiの周波数成分F(xi)と、楽音信号sの周波数成分F(s)と、パラメータλと、残差重みwrと、残差rの周波数成分F(r)に関する最小二乗問題
【００３５】
【数８】
2019159018.tif 000011
【００３６】
の解を、完全再構成性と因果性を充たすこと、すなわち
【００３７】
【数９】
2019159018.tif 000012
【００３８】
を制約として求めることにより、楽音信号sを各モードx1,...,xNに分解する（Ｓ１５）。
【００３９】
＜実験＞
一例として、The University of Iowa Musical Instrument Samples (MIS)（参考非特許文献１）のbellsC5stereobrassを提案手法により解析した。ARモデルの次数を300とし、楽器音sの近似にはBurg法（参考非特許文献２）を用いた。また、重みwi(1)のパラメータをμ=10-3とした。制約付き最小二乗問題（式（４））は、λ=3×10-3とし、Over-relaxation Monotone Fast Iterative Shrinkage-Thresholding Algorithm (OMFISTA)（参考非特許文献３）により解いた。第５モードに対応する共振フィルタの振幅応答|h5|,第１重みw5(1),第２重みw5(2)，第３重みwi，残差重みwrを図３に示す。
【００４０】
提案手法で得られた各モード及び残差の波形x1，x2，x3，x4，rと、式（１）で表される位相遅延のない共振フィルタを用いた場合の各モード及び残差の波形を図４に示す。共振フィルタでは，フィルタの特性により各モードの音の立ち上がりが鈍くなり、プリリンギングも発生している。一方、提案手法ではプリリンギングが発生せず、残差が発音時の非調波成分に集中している。また、OMFISTAを用いることで式（４）下段の制約は常に満たされるため、完全再構成性も必ず満たされる。
（参考非特許文献１：Electronic Music Studios School of Music University of Iowa，"the University of Iowa Electronic Music Studios," Samples, http://theremin.music.uiowa.edu/MIS.html）
（参考非特許文献２：S. Kay, Modern spectral estimation: Theory and application, Prentice Hall, New Jersey, 1988.）
（参考非特許文献３：M. Yamagishi and I. Yamada, "Over-relaxation of the fast iterative shrinkage-thresholding algorithm with variable stepsize" IOP Sci. Inv. Probl., vol.27, no. 10, p. 105008, Sep. 2011．）
【実施例２】
【００４１】
＜理論＞
以下、実施例２のモード分解装置２の基礎となる理論を説明する。最小二乗問題（式（３））は各モードごとの線形なフィルタ処理に対応しているが、完全再構成性、因果性は必ずしも考慮されていない。そこで、完全再構成性及び因果性を制約とした上で、フィルタ処理に対応する最小二乗問題を解くことで、プリリンギングが発生せず位相遅延のないモード分解を実現した。
【００４２】
すなわち、以下の制約付き最小二乗問題
【００４３】
【数１０】
2019159018.tif 000013
【００４４】
を解くことで、楽器音sをN個のモードxiへ分解する。実施例１と同様、nは時間インデックスを、τは楽器の発音時刻に対応するサンプルを表す。式（８）上段のコスト関数における周波数領域での重みwiをフィルタの周波数特性の逆hi-1の拡張として考えると、上段のコスト関数の最小化はフィルタ処理の一般化に対応する。特に、重みwiが非負の実数ベクトルの時、零位相のフィルタによるフィルタ処理の拡張に対応する。
【００４５】
重みwiを非負の実数ベクトルにし、式（８）下段で完全再構成性を満たす制約を導入したことで位相遅延のないモード分解を実現する。また、式（８）下段の制約は各モードが発音時刻τより前に成分を持たないことを示している。これにより、プリリンギングの発生を防ぎ、各モードにおいて因果性を満たすことを示している。
【００４６】
式（８）における重みwiは各モードxiが所望の周波数特性を持つように任意に作成することができる。楽器音のモード分解を考えているため、各モードxiが狭帯域かつ他のモードの成分の混入が少ない信号にする重みwiを作成する。
【００４７】
まず、ARモデルを用いて楽器音sを近似し、ARモデルのパラメータから楽器音のモードmiに対応する極piを選ぶ。この時、極piの偏角は各モードの中心周波数fiを正規化したものに対応する。極piごとに極piとその複素共役pi−のみを極に持つ共振フィルタを設計し、この共振フィルタの振幅応答を正規化したものを|hi~|とする。また、ARモデルで近似した楽器音sスペクトルで各モードの中心周波数fiにあたるピークの大きさをgiとする。
【００４８】
これらの情報を利用し、各モードxiを中心周波数fiの狭帯域な信号にする第１の重みwi(1)と、他のモードの混入を防ぐ第２の重みwi(2)を作成する。これら２つの重みのアダマール積wi(1)○wi(2)を後述する式（１４）の第３の重みwiとした。第１の重みwi(1)は
【００４９】
【数１１】
2019159018.tif 000014
【００５０】
とした。これにより、第１の重みwi(1)のうち、モードの中心周波数fiへの重みは0になり急峻なディップを持つ。また、第２の重みwi(2)を
【００５１】
【数１２】
2019159018.tif 000015
【００５２】
とする。これにより第２の重みwi(2)は他のモードの中心周波数fj≠iに対しピークを持つ。そのため、第３の重みwiはモードの中心周波数fiでディップを持ち、他のモード中心周波数fj≠iでピークを持つ。
【００５３】
式（８）の手法はフィルタ処理の一般化にあたるため、コスト関数でデータ忠実性を考慮している。これにより、各モードに他のモードの成分が混在するモードミキシングが発生してしまう。そこで本実施例では、フィルタ処理の一般化にあたるコスト関数ではなく、各モードの周波数領域で重み付けしたエネルギーをコスト関数にした問題を考える。そこで、楽器音のスペクトルs^、推定するモードのスペクトルをxi^を利用し、周波数領域での完全再構成性の制約のもと重み付けしたエネルギーの最小化問題
【００５４】
【数１３】
2019159018.tif 000016
【００５５】
を解くことで楽器音のスペクトルs^を各モードのスペクトルをxi^へ分解することを考える。但し、wi∈RLは各モードに対する重みを表す。離散フーリエ変換Fはユニタリ作用素であるため、周波数領域での完全再構成性が満たされる時、式（１１）下段の周波数領域の完全再構成性制約の両辺に逆離散フーリエ変換F-1をかけた時間領域での完全再構成性も満たされる。式（１１）上段のコスト関数は周波数インデックスξごとに独立であり、周波数インデックスξごとに考えることができる。式（１１）を周波数インデックスξごとに考えると
【００５６】
【数１４】
2019159018.tif 000017
【００５７】
と書き換えることができる。ここで、全ての周波数サンプルξに対し複数のモードの重みwi[ξ]が同時に0でない時、式（１２）の解は
【００５８】
【数１５】
2019159018.tif 000018
【００５９】
となる。式（１３）は、式（１１）の解xiが、周波数応答がGiとなるフィルタを楽器音sにかけたものと一致することを示している。つまり、式（１１）はモードの周波数領域で重み付けしたエネルギーの最小化を考えていたが、完全再構成性を考慮した位相遅延のないフィルタバンクG1,...,Nを作成し、そのフィルタバンクを楽器音sにかけることと一致している。
【００６０】
Gi[ξ]は重みw1,...,N[ξ]どうしの比によって決まり、重みwi[ξ]=0の時、Gi[ξ]=1となり、Gi≠j[ξ]=0となる。そのため、周波数応答がGiで与えられるフィルタは各モードの成分を排他的にする性質がある。そのため、式（１１）を解くことで、モードミキシングの少ないモード分解を実現することができる。
【００６１】
フィルタ処理を拡張した最適化問題によるモード分解より、周波数領域で重み付けした各モードのエネルギーの最小化によるモード分解（式（１１））の方がモードミキシングが少ないことが式（１１）の解析解（１３）より分かった。しかし、式（１１）では因果性を考慮していなかった。そこで、各モードの時間波形xiを変数とし、各モードの因果性を制約に加えた上で、周波数領域で重み付けした各モードのエネルギーの最小化問題
【００６２】
【数１６】
2019159018.tif 000019
【００６３】
を解くことによる楽器音のモード分解を考える。この問題は凸なので、適当な凸最適化アルゴリズムで解くことができる。
【００６４】
変数を各モードの時間波形xiで持ち、式（１４）下段の各モードの因果性を表す制約を導入することで、各モードxiにおけるプリリンギングの発生を防ぐ。これにより、各モードxiが完全再構成性、因果性を同時に満たした位相遅延のないモード分解が可能である。
【００６５】
本実施例では式（１４）をFast Iterative Shrinkage-Thresholding Algorithm (FISTA) を用いて解いた。FISTAを用いると式（１４）の解は
【００６６】
【数１７】
2019159018.tif 000020
【００６７】
を反復することで求められる。但し、kは反復のインデックス、μはステップサイズを表し、作用素Tは
【００６８】
【数１８】
2019159018.tif 000021
【００６９】
を表している。また、T*は作用素Tの随伴作用素を表す。κCは式（１４）の下段の完全再構成性と因果性表す制約を同時に満たすアフィン部分空間Cへの射影を表し、
【００７０】
【数１９】
2019159018.tif 000022
【００７１】
となる。式（１７）の上段は各モードxiの足し合わせを楽器音sに一致させる効果があり、完全再構成性の制約に対応している。一方、下段はプリリンギングの発生を防ぐ効果があり、因果性の制約に対応している。FISTAでは、式（１５）に示す通り、各反復ごとに各モードの時間波形xiに対し上記の制約を同時に満たすアフィン部分空間Cへの射影κCを行う。そのため、制約付き最小二乗問題（式（１４））をFISTAを用いて解くことで、各反復においてモードxiは必ず完全再構成性と因果性を満たす。
【００７２】
＜モード分解装置２＞
図５を参照して、本実施例のモード分解装置２の構成について説明する。同図に示すように、本実施例のモード分解装置２は、極取り出し部２１と、第１重み推定部２２と、第２重み推定部２３と、モード分解部２４を含む。
【００７３】
図６を参照して各部の動作について説明する。極取り出し部２１は、楽音信号sまたはその近似音から極を取り出す（Ｓ２１）。
【００７４】
第１重み推定部２２は、式（９）に基づいて、各モードxiを中心周波数fiの狭帯域な信号にする第１の重みwi(1)を推定する（Ｓ２２）。
【００７５】
第２重み推定部２３は、式（１０）に基づいて、他のモードの混入を防ぐ第２の重みwi(2)を推定する（Ｓ２３）。
【００７６】
モード分解部２４は、式（１４）の上段で表される最小二乗問題を、式（１４）の下段に表す完全再構成性、因果性を充たすことを制約として、FISTAを用いて求めることにより、楽音信号sを各モードx1,...,xNに分解する（Ｓ２４）。
【実施例３】
【００７７】
＜理論＞
以下、実施例３のモード分解装置３の基礎となる理論を説明する。一般に、楽器音sから周波数応答がhn^である線形フィルタを用いてモードxnを推定することは、
【００７８】
【数２０】
2019159018.tif 000023
【００７９】
と表される。ただし、s^=Fsは楽器音sのスペクトルを表し、F、○はそれぞれ離散フーリエ変換、アダマール積を表す。フィルタの周波数応答hn^に逆特性hn-1^が存在する時、式（１８）の解は最小二乗問題:
【００８０】
【数２１】
2019159018.tif 000024
【００８１】
の解と一致する。一般に知られている通り、線形なフィルタ処理では各モードの位相遅延とプリリンギングを同時に解消することができない。これらは分解されたモードの立ち上がりを変化させ、楽器音の立ち上がりの解析を困難にする。
【００８２】
そこで、式（１９）に対し、(I)完全再構成・因果性を表す制約の導入、(II)制約のみによるデータ忠実性の考慮、(III)アタック成分の変数への組み込み、という３点の変更を加えたモード分解手法を提案した。(I)により位相遅延及びプリリンギングを明示的に解消し、(III)においてモード分解の残差をアタック成分とする。上記手法において、凸最適化問題:
【００８３】
【数２２】
2019159018.tif 000025
【００８４】
を解くことで、楽器音sをN個のモードxnとアタック成分xN+1に分解する。ただし、実施例１と同様、F-1は逆離散フーリエ変換、tは時間インデックス、τは楽器の発音時刻に対応するサンプルを表す。しかし、上記手法では、モード数が多く複数のモードの中心周波数が近い実楽器音に適用する場合、計算時間と分解性能が課題となっていた。そこで本実施例では、式（２０）のADMMを用いた高速な解法（後述）を利用し、分解性能向上のための重みwiを改善することで（後述）、上記の課題を解決する。
【００８５】
ADMMを用いることで、式（２０）の高速な解法を実現する。式（２０）の解はADMMを用いて任意の初期値z[0]^、u[0]^から
【００８６】
【数２３】
2019159018.tif 000026
【００８７】
を反復することで求められる。ただし、
【００８８】
【数２４】
2019159018.tif 000027
【００８９】
とし、xT^はx^の転置を、kは反復インデックスを表す。また、
【００９０】
【数２５】
2019159018.tif 000028
【００９１】
とする。ただし、C1、C2はそれぞれ完全再構成性と因果性を表す制約を満たす空間、χCは空間Cに対する指示関数、PC2は空間C2への距離射影、proxρ,fは関数fの近接作用素を表し、それぞれ
【００９２】
【数２６】
2019159018.tif 000029
【００９３】
で与えられる。式（２８）は解析的に解くことができ
【００９４】
【数２７】
2019159018.tif 000030
【００９５】
と表せる。ただし、ynξ^=znξ^-xnξ^、ξは周波数インデックスを表す。すなわち、式（２９）〜（３３）を用いて式（２１）を計算し、式（２６）、（２７）を用いて式（２２）を計算することが可能である。これにより、式（２０）の解を効率的に求めることができる。
【００９６】
式（２０）における重みwnは各モード及びアタック成分xnが所望の周波数特性を持つよう任意に作成することができる。本実施例では、自己回帰モデルを用いて楽器音sを近似し、各モードに対応する極の情報から共振フィルタを設計する。共振フィルタの正規化した振幅応答|hn~|及び振幅応答の最大値anを利用し、各モードxnを狭帯域な信号にする第１の重みwndipと、他のモードの混入を防ぐ第２の重みwnpeaksを作成する。式（２０）の第３の重みwnを２つの重みのアダマール積wndip○wnpeaksで表し、第１の重みwndip、第２の重みwnpeaksをそれぞれ
【００９７】
【数２８】
2019159018.tif 000031
【００９８】
とする。ただし、μ>0は第１の重みwndipによって各モードxnをどれだけ狭帯域にするか調整するパラメータを表す。本実施例では第１の重みwndip、第２の重みwnpeaksはそれぞれ目的のモードの中心周波数に対する共振フィルタの振幅応答、他のモードの中心周波数に対するノッチフィルタの振幅応答の和に対応する。また、アタック成分に対する重みwN+1を
【００９９】
【数２９】
2019159018.tif 000032
【０１００】
とし、λ>0はアタック成分の大きさを調整するパラメータを表す。
【０１０１】
＜モード分解装置３＞
図７を参照して、本実施例のモード分解装置３の構成について説明する。同図に示すように、本実施例のモード分解装置３は、極取り出し部３１と、第１重み推定部３２と、第２重み推定部３３と、アタック成分重み推定部３４と、モード分解部３５を含む。
【０１０２】
図８を参照して各部の動作について説明する。極取り出し部３１は、楽音信号sまたはその近似音から極を取り出す（Ｓ３１）。
【０１０３】
第１重み推定部３２は、式（３４）に基づいて、第１の重みwndipを推定する（Ｓ３２）。
【０１０４】
第２重み推定部３３は、式（３５）に基づいて、第２の重みwnpeaksを推定する（Ｓ３３）。
【０１０５】
アタック成分重み推定部３４は、式（３６）に基づいて、アタック成分重みwN+1を推定する（Ｓ３４）。
【０１０６】
モード分解部３５は、式（２０）の上段で表される最小二乗問題を、式（２０）の下段に表す完全再構成性、因果性を充たすことを制約として、ADMMを用いて求めることにより、楽音信号sを各モードx1,...,xNとアタック成分xN+1に分解する（Ｓ３５）。
【０１０７】
＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ−ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【０１０８】
ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。
【０１０９】
ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。
【０１１０】
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１１１】
既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１１２】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１１３】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１１４】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１１５】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
2019159018.tif 000041
【図２】
2019159018.tif 000042
【図３】
2019159018.tif 000043
【図４】
2019159018.tif 000044
【図５】
2019159018.tif 000045
【図６】
2019159018.tif 000046
【図７】
2019159018.tif 000047
【図８】
2019159018.tif 000048
【手続補正書】
【提出日】平成31年2月5日(2019.2.5)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】００８４
【補正方法】変更
【補正の内容】
【００８４】
を解くことで、楽器音sをN個のモードxnとアタック成分xN+1に分解する。ただし、実施例１と同様、F-1は逆離散フーリエ変換、tは時間インデックス、τは楽器の発音時刻に対応するサンプルを表す。しかし、上記手法では、モード数が多く複数のモードの中心周波数が近い実楽器音に適用する場合、計算時間と分解性能が課題となっていた。そこで本実施例では、式（２０）のADMM（Alternating Direction Method of Multipliers）を用いた高速な解法（後述）を利用し、分解性能向上のための重みwiを改善することで（後述）、上記の課題を解決する。
【手続補正２】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、前記楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、前記楽音信号を各モードに分解するモード分解装置。
【請求項２】
請求項１に記載のモード分解装置であって、
楽音信号sまたはその近似音から極を取り出す極取り出し部と、
iを前記楽音信号sの各モードのインデックスとし、前記極の情報に基づいて定まるi番目の共振周波数fiの共振フィルタの振幅応答|hi|に基づいて、i番目の第１の重みwi(1)を推定する動作を各モードについて繰り返し実行する第１重み推定部と、
jをj≠iを満たす楽音信号sのモードのインデックスとし、j番目の共振周波数fjの共振フィルタの振幅応答|hj|に基づいて、i番目の第２の重みwi(2)を推定する動作を各モードについて繰り返し実行する第２重み推定部と、
残差rに対する重みである残差重みwrを推定する残差重み推定部と、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積であるi番目の第３の重みwiと、i番目の前記モードxiの周波数成分F(xi)と、前記楽音信号sの周波数成分F(s)と、パラメータλと、前記残差重みwrと、前記残差rの周波数成分F(r)に関する最小二乗問題
【数３０】
2019159018.tif 000049
の解を、完全再構成性と因果性を充たすことを制約として求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解部
を含むモード分解装置。
【請求項３】
請求項２に記載のモード分解装置であって、
前記第１重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、μをパラメータとし、i番目の第１の重みwi(1)[k]を、
【数３１】
2019159018.tif 000050
と推定する
モード分解装置。
【請求項４】
請求項２または３に記載のモード分解装置であって、
前記第２重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、i番目の第２の重みwi(2)[k]を
【数３２】
2019159018.tif 000051
と推定する
モード分解装置。
【請求項５】
請求項２から４の何れかに記載のモード分解装置であって、
前記残差重み推定部は、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、残差重みwr[k]を
【数３３】
2019159018.tif 000052
と推定する
モード分解装置。
【請求項６】
i,jを楽音信号sの各モードのインデックスとし、
ARモデルを用いて近似した楽音信号sから極piを取り出す極取り出し部と、
極piとその複素共役pi−のみを極に持つ共振フィルタの振幅応答を正規化したものを|hi~|とし、ξを周波数インデックスとして、
【数３４】
2019159018.tif 000053
に基づいて、各モードxiを中心周波数fiの狭帯域な信号にする第１の重みwi(1)を推定する第１重み推定部と、
前記ARモデルで近似した楽器音sスペクトルで各モードの中心周波数fiにあたるピークの大きさをgiとし、j番目の共振周波数fjの共振フィルタの振幅応答を|hj|とし、
【数３５】
2019159018.tif 000054
に基づいて、他のモードの混入を防ぐ第２の重みwi(2)を推定する第２重み推定部と、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積をi番目の第３の重みwiとし、i番目の前記モードxiの周波数成分F(xi)に関する最小二乗問題
【数３６】
2019159018.tif 000055
の解を、完全再構成性と因果性を充たすことを制約として、FISTAを用いて求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解部
を含むモード分解装置。
【請求項７】
n,mを楽音信号sの各モードのインデックスとし、
楽音信号sまたはその近似音から極を取り出す極取り出し部と、
共振フィルタの正規化した振幅応答を|hn~|とし、振幅応答の最大値をanとし、μは第１の重みwndipによって各モードxnをどれだけ狭帯域にするかを調整するパラメータとし、ξを周波数インデックスとして、
【数３７】
2019159018.tif 000056
に基づいて、各モードxnを狭帯域な信号にする第１の重みwndipを推定する第１重み推定部と、
【数３８】
2019159018.tif 000057
に基づいて、他のモードの混入を防ぐ第２の重みwnpeaksを推定する第２重み推定部と、
λをアタック成分の大きさを調整するパラメータとし、
【数３９】
2019159018.tif 000058
に基づいて、アタック成分に対する重みwN+1を推定するアタック成分重み推定部と、
前記第１の重みwndipと前記第２の重みwnpeaksのアダマール積を第３の重みwnとし、
【数４０】
2019159018.tif 000059
の最小二乗問題を、完全再構成性、因果性を充たすことを制約として、ADMMを用いて求めることにより、楽音信号sを各モードx1,...,xNとアタック成分xN+1に分解するモード分解部
を含むモード分解装置。
【請求項８】
モード分解装置が実行するモード分解方法であって、
楽音信号を各モードに分解するフィルタの振幅応答の逆特性と各モードの周波数成分のアダマール積と、前記楽音信号の周波数成分に関する最小二乗問題であって、完全再構成性と因果性を充たすことを制約とした制約付き最小二乗問題を解くことにより、前記楽音信号を各モードに分解するモード分解方法。
【請求項９】
請求項８に記載のモード分解方法であって、
楽音信号sまたはその近似音から極を取り出す極取り出しステップと、
iを前記楽音信号sの各モードのインデックスとし、前記極の情報に基づいて定まるi番目の共振周波数fiの共振フィルタの振幅応答|hi|に基づいて、i番目の第１の重みwi(1)を推定する動作を各モードについて繰り返し実行する第１重み推定ステップと、
jをj≠iを満たす楽音信号sのモードのインデックスとし、j番目の共振周波数fjの共振フィルタの振幅応答|hj|に基づいて、i番目の第２の重みwi(2)を推定する動作を各モードについて繰り返し実行する第２重み推定ステップと、
残差rに対する重みである残差重みwrを推定する残差重み推定ステップと、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積であるi番目の第３の重みwiと、i番目の前記モードxiの周波数成分F(xi)と、前記楽音信号sの周波数成分F(s)と、パラメータλと、前記残差重みwrと、前記残差rの周波数成分F(r)に関する最小二乗問題
【数４１】
2019159018.tif 000060
の解を、完全再構成性と因果性を充たすことを制約として求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解ステップ
を含むモード分解方法。
【請求項１０】
請求項９に記載のモード分解方法であって、
前記第１重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、μをパラメータとし、i番目の第１の重みwi(1)[k]を、
【数４２】
2019159018.tif 000061
と推定する
モード分解方法。
【請求項１１】
請求項９または１０に記載のモード分解方法であって、
前記第２重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、i番目の第２の重みwi(2)[k]を
【数４３】
2019159018.tif 000062
と推定する
モード分解方法。
【請求項１２】
請求項９から１１の何れかに記載のモード分解方法であって、
前記残差重み推定ステップは、
kをインデックスを表すものとし、ξ[k]をインデックスkに対応する周波数とし、gjをj番目のモードxjに対応するピークの大きさとし、残差重みwr[k]を
【数４４】
2019159018.tif 000063
と推定する
モード分解方法。
【請求項１３】
モード分解装置が実行するモード分解方法であって、
i,jを楽音信号sの各モードのインデックスとし、
ARモデルを用いて近似した楽音信号sから極piを取り出す極取り出しステップと、
極piとその複素共役pi−のみを極に持つ共振フィルタの振幅応答を正規化したものを|hi~|とし、ξを周波数インデックスとして、
【数４５】
2019159018.tif 000064
に基づいて、各モードxiを中心周波数fiの狭帯域な信号にする第１の重みwi(1)を推定する第１重み推定ステップと、
前記ARモデルで近似した楽器音sスペクトルで各モードの中心周波数fiにあたるピークの大きさをgiとし、j番目の共振周波数fjの共振フィルタの振幅応答を|hj|とし、
【数４６】
2019159018.tif 000065
に基づいて、他のモードの混入を防ぐ第２の重みwi(2)を推定する第２重み推定ステップと、
前記i番目の第１の重みwi(1)と前記i番目の第２の重みwi(2)のアダマール積をi番目の第３の重みwiとし、i番目の前記モードxiの周波数成分F(xi)に関する最小二乗問題
【数４７】
2019159018.tif 000066
の解を、完全再構成性と因果性を充たすことを制約として、FISTAを用いて求めることにより、前記楽音信号sを各モードx1,...,xNに分解するモード分解ステップ
を含むモード分解方法。
【請求項１４】
モード分解装置が実行するモード分解方法であって、
n,mを楽音信号sの各モードのインデックスとし、
楽音信号sまたはその近似音から極を取り出す極取り出しステップと、
共振フィルタの正規化した振幅応答を|hn~|とし、振幅応答の最大値をanとし、μは第１の重みwndipによって各モードxnをどれだけ狭帯域にするかを調整するパラメータとし、ξを周波数インデックスとして、
【数４８】
2019159018.tif 000067
に基づいて、各モードxnを狭帯域な信号にする第１の重みwndipを推定する第１重み推定ステップと、
【数４９】
2019159018.tif 000068
に基づいて、他のモードの混入を防ぐ第２の重みwnpeaksを推定する第２重み推定ステップと、
λをアタック成分の大きさを調整するパラメータとし、
【数５０】
2019159018.tif 000069
に基づいて、アタック成分に対する重みwN+1を推定するアタック成分重み推定ステップと、
前記第１の重みwndipと前記第２の重みwnpeaksのアダマール積を第３の重みwnとし、
【数５１】
2019159018.tif 000070
の最小二乗問題を、完全再構成性、因果性を充たすことを制約として、ADMMを用いて求めることにより、楽音信号sを各モードx1,...,xNとアタック成分xN+1に分解するモード分解ステップ
を含むモード分解方法。
【請求項１５】
コンピュータを請求項１から７の何れかに記載のモード分解装置として機能させるプログラム。