(57)【特許請求の範囲】
【請求項１】
複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力する混合信号時間周波数変換部と、
前記観測時間周波数成分Ｙに基づいて、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk,tを表すアクティベーションＵとの積Ｘと、前記観測時間周波数成分Ｙとの誤差、及び前記アクティベーションＵと、各基底ｋのパワーＵk,tに周期性を仮定して表現されるアクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨと、前記アクティベーションＵを推定する教師なし音源分離適用部と、
を含む音源分離装置。
【請求項２】
複数の音源ｉの各々に対して、前記音源ｉからの音源信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｓiω,tを表す観測時間周波数成分Ｓiを出力する分離信号時間周波数変換部と、
前記複数の音源ｉの各々に対して、前記観測時間周波数成分Ｓiに基づいて、各基底ｋiの、各周波数ωにおけるパワースペクトルＨiω,kを表す基底スペクトルＨiと、各基底ｋiの、各時刻ｔにおけるパワーＵik,tを表すアクティベーションＵiとの積Ｘiと、前記観測時間周波数成分Ｓiとの誤差、及び各基底ｋiのパワーＵik,tに周期性を仮定したときの、前記アクティベーションＵiとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨiと、前記アクティベーションＵiを推定する分離信号音源分離適用部と、
前記複数の音源ｉの各々に対して推定された前記基底スペクトルＨiを統合して、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを生成する基底スペクトル統合部と、
複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力する混合信号時間周波数変換部と、
前記観測時間周波数成分Ｙに基づいて、前記基底スペクトル統合部によって生成された前記基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、及び前記アクティベーションＵと、各基底ｋのパワーＵk,tに周期性を仮定して表現されるアクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記アクティベーションＵを推定する混合信号音源分離適用部と、
を含む音源分離装置。
【請求項３】
前記目的関数は、各基底ｋのパワーＵk,tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、前記周期関数のスパース化のためのスパース化項を更に含む請求項１又は２記載の音源分離装置。
【請求項４】
前記目的関数に含まれる前記スパース化項は、各基底ｋのパワーＵk,tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、同じ基本周波数で表わされる前記周期関数のグループのスパース化のためのものである請求項３記載の音源分離装置。
【請求項５】
混合信号時間周波数変換部及び教師なし音源分離適用部を含む音源分離装置における音源分離方法であって、
前記混合信号時間周波数変換部が、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力し、
前記教師なし音源分離適用部が、前記観測時間周波数成分Ｙに基づいて、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk,tを表すアクティベーションＵとの積Ｘと、前記観測時間周波数成分Ｙとの誤差、及び前記アクティベーションＵと、各基底ｋのパワーＵk,tに周期性を仮定して表現されるアクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨと、前記アクティベーションＵを推定する
音源分離方法。
【請求項６】
分離信号時間周波数変換部、分離信号音源分離適用部、基底スペクトル統合部、混合信号時間周波数変換部、及び混合信号音源分離適用部を含む音源分離装置における音源分離方法であって、
前記分離信号時間周波数変換部が、複数の音源ｉの各々に対して、前記音源ｉからの音源信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｓiω,tを表す観測時間周波数成分Ｓiを出力し、
前記分離信号音源分離適用部が、前記複数の音源ｉの各々に対して、前記観測時間周波数成分Ｓiに基づいて、各基底ｋiの、各周波数ωにおけるパワースペクトルＨiω,kを表す基底スペクトルＨiと、各基底ｋiの、各時刻ｔにおけるパワーＵik,tを表すアクティベーションＵiとの積Ｘiと、前記観測時間周波数成分Ｓiとの誤差、及び各基底ｋiのパワーＵik,tに周期性を仮定したときの、前記アクティベーションＵiとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨiと、前記アクティベーションＵiを推定し、
前記基底スペクトル統合部が、前記複数の音源ｉの各々に対して推定された前記基底スペクトルＨiを統合して、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを生成し、
前記混合信号時間周波数変換部が、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力し、
前記混合信号音源分離適用部が、前記観測時間周波数成分Ｙに基づいて、前記基底スペクトル統合部によって生成された前記基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk,tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、及び前記アクティベーションＵと、各基底ｋのパワーＵk,tに周期性を仮定して表現されるアクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記アクティベーションＵを推定する
音源分離方法。
【請求項７】
前記目的関数は、各基底ｋのパワーＵk,tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、前記周期関数のスパース化のためのスパース化項を更に含む請求項５又は６記載の音源分離方法。
【請求項８】
前記目的関数に含まれる前記スパース化項は、各基底ｋのパワーＵk,tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、同じ基本周波数で表わされる前記周期関数のグループのスパース化のためのものである請求項７記載の音源分離方法。
【請求項９】
コンピュータを、請求項１〜請求項４の何れか１項記載の音源分離装置を構成する各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音源分離装置、方法、及びプログラムに係り、特に、複数の音源からの音源信号が混合された混合信号を分離する音源分離装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
音楽音響信号の多くは、複数の例えば楽器などの音源から構成されており、混合信号の分離は音楽構造理解や自動採譜、楽器やピッチの部分的な編集、雑音除去、欠損値補間など多くの場面で必要になる。近年、音響信号処理分野において、自動採譜や音源分離などへの有効性から非負値行列分解（Non-negative Matrix Factorization, NMF）への注目が高まっている。NMFでは、混合信号（観測信号）のスペクトログラムを非負値行列Ｙで表現した上で、2つの非負値行列Ｈ、Ｕの積の形に分解する（図１参照）。これは、各時間窓で観測されたスペクトルを、時間変化する係数行列U（アクティベーション）の値でスケールされた基底スペクトルＨの線形和で近似していることに該当し、観測行列はランク1 のスペクトログラムＨω,kＵk,tの和の形で表現される。
【０００３】
NMFの重要な特徴として、非負値制約により分解行列がスパースな値を持つことが挙げられる。例えば、Ｕには多くの０要素が含まれることが多い。この特徴により、各時間窓の観測スペクトルは、多くの場合少数の基底スペクトルの和で表現される。従って、各基底スペクトルには観測信号全体において複数回出現する典型的なスペクトル（多くの場合は頻出する楽音）が抽出され、効率の良い分解表現が可能になる。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】S.A. Raczynski, N.Ono, and S.Sagayama. Multipitch Analysis with Harmonic Nonnegative Matrix Approximation. ISMIR 2007, 8th International Conference on Music Information Retrieval. Citeseer, 2007.
【非特許文献２】E.Vincent, N.Bertin, and R.Badeau. Harmonic and Inharmonic Nonnegative Matrix Factorization for Polyphonic Pitch Transcription. IEEE International Conference on Acoustics, Speech and Signal Processing, 2008. ICASSP 2008., pages 109-112. IEEE, 2008.
【非特許文献３】H.Kameoka, M.Nakano, K.Ochiai, Y.Imoto, K.Kashino, and S.Sagayama. Constrained and Regularized Variants of Non-negative Matrix Factorization Incorporating Music-specic Constraints. ICASSP, pages 5365-5368, 2012.
【非特許文献４】T.Virtanen. Monaural Sound Source Separation by Nonnegative Matrix Factorization with Temporal Continuity and Sparseness Criteria. IEEE Transactions on Audio, Speech, and Language Processing, 15(3):1066-1074, 2007.
【非特許文献５】M.Nakano, J.LeRoux, H.Kameoka, Y.Kitano, N.Ono, and S.Sagayama. Nonnegative Matrix Factorization with Markov-chained Bases for Modeling Time-varying Patterns in Music Spectrograms. Latent Variable Analysis and Signal Separation, pages 149-156. Springer, 2010.
【非特許文献６】K.Ochiai, H.Kameoka, and S.Sagayama. Explicit Beat Structure Modeling for Non-negative Matrix Factorization-based Multipitch Analysis. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 133-136. IEEE, 2012.
【非特許文献７】Rafii Zafar and Bryan Pardo. Repeating pattern extraction technique (REPET): A simple method for music/voice separation, IEEE Transactions on Audio, Speech, and Language Processing, Vol. 21, no.1, pp. 73-84, 2013.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、音楽スペクトログラムにNMFを適用する際には、以下に示す２つの課題がある。第１に、精度と、分離信号の有意性との間にトレードオフが発生する。NMFではＸ＝ＨＵと観測行列Ｙの間の近似誤差を目的関数として定義し、その目的関数を最小化する。基底数を増やすことにより近似誤差は減少するが、基底数を増やしすぎると分解行列の解釈性が下がる。例えば、極端な例としてＹ＝ＨＩやＹ＝ＩＵなどの無意味な分解結果が出力されてしまう場合がある。ここでＩは単位行列を示す。したがって、多くの基底を用いて精度（少ない近似誤差）を保ちながら、有意性の高い分離を行うためには、非負値性以外の有用な制約をＨＵに課したモデルが必要になる。
【０００６】
NMFの第２の課題として、局所的最適解問題が上げられる。既存のNMFの最適化アルゴリズムの多くは収束が保証されているものの、ある限られた領域においてのみ最適な局所的最適解へ収束してしまう場合が多く、必ずしも全ての領域において最適な大域的最適解が求まるとは限らない。音楽音響信号への適用において、この局所的最適解は、例えば分離信号が聴覚で捉えられる音楽事象に対応しないなど、音楽的に解釈不可能な解となり得る。
【０００７】
上記の問題を解決するには、データの特徴を反映した制約をモデルＨＵや目的関数に組み込むことが考えられる。例えば、音楽の和声構造を各基底スペクトルに仮定した手法（非特許文献１、２）や、音色の類似性を考慮して基底をクラスタリングする枠組みを導入した手法（非特許文献３）が提案されている。音楽の時間変化特徴を考慮した制約付NMF も提案されており、発音が一定時間持続することを仮定したモデル（非特許文献４）や、「アタック」「ディケイ」「サステイン」「リリース」などの楽音の状態変化を基底スペクトルに仮定したモデル（非特許文献５）、ビート（拍）に応じたリズムをアクティベーションに仮定した制約付NMF（非特許文献６）が提案されている。単純かつ単一な周期成分を抽出する手法（非特許文献７）も提案されているが、周期性は単一であるとは限らず，複数の音源で異なる複数の周期性が存在する可能性がある。
【０００８】
本発明では、上記問題を解決するために成されたものであり、周期的な発音パターンを抽出することができる音源分離装置、方法、及びプログラムを提供することを目的とする。周期性を想定した一部の基底で、どの音源がいつ発音するかを示すアクティベーションに基底毎に異なる周期性を仮定し、目的関数に制約として加えることにより，周期的な発音パターンの抽出を促進する。
【課題を解決するための手段】
【０００９】
上記目的を達成するために、第１の発明に係る音源分離装置は、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力する混合信号時間周波数変換部と、前記観測時間周波数成分Ｙに基づいて、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、前記観測時間周波数成分Ｙとの誤差、及び各基底ｋのパワーＵk.tに周期性を仮定したときの、前記アクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨと、前記アクティベーションＵを推定する教師なし音源分離適用部と、を含んで構成されている。
【００１０】
第２の発明に係る音源分離方法は、混合信号時間周波数変換部及び教師なし音源分離適用部を含む音源分離装置における音源分離方法であって、前記混合信号時間周波数変換部が、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力し、前記教師なし音源分離適用部が、前記観測時間周波数成分Ｙに基づいて、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、前記観測時間周波数成分Ｙとの誤差、及び各基底ｋのパワーＵk.tに周期性を仮定したときの、前記アクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨと、前記アクティベーションＵを推定する。
【００１１】
第３の発明に係る音源分離装置は、複数の音源ｉの各々に対して、前記音源ｉからの音源信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｓiω,tを表す観測時間周波数成分Ｓiを出力する分離信号時間周波数変換部と、前記複数の音源ｉの各々に対して、前記観測時間周波数成分Ｓiに基づいて、各基底ｋiの、各周波数ωにおけるパワースペクトルＨiω,kを表す基底スペクトルＨiと、各基底ｋiの、各時刻ｔにおけるパワーＵik.tを表すアクティベーションＵiとの積Ｘiと、前記観測時間周波数成分Ｓiとの誤差、及び各基底ｋiのパワーＵik.tに周期性を仮定したときの、前記アクティベーションＵiとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨiと、前記アクティベーションＵiを推定する分離信号音源分離適用部と、前記複数の音源ｉの各々に対して推定された前記基底スペクトルＨiを統合して、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを生成する基底スペクトル統合部と、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力する混合信号時間周波数変換部と、前記観測時間周波数成分Ｙに基づいて、前記基底スペクトル統合部によって生成された前記基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、及び各基底ｋのパワーＵk.tに周期性を仮定したときの、前記アクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記アクティベーションＵを推定する混合信号音源分離適用部と、を含んで構成されている。
【００１２】
第４の発明に係る音源分離方法は、分離信号時間周波数変換部、分離信号音源分離適用部、基底スペクトル統合部、混合信号時間周波数変換部、及び混合信号音源分離適用部を含む音源分離装置における音源分離方法であって、前記分離信号時間周波数変換部が、複数の音源ｉの各々に対して、前記音源ｉからの音源信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｓiω,tを表す観測時間周波数成分Ｓiを出力し、前記分離信号音源分離適用部が、前記複数の音源ｉの各々に対して、前記観測時間周波数成分Ｓiに基づいて、各基底ｋiの、各周波数ωにおけるパワースペクトルＨiω,kを表す基底スペクトルＨiと、各基底ｋiの、各時刻ｔにおけるパワーＵik.tを表すアクティベーションＵiとの積Ｘiと、前記観測時間周波数成分Ｓiとの誤差、及び各基底ｋiのパワーＵik.tに周期性を仮定したときの、前記アクティベーションＵiとの誤差を用いて表わされた目的関数を最適化するように、前記基底スペクトルＨiと、前記アクティベーションＵiを推定し、前記基底スペクトル統合部が、前記複数の音源ｉの各々に対して推定された前記基底スペクトルＨiを統合して、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを生成し、前記混合信号時間周波数変換部が、複数の音源からの音源信号が混合された混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力し、前記混合信号音源分離適用部が、前記観測時間周波数成分Ｙに基づいて、前記基底スペクトル統合部によって生成された前記基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、及び各基底ｋのパワーＵk.tに周期性を仮定したときの、前記アクティベーションＵとの誤差を用いて表わされた目的関数を最適化するように、前記アクティベーションＵを推定する。
【００１３】
また、本発明のプログラムは、コンピュータを、上記の音源分離装置を構成する各部として機能させるためのプログラムである。
【発明の効果】
【００１４】
以上説明したように、本発明の音源分離装置、方法、及びプログラムによれば、基底スペクトルＨと、アクティベーションＵとの積Ｘと、観測時間周波数成分との誤差、各基底ｋのパワーＵk.tに周期性を仮定したときの、アクティベーションＵとの誤差、及び周期関数のスパース化のためのスパース化項を用いて表わされた目的関数を最適化するように、基底スペクトルＨと、アクティベーションＵを推定することにより、周期的な発音パターンを抽出することができる。
【図面の簡単な説明】
【００１５】
【図１】ＮＭＦの表現を示す図である。
【図２】アクティベーションＵk,tとα、βの推定結果の例を示す図である。
【図３】本発明の第１の実施の形態に係る音源分離装置の機能的構成を示すブロック図である。
【図４】本発明の第１の実施の形態に係る音源分離装置の教師あり音源分離適用部の構成を示すブロック図である。
【図５】本発明の第１、第２の実施の形態に係る音源分離装置の分離信号音源分離適用部又は混合信号音源分離適用部の構成を示すブロック図である。
【図６】本発明の第１の実施の形態に係る音源分離装置における音源分離処理ルーチンを示すフローチャート図である。
【図７】本発明の第２の実施の形態に係る音源分離装置の機能的構成を示すブロック図である。
【図８】本発明の第２の実施の形態に係る音源分離装置の教師なし音源分離適用部の構成を示すブロック図である。
【図９】本発明の第２の実施の形態に係る音源分離装置における音源分離処理ルーチンを示すフローチャート図である。
【図１０】提案手法による音源分離結果を示す図である。
【図１１】従来のＮＭＦによる音源分離結果を示す図である。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１７】
＜本発明の実施の形態の概要＞
まず、本実施の形態における概要について説明する。本実施の形態では、各音符やドラム音の発音タイミングにおける、リズムの周期性を仮定して、制約付NMFを提案する。特にベースやドラムなどの楽器で周期性が想定される。どの音源がいつ発音するかを示すアクティベーションに周期性を仮定し、目的関数に制約として加えることにより、周期的な発音パターンの抽出を促進する。周期性を仮定する制約を加えた目的関数において、最適解を求めるアルゴリズムを定式化した上で、教師あり・教師なし両方の枠組みにおいて音源分離に適用する方法を提案する。
【００１８】
＜提案モデルの定式化＞
本実施の形態では、周期を考慮したパターン抽出方法（提案手法）を非負値周期成分分析手法（Non-negative Periodic Component Analysis, NPCA）と呼ぶ。NPCAの目的関数は以下のように定義される。
【００１９】
【数１】
0006535611.tif 000002
【００２０】
w1,w2＞0 は正則化パラメータである。
【００２１】
Ｉ（Ｈ，Ｕ）は通常のNMFの目的関数で考慮されるものと同じ、観測行列Ｙと、NMFモデルによる分解行列の積Ｘ＝ＨＵとの近似誤差である。
【００２２】
以下のようなNMFモデルを考える。
【００２３】
【数２】
0006535611.tif 000003
ここで、kは基底インデックス、Kは基底数を示し、ωとtはそれぞれ周波数と時間窓を示す。ＹとＸとの間の距離尺度としてＩダイバージェンス（参考文献１を参照）を用いると、Ｉ(Ｈ,Ｕ)は以下のように書ける。
【００２４】
【数３】
0006535611.tif 000004
［参考文献１］：I.Csiszar. I-divergence Geometry of Probability Distributions and Minimization Problems. The Annals of Probability, pages 146-158, 1975.
【００２５】
Ｊ（Ｕ, α, β）はアクティベーションＵk,1, ... , Ｕk,T に周期性を仮定するための制約を表す項である。T は時間-周波数表現における時間窓数を示す。提案法では、Ｕの非負性を保つために、Ｕk,tの対数を取ったものが、倍音構造を考慮した正弦曲線の重み付き和で近似されること、すなわち
【００２６】
【数４】
0006535611.tif 000005
と表現できることを仮定する。提案手法では、
【００２７】
【数５】
0006535611.tif 000006
【００２８】
が、周期1/Pmの繰り返しパターンの形状を決定する。
【００２９】
ここで、各正弦曲線の周期は基本周波数Pmの整数倍に限定される形を取っており、nは各周波数のn倍音を表現するための変数である。この制約により、
【００３０】
【数６】
0006535611.tif 000007
【００３１】
がちょうど周波数Pmの周期構造を持つことが保証される。以上により、log Ｕk,t はM 個の周期関数の和で表現される。なお、式(4) はＵを正弦曲線の重み付き和の指数の形で近似していると捉えられる。
【００３２】
周期{Pm}について、例えば
【００３３】
【数７】
0006535611.tif 000008
【００３４】
などと定義することが考えられる。また、倍音を表す定数の集合{n} について、{1, 2, ..., 9, 10} などの定義が考えられる。上述の例では、想定できる周期の数Mは20 となる。
【００３５】
提案モデルでは、自動的にＵk,tを最適に表現する周期を推定することができる。図2に、ある短い音楽データにおけるＵk,tと、Ｕk,tの太枠で囲まれた時間窓内に対応するαkm,n、βkm,nの推定例を示す。例においては、図２（Ａ）〜（Ｄ）で示されるいくつかのmとnの組み合わせに該当するαkm,nとβkm,nの値が他のmやnに比べて、大きい値を持つと推定された。下の周期関数の図は、抽出されたm, nに対応するｃｏｓ nＰmt またはｓｉｎ nＰmt のＵの太枠内の時間窓における離散曲線を示す。Sum と書かれた曲線は、図２（Ａ）〜（Ｄ）に該当する曲線の重み付き和を示し、Exp と書かれた尖った曲線は、Sum の指数を取ったものである。Exp の曲線は実際に推定されたＵk,tに近い形をしていることが分かる。
なお、図２（Ａ）、（Ｂ）などで発生しているエーリアシングを回避して、より明示的な周期を推定するためには、Ｐmの定義を式（５）から変更することも考えられる。最大周波数ＰMがナイキスト周波数以下になるように以下の式(５Ａ)でＰmを定義する。ここで
a=2π/T、b=(π-0.001)/Nである。加えて、式（５）や式（５Ａ）では周波数が少ないほど密に値を取るように設定していたが、これを以下の式（５Ｂ）のようにエーリアシングを回避した上で均等に分割することも考えられる。なお、式（５Ａ）、式（５Ｂ）でＰmを定義する際には、式（５）に比べてＰmの範囲が狭くなるため、Ｍの値を上述の20よりも少なくする方がよいと考えられる。
【数８】
0006535611.tif 000009
（５Ａ）
0006535611.tif 000010
（５Ｂ）
【００３６】
提案手法では、アクティベーションＵに周期性を仮定するための目的関数Ｊ(Ｕ, α, β) を、式(4) の両辺間のユークリッド距離を用いて、以下のように定義する。
【００３７】
【数９】
0006535611.tif 000011
【００３８】
ここで、
【００３９】
【数１０】
0006535611.tif 000012
【００４０】
と置いた。
【００４１】
式(4) でＵに周期性を仮定すると、m, n の組み合わせによっては、想定する正弦曲線が冗長になる場合がある。Ｕk,tを最小限の周期関数の重み付き和で表現するために、本発明の実施の形態ではα とβ に、以下に示すgroup sparsity を考慮する。
【００４２】
ここで、係数αとβを、同じ基本周波数を用いるＭ個の異なるグループに分け、group sparsity を誘導するための項
【００４３】
【数１１】
0006535611.tif 000013
【００４４】
を定義する。Ｌ(α, β) は、高々数個のα とβ のグループのみが非ゼロの値を取り、他の全てのグループでは値が0になっている場合、小さい値を取る。従って、Ｌ(α, β)を目的関数に含めて最小化することにより、Ｕk,tを最小限の周期で表現できる。図２はgroup sparsity を考慮したα とβ の推定結果を示している。
【００４５】
数種類のmに関連するα, β のグループでのみ、αとβの値が大きくなっている。グループｍごとに、倍音n方向の二乗和を最小化することにより、n方向の二乗和が0に近くなった場合、全てのnで該当するα, βの値が強制的に0に変換される。スパース性を保ちつつ、n方向では複数の係数が非ゼロになることを許容している。
【００４６】
なお、α, βをスパース化する際には、上述したgroup sparsity だけでなく、倍音方向、すなわちn方向にもスパース性を仮定する方法も考えられる。倍音成分も最小限に抑えることが可能になる。m, n両方向にスパース性を仮定する場合、sparsity 誘導項は
【００４７】
【数１２】
0006535611.tif 000014
【００４８】
となる。
【００４９】
本発明の実施の形態では、以上のように各基底スペクトルのアクティベーションに対して、十分な数の周期候補の中から適切に選択された、最小限の周期関数の和で表現される任意の周期を持つ周期的イベントを検出・推定した上で、その周期性を行列分解の枠組みで考慮することが可能になる。
【００５０】
＜最適化アルゴリズム＞
目的関数である式(1)を最小化するＨ,Ｕ,α,βの値を求めるアルゴリズムについて説明する。繰り返し計算により目的関数を最小化することができる。各ループでは、それぞれＨ、Ｕ、αとβについて最適化を行う3段階のステップを実施する。最適化を解析的に行うことは不可能であるが、補助関数法を用いて、目的関数の各項に対して設定した上界関数Ｉ+、Ｊ+、Ｌ+を最小化することにより、間接的にＨ,Ｕ,α,βの値を最適化する更新式を導出することができる。上界関数をＨ,Ｕ,α,βに関して最小化するステップと、目的関数と上界関数の値が等しくなるように補助変数を更新するステップを繰り返すことで、目的関数を単調に降下させることができる。
【００５１】
まず、Ｉ(Ｈ，Ｕ)の上界関数Ｉ+ については、通常のNMFと同様に、Jensen の不等式を用いて以下のように設定することができる。
【００５２】
【数１３】
0006535611.tif 000015
補助変数λは
【００５３】
【数１４】
0006535611.tif 000016
【００５４】
を満たす正の重み定数であり、等号成立条件は
【００５５】
【数１５】
0006535611.tif 000017
【００５６】
である。
【００５７】
次に、Ｊ(Ｕ, α, β) の上界関数Ｊ+ は以下のように設計できる。
【００５８】
【数１６】
0006535611.tif 000018
【００５９】
ここで、
【００６０】
【数１７】
0006535611.tif 000019
【００６１】
である。等号成立条件はξ＝Ｕk,t である。なお、上記の上界関数を設計するにあたって非特許文献８で証明されている以下の不等式を使用した。
【００６２】
【数１８】
0006535611.tif 000020
【００６３】
[非特許文献８]：樋口卓哉， 亀岡弘和. ケプストラム距離正則化に基づく多重音解析. 情報処理学会研究報告[音楽情報科学], 2014(10):1-6, 2014.
【００６４】
最後に、f(x)＝ｘp/2 に接する放物線がx＞0 のとき常にf(x) 以上の値を取ることを利用し、Ｌ(α, β) に関する上界関数Ｌ+ を以下のように設計する。
【００６５】
【数１９】
0006535611.tif 000021
0006535611.tif 000022
【００６６】
であり、等号成立条件は
【００６７】
【数２０】
0006535611.tif 000023
である。
【００６８】
続いて、上述した上界関数を用いてＨ，Ｕ， α, β の更新式を導出する。更新式は、上界関数を偏微分したものを0と置いて、各変数について解くことにより得られる。
【００６９】
Ｈ、Ｕの更新式は以下のように導出される。
【００７０】
【数２１】
0006535611.tif 000024
【００７１】
ここで、
【００７２】
【数２２】
0006535611.tif 000025
【００７３】
である。
【００７４】
α、β は、以下のようなαkm,n とβkm,n に関する連立方程式を解くことにより更新できる。
【００７５】
【数２３】
0006535611.tif 000026
（２２）
0006535611.tif 000027
（２３）
【００７６】
連立方程式
【００７７】
【数２４】
0006535611.tif 000028
【００７８】
の解は、例えばMatlab では
【００７９】
【数２５】
0006535611.tif 000029
【００８０】
、Python 言語のNumPy では
【００８１】
【数２６】
0006535611.tif 000030
【００８２】
、R 言語では
【００８３】
【数２７】
0006535611.tif 000031
【００８４】
とすることにより解くことができる。
【００８５】
＜本発明の第１の実施の形態に係る音源分離装置の構成＞
次に、本発明の第１の実施の形態に係る音源分離装置の構成について説明する。第１の実施の形態では、教師データとして分離信号が入手できる状況を想定して、教師あり音源分離を行う。図３に示すように、本発明の第１の実施の形態に係る音源分離装置１００は、ＣＰＵと、ＲＡＭと、後述する音源分離処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することができる。この音源分離装置１００は、機能的には図３に示すように入力部１０と、演算部２０と、出力部９０と、を含んで構成されている。
【００８６】
入力部１０は、複数の音源ｉの各々に対して、当該音源ｉからの音源信号の時系列データを受け付け、分離信号記憶部２２に記憶する。また、入力部１０は、複数の音源からの音源信号が混合された混合信号の時系列データを受け付け、混合信号記憶部２４に記憶する。
【００８７】
演算部２０は、分離信号記憶部２２と、混合信号記憶部２４と、分離信号入力受付部２６と、分離信号時間周波数変換部２８と、混合信号入力受付部３０と、混合信号時間周波数変換部３２と、教師あり音源分離適用部３４と、Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部３６と、ＩＦＦＴ適用部３８と、分離信号描画部４０と、分離信号再生部４２とを含んで構成されている。
【００８８】
分離信号記憶部２２には、入力部１０において受け付けた複数の音源の各々の分離信号の時系列データが記憶されている。なお、音源の一例は、楽器である。
【００８９】
混合信号記憶部２４には、入力部１０において受け付けた混合信号の時系列データが記憶されている。
【００９０】
分離信号入力受付部２６は、分離信号記憶部２２から、複数の音源の各々の分離信号の時系列データを受け付ける。
【００９１】
分離信号時間周波数変換部２８は、複数の音源ｉの各々に対して、分離信号入力受付部２６により受け付けた当該音源ｉからの音源信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｓiω,tを表す観測時間周波数成分Ｓiを出力する。なお、本実施の形態においては、音源ｉからの音源信号をフーリエ変換した複素スペクトログラムの絶対値を格納した観測時間周波数成分Ｓiを求める。
【００９２】
混合信号入力受付部３０は、混合信号記憶部２４から、混合信号の時系列データを受け付ける。
【００９３】
混合信号時間周波数変換部３２は、混合信号入力受付部３０により受け付けた混合信号の時系列データを入力として、各時刻ｔにおける各周波数ωの観測時間周波数成分Ｙω,tを表す観測時間周波数成分Ｙを出力する。
【００９４】
教師あり音源分離適用部３４は、図４に示すように、入力受付部５０、分離信号音源分離適用部５２、基底スペクトル統合部５４、混合信号音源分離適用部５６、及び結果出力部５８を含んで構成されている。
【００９５】
入力受付部５０は、分離信号時間周波数変換部２８によって、複数の音源ｉの各々に対して出力された、観測時間周波数成分Ｓiを受け付け、混合信号時間周波数変換部３２によって出力された観測時間周波数成分Ｙを受け付ける。
【００９６】
分離信号音源分離適用部５２は、複数の音源ｉの各々に対して、観測時間周波数成分Ｓiに基づいて、各基底ｋiの、各周波数ωにおけるパワースペクトルＨiω,kを表す基底スペクトルＨiと、各基底ｋiの、各時刻ｔにおけるパワーＵik.tを表すアクティベーションＵiとの積Ｘiと、前記観測時間周波数成分Ｓiとの誤差、各基底ｋiのパワーＵik.tに周期性を仮定したときの、前記アクティベーションＵiとの誤差、及び各基底ｋiのパワーＵik.tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、同じ基本周波数で表わされる前記周期関数のグループのスパース化のためのスパース化項を用いて表わされた、上記（１）式と同様の式で示される目的関数を最適化するように、基底スペクトルＨiと、アクティベーションＵi、周期性を表す変数α、βを推定する。なお、本実施の形態では、目的関数に含まれるＬ(α, β)として、上記（８）式で表わされるgroup sparsity を誘導するための項を用いる。
【００９７】
分離信号音源分離適用部５２は、図５に示すように、初期値設定部６０、変数更新部６２、補助変数更新部６４、及び収束判定部６６を含んで構成されている。
【００９８】
初期値設定部６０は、音源ｉに対して、基底スペクトルＨiと、アクティベーションＵiと、周期性を表す変数α、βと、補助変数λとに初期値を設定する。Ｈi，Ｕiの初期値には、例えば、通常のＮＭＦによる分解結果を使用すればよい。α、βの初期値については、ＮＰＣＡの目的関数のうち、α、βの値がないと計算できないＬ（α、β）に関する項を取り除いたものを用いて算出すればよい。
【００９９】
変数更新部６２は、音源ｉに対して、観測時間周波数成分Ｓiと、初期値設定部６０により初期値が設定された、又は変数更新部６２及び補助変数更新部６４により前回更新された、基底スペクトルＨiと、アクティベーションＵiと、周期性を表す変数α、βと、補助変数λとに基づいて、上記（１８）式、（１９）式、（２２）式、（２３）式と同様の式に従って、基底スペクトルＨiと、アクティベーションＵiと、周期性を表す変数α、βとを更新する。
【０１００】
補助変数更新部６４は、音源ｉに対して、変数更新部６２により更新された、基底スペクトルＨiと、アクティベーションＵiとに基づいて、上記（１１）式と同様の式に従って、補助変数λを更新する。
【０１０１】
収束判定部６６は、音源ｉに対して、予め定められた収束条件を満たすまで、変数更新部６２における更新処理と、補助変数更新部６４における更新処理とを繰り返させる。
【０１０２】
例えば、予め定めた繰り返し回数に到達した場合に、収束条件を満たすと判定する。
【０１０３】
音源信号音源分離適用部５２は、複数の音源ｉの各々に対して、上述した、初期値設定部６０、変数更新部６２、補助変数更新部６４、及び収束判定部６６の一連の処理を繰り返す。
【０１０４】
基底スペクトル統合部５４は、音源信号音源分離適用部５２によって複数の音源ｉの各々に対して推定された基底スペクトルＨiを統合して、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを生成する。
【０１０５】
例えば、基底スペクトル統合部５４は、音源ごとに推定された基底スペクトルＨiω,kをひとつの行列に統合し、混合信号の分離用の基底スペクトルＨω,kを作成する。この時、統合された基底スペクトルの基底数Ｋは
【０１０６】
【数２８】
0006535611.tif 000032
【０１０７】
個となる。
【０１０８】
混合信号音源分離適用部５６は、観測時間周波数成分Ｙに基づいて、基底スペクトル統合部５４によって生成された基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、各基底ｋのパワーＵk.tに周期性を仮定したときの、アクティベーションＵとの誤差、及び各基底ｋのパワーＵk.tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、同じ基本周波数で表わされる前記周期関数のグループのスパース化のためのスパース化項を用いて表わされた、上記（１）式に示す目的関数を用いて表わされた目的関数を最適化するように、アクティベーションＵ、周期性を表す変数α、βを推定する。
【０１０９】
この時、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨを更新せずに固定したまま、アクティベーションＵ、周期性を表す変数α、βのみを最適化する。
【０１１０】
混合信号音源分離適用部５６は、音源信号音源分離適用部５２と同様に、上記図５に示すように、初期値設定部６０、変数更新部６２、補助変数更新部６４、及び収束判定部６６を含んで構成されている。
【０１１１】
混合信号音源分離適用部５６の初期値設定部６０は、基底スペクトルＨとして、基底スペクトル統合部５４によって生成された基底スペクトルＨを設定すると共に、アクティベーションＵと、周期性を表す変数α、βと、補助変数λとに初期値を設定する。
【０１１２】
混合信号音源分離適用部５６の変数更新部６２は、観測時間周波数成分Ｙと、初期値設定部６０により初期値が設定された、又は変数更新部６２及び補助変数更新部６４により前回更新された、基底スペクトルＨと、アクティベーションＵと、周期性を表す変数α、βと、補助変数λとに基づいて、上記（１９）式、（２２）式、（２３）式に従って、アクティベーションＵと、周期性を表す変数α、βとを更新する。
【０１１３】
混合信号音源分離適用部５６の補助変数更新部６４は、基底スペクトルＨと、変数更新部６２により更新されたアクティベーションＵとに基づいて、上記（１１）式に従って、補助変数λを更新する。
【０１１４】
混合信号音源分離適用部５６の収束判定部６６は、予め定められた収束条件を満たすまで、変数更新部６２における更新処理と、補助変数更新部６４における更新処理とを繰り返させる。
【０１１５】
結果出力部５８は、基底スペクトル統合部５４によって生成された基底スペクトルＨと、混合信号音源分離適用部５６によって得られたアクティベーションＵとを、Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部３６へ出力する。
【０１１６】
Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部３６は、結果出力部５８から出力された基底スペクトルＨと、アクティベーションＵとに基づいて、Ｗｉｎｎｅｒ  ｆｉｌｔｅｒを作成した上で、混合信号から得られた複素スペクトログラムである観測時間周波数成分ＹにＷｉｎｎｅｒ  ｆｉｌｔｅｒをかけて分離信号の複素スペクトログラムを計算する。
【０１１７】
ここで、分離信号の複素スペクトログラム^siω,tが以下のように計算される。
【０１１８】
【数２９】
0006535611.tif 000033
ここで、ｙω,tは入力された観測時間周波数成分Ｙの複素スペクトログラムを表す。k ∈ i は、i番目の音源（楽器）に割り当てられた基底インデックスの集合を示す。
【０１１９】
ＩＦＦＴ適用部３８は、Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部３６によって計算された音源ｉの各々の分離信号の複素スペクトログラム^siω,tに対して、逆フーリエ変換することにより、音源ｉの各々の分離信号を作成する。
【０１２０】
分離信号描画部４０は、ＩＦＦＴ適用部３８によって作成された音源ｉの各々の分離信号を描画した結果を、出力部９０により出力する。
【０１２１】
分離信号再生部４２は、ＩＦＦＴ適用部３８によって作成された音源ｉの各々の分離信号を再生して、出力部９０により出力する。
【０１２２】
＜本発明の第１の実施の形態に係る音源分離装置の作用＞
次に、本発明の第１の実施の形態に係る音源分離装置１００の作用について説明する。まず、入力部１０において複数の音源ｉの各々の分離信号の時系列データを受け付け、分離信号記憶部２２に記憶する。次に、入力部１０において、混合信号の時系列データを受け付け、混合信号記憶部２４に記憶する。そして、音源分離装置１００は、図６に示す音源分離処理ルーチンを実行する。
【０１２３】
まず、ステップＳ１００では、分離信号記憶部２２に記憶されている複数の音源ｉの各々の分離信号の時系列データに基づいて、音源ｉの各々について、観測時間周波数成分Ｓiに変換する。
【０１２４】
次に、ステップＳ１０２では、混合信号記憶部２４に記憶されている混合信号の時系列データに基づいて、観測時間周波数成分Ｙに変換する。
【０１２５】
そして、ステップＳ１０４において、音源ｉについて、ステップＳ１００において取得した観測時間周波数成分Ｓiに基づいて、基底スペクトルＨi、アクティベーションＵi、周期性を表す変数α、βを推定する。
【０１２６】
ステップＳ１０４では、上述したように、音源ｉについて、基底スペクトルＨi、アクティベーションＵi、及び周期性を表す変数α、βの更新と、補助変数λの更新とを繰り返すことにより、基底スペクトルＨi、アクティベーションＵi、周期性を表す変数α、βが推定される。
【０１２７】
ステップＳ１０６では、全ての音源について、上記ステップＳ１０４の処理が実行されたか否かを判定する。上記ステップＳ１０４の処理が実行されていない音源が存在する場合には、上記ステップＳ１０４へ戻り、当該音源について、ステップＳ１０４の処理が実行する。一方、全ての音源について、上記ステップＳ１０４の処理が実行された場合には、ステップＳ１０８へ進む。
【０１２８】
ステップＳ１０８では、上記ステップＳ１０４で音源ｉの各々について推定された基底スペクトルＨiを統合することにより、基底スペクトルＨを生成する。
【０１２９】
そして、ステップＳ１１０では、上記ステップＳ１０２において取得した観測時間周波数成分Ｙと、上記ステップＳ１０８で生成された基底スペクトルＨとに基づいて、アクティベーションＵ、及び周期性を表す変数α、βを推定する。
【０１３０】
ステップＳ１１０では、上述したように、アクティベーションＵ、及び周期性を表す変数α、βの更新と、補助変数λの更新とを繰り返すことにより、アクティベーションＵ、周期性を表す変数α、βが推定される。
【０１３１】
そして、ステップＳ１１２では、上記ステップＳ１０２において取得した観測時間周波数成分Ｙと、上記ステップＳ１０８で生成された基底スペクトルＨと、上記ステップＳ１１０で推定されたアクティベーションＵとに基づいて、音源ｉの各々の分離信号の複素スペクトログラム^siω,tを計算する。
【０１３２】
次のステップＳ１１４では、上記ステップＳ１１２で計算された音源ｉの各々の分離信号の複素スペクトログラム^siω,tに対して、逆フーリエ変換することにより、音源ｉの各々の分離信号を作成する。
【０１３３】
ステップＳ１１６では、上記ステップＳ１１４で作成された音源ｉの各々の分離信号を描画した結果を、出力部９０により出力する。
【０１３４】
ステップＳ１１８では、上記ステップＳ１１４作成された音源ｉの各々の分離信号を再生して、出力部９０により出力し、音源分離処理ルーチンを終了する。
【０１３５】
以上説明したように、第１の実施の形態に係る音源分離装置によれば、基底スペクトルＨと、アクティベーションＵとの積Ｘと、観測時間周波数成分との誤差、各基底ｋのパワーＵk.tに周期性を仮定したときの、アクティベーションＵとの誤差、及び周期関数のスパース化のためのスパース化項を用いて表わされた目的関数を最適化するように、基底スペクトルＨとアクティベーションＵを推定することにより、周期的な発音パターンを抽出することができる。
【０１３６】
また、ＮＭＦにおいてどの音源がいつ発音するかを示すアクティベーションに任意の周期性を仮定し、その周期を推定しながら行列分解を行うことにより、周期的な発音パターンの抽出を促進できる。この際、音楽音響信号の特性を考慮し、周期性を最小限の周期とその倍音成分にあたる正弦曲数の和で表現することにより、音源ごとに異なる任意の周期性を自動推定することができる。
【０１３７】
また、周期的リズムが含まれる楽曲に対する、教師ありの音源分離性能の向上が期待される。また、周期を最小限の周期関数とその倍音にあたる周期関数の和で表現することにより、周期性形状の把握が可能になり、自動採譜の効率化や音楽におけるリズム構造の理解促進が期待される。
【０１３８】
また、一部の音源にのみ周期性を仮定することにより、周期成分と非周期成分の分離が可能になる。
【０１３９】
＜本発明の第２の実施の形態に係る音源分離装置の構成＞
次に、本発明の第２の実施の形態に係る音源分離装置の構成について説明する。なお、第１の実施の形態と同様の構成となる部分については、同一符号を付して説明を省略する。
【０１４０】
第２の実施の形態では、教師データとしての分離信号が入力されずに、教師なしの音源分離を行う点が、第１の実施の形態と主に異なっている。
【０１４１】
図７に示すように、本発明の第２の実施の形態に係る音源分離装置２００は、入力部１０と、演算部２２０と、出力部９０と、を含んで構成されている。
【０１４２】
入力部１０は、複数の音源からの音源信号が混合された混合信号の時系列データを受け付け、混合信号記憶部２４に記憶する。
【０１４３】
演算部２２０は、混合信号記憶部２４と、混合信号入力受付部３０と、混合信号時間周波数変換部３２と、教師なし音源分離適用部２３４と、Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部３６と、ＩＦＦＴ適用部３８と、分離信号描画部４０と、分離信号再生部４２とを含んで構成されている。
【０１４４】
教師なし音源分離適用部２３４は、図８に示すように、入力受付部２５０、混合信号音源分離適用部２５６、及び結果出力部５８を含んで構成されている。
【０１４５】
入力受付部５０は、混合信号時間周波数変換部３２によって出力された観測時間周波数成分Ｙを受け付ける。
【０１４６】
混合信号音源分離適用部２５６は、観測時間周波数成分Ｙに基づいて、各基底ｋの、各周波数ωにおけるパワースペクトルＨω,kを表す基底スペクトルＨと、各基底ｋの、各時刻ｔにおけるパワーＵk.tを表すアクティベーションＵとの積Ｘと、観測時間周波数成分Ｙとの誤差、各基底ｋのパワーＵk.tに周期性を仮定したときの、前記アクティベーションＵとの誤差、及び各基底ｋのパワーＵk.tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、同じ基本周波数で表わされる前記周期関数のグループのスパース化のためのスパース化項を用いて表わされた、上記（１）式で示される目的関数を最適化するように、基底スペクトルＨと、アクティベーションＵ、周期性を表す変数α、βを推定する。
【０１４７】
混合信号音源分離適用部２５６は、上記第１の実施の形態の分離信号音源分離適用部５２と同様に、初期値設定部６０、変数更新部６２、補助変数更新部６４、及び収束判定部６６を含んで構成されている。
【０１４８】
混合信号音源分離適用部２５６の初期値設定部６０は、基底スペクトルＨと、アクティベーションＵと、周期性を表す変数α、βと、補助変数λとに初期値を設定する。
【０１４９】
混合信号音源分離適用部２５６の変数更新部６２は、観測時間周波数成分Ｙと、初期値設定部６０により初期値が設定された、又は変数更新部６２及び補助変数更新部６４により前回更新された、基底スペクトルＨと、アクティベーションＵと、周期性を表す変数α、βと、補助変数λとに基づいて、上記（１８）式、（１９）式、（２２）式、（２３）式に従って、基底スペクトルＨと、アクティベーションＵと、周期性を表す変数α、βとを更新する。
【０１５０】
混合信号音源分離適用部２５６の補助変数更新部６４は、変数更新部６２により更新された、基底スペクトルＨと、アクティベーションＵとに基づいて、上記（１１）式に従って、補助変数λを更新する。
【０１５１】
混合信号音源分離適用部２５６の収束判定部６６は、予め定められた収束条件を満たすまで、変数更新部６２における更新処理と、補助変数更新部６４における更新処理とを繰り返させる。
【０１５２】
＜本発明の第２の実施の形態に係る音源分離装置の作用＞
次に、本発明の第２の実施の形態に係る音源分離装置２００の作用について説明する。まず、入力部１０において、混合信号の時系列データを受け付け、混合信号記憶部２４に記憶する。そして、音源分離装置２００は、図９に示す音源分離処理ルーチンを実行する。なお、第１の実施の形態と同様の処理については、同一符号を付して詳細な説明を省略する。
【０１５３】
まず、ステップＳ１０２では、混合信号記憶部２４に記憶されている混合信号の時系列データに基づいて、観測時間周波数成分Ｙに変換する。
【０１５４】
そして、ステップＳ２１０において、ステップＳ１０２において取得した観測時間周波数成分Ｙに基づいて、基底スペクトルＨ、アクティベーションＵ、周期性を表す変数α、βを推定する。
【０１５５】
ステップＳ２１０では、上述したように、基底スペクトルＨ、アクティベーションＵ、及び周期性を表す変数α、βの更新と、補助変数λの更新とを繰り返すことにより、基底スペクトルＨ、アクティベーションＵ、周期性を表す変数α、βが推定される。
【０１５６】
そして、ステップＳ１１２では、上記ステップＳ１０２において取得した観測時間周波数成分Ｙと、上記ステップＳ２１０で推定された基底スペクトルＨ、及びアクティベーションＵとに基づいて、音源ｉの各々の分離信号の複素スペクトログラム^siω,tを計算する。
【０１５７】
次のステップＳ１１４では、上記ステップＳ１１２で計算された音源ｉの各々の分離信号の複素スペクトログラム^siω,tに対して、逆フーリエ変換することにより、音源ｉの各々の分離信号を作成する。
【０１５８】
ステップＳ１１６では、上記ステップＳ１１４で作成された音源ｉの各々の分離信号を描画した結果を、出力部９０により出力する。
【０１５９】
ステップＳ１１８では、上記ステップＳ１１４作成された音源ｉの各々の分離信号を再生して、出力部９０により出力し、音源分離処理ルーチンを終了する。
【０１６０】
なお、第２の実施の形態に係る音源分離装置２００の他の構成及び作用については、第１の実施の形態と同様であるため、説明を省略する。
【０１６１】
このように、第２の実施の形態に係る音源分離装置によれば、周期的リズムが含まれる楽曲に対する、教師なしの音源分離性能の向上が期待される。
【０１６２】
＜実施例＞
第１、第２の実施の形態における音源分離装置１００、２００による教師なし音源分離、教師あり音源分離の実施例をそれぞれ示す。周波数Ｐmの定義としては式(５)のものを用い、M=20とした。音楽データセットとしては、SiSEC2015(https://sisec.inria.fr/professionally-produced-music-recordings/)で”Professionally-produced music recordings (MUS) ”として公開されている”The Mixing Secret Dataset 100 (MSD100) ” を用いた。
【０１６３】
データセットは異なるジャンルを含む100 のヴォーカル付き楽曲で構成されている。混合信号に加え、歌手または楽器毎に分離された分離信号も公開されている。
【０１６４】
実験データとして、下記の３楽曲のデータを選出した。
【０１６５】
・ musicA - ”AM Contra - Heart Peripheral (Electronic Dance Pop) ”
・ musicB - ”Actions - Devil's Words (Power Pop) ”
・ musicC - ”Actions - South of the Water (Power Pop) ”
【０１６６】
上記の3 楽曲の分離信号は、全てベース、ドラム、その他の楽器、ヴォーカルの4 トラックに分離されていた。
【０１６７】
本実験では、musicA1、musicA2、musicB、musicC として、それぞれ上記楽曲の一部を抽出した。抽出区間は、順番に1:50-2:00, 0:50-0:60, 1:47-1:57, 1:07-1:17 である。各抽出区間は、大きな旋律の変化を避けて、人手で決定した。音楽データはモノラルでサンプリング周波数は22.050kHz とした。短時間フーリエ変換におけるフレーム長は32ms として、フレームシフト長は16ms とした。
【０１６８】
NPCAにおけるＨ，Ｕの初期値には、通常のNMFによる分解結果を使用した。
【０１６９】
NMF におけるＨ，Ｕの初期値には乱数を使用した。正則化パラメータw1、w2はともに0.2 とし、p = 1.0 とした。
【０１７０】
まず、教師なし音源分離において、提案手法NPCA が分離信号を使用することなく周期的に発生する楽音を抽出できることを検証する。図１０はmusicA1 に対するアクティベーションＵk,t の推定結果を示す。基底数Kは24とした。図１０中下部の太枠内の12 個の基底にのみ、周期性制約を適用した。以下、この基底を周期基底と呼ぶ。α、β の推定は、周期基底のみで行った。図１０中で各基底のアクティベーションの左に付けられた楽器ラベルは、分離信号を人手で聴いて付加した。周期的に発音される音と、非周期的に発音される音を分離して、主に周期基底において周期的に発音される音が抽出されていることが確認できる。周期的なリズムが想定されるドラムやベースの音は主に周期基底に、ランダムなリズムが想定されるボーカルは主にそれ以外の基底に割り当てられている。その他の楽器についても周期成分と非周期成分がそれぞれ周期基底、それ以外の基底に分類されていることが分かる。
【０１７１】
図１１に、同じデータに対して通常のNMF を適用した場合の結果を示す。多くの基底のアクティベーションが、ドラムやベースの周期的なリズムに影響され、周期的でないその他の楽器などにおいても、ドラムやベースのリズムに該当するアクティベーションの微小な増加がみられる。アクティベーションの左側に示された楽器ラベルにおいても、多くの基底がベース（11 基底）とドラム（8 基底）に割り当てられており、ヴォーカルは部分的に1 基底に抽出されているのみである。これに対して、提案手法NPCA では、周期成分と非周期成分が分離され、4 基底でヴォーカルの大部分を表現できており、分離性能もNMFと比べて向上していることが確認できた。実験を通して、提案モデルが周期的にアクティベートされる基底を限られた基底数の中にまとめて扱うことができることが示され、この特長は精度と分離信号の解釈性とのトレードオフを解決するのに有効であることが示唆される。加えて、提案モデルを使用すれば、例えば、上記図２のような形でα とβ の値を参照することにより、各基底にどのような周期と位相を持つ周期性が推定されたのかを把握することもできる。このようなリズム構造の把握は、自動採譜や楽曲理解のための分析において重要である。
【０１７２】
続いて、教師あり音源分離においてNPCAの性能を通常のNMFと比較する。
【０１７３】
ここで、NPCA・NMF ともに全ての楽器iに対して基底数Ki= 6 とした。したがって、NPCAにおける、４楽器による混合信号の分離時の基底数Kは24となる。
【０１７４】
定量的に分離性能の比較を行うために、signal-to-noise ratio (SNR) を評価指標として用いる。SNR は
【０１７５】
【数３０】
0006535611.tif 000034
で定義される。ここで、siω,t は入力された楽器i の分離信号の複素スペクトログラムを、^ siω,tは式(24) でWinnerfilter により得られる分離信号の複素スペクトログラムを示す。SNR が高いほど、分解性能が高いことを示す。以下、本実施例では、SNR improvement、すなわちSNR とSNR-all の差を評価のために用いる。SNR-all は^siω,tの代わりに混合信号の複素スペクトログラムを用いて算出したSNR の値である。表1 はNMF とNPCA による音源分離結果から算出したSNR improvement の値の一覧である。
【０１７６】
【表１】
0006535611.tif 000035
【０１７７】
結果から、特にベースやドラムなどの周期性が想定される楽器における分離性能の向上が確認できる。musicC においては、ヴォーカルパートにおいてもSNR improvement の値が向上している。実際、musicC のヴォーカルの旋律は抽出した10 秒間で2 回繰り返されており、提案法が多くの周期的なリズムを奏でる楽器に対して分離性能を向上できることが分かる。楽器によってはSNR improvement の値が低下している場合があるものの、周期性という強い制約を与えているにも限らずtotal のSNR improvement の値がNMF よりも増加していることが確認でき、教師あり音源分離の枠組みにおいても提案手法が音源分離性能を向上できることが示せた。
【０１７８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１７９】
例えば、上記の実施の形態においては、目的関数のスパース化項として、gruop sparsityを誘導するための項を用いる場合を例に説明したが、これに限定されるものではなく、目的関数のスパース化項として、各基底ｋのパワーＵk.tに対して、複数の基本周波数の各々の倍音で表わされる各周期関数の和を用いて周期性を仮定したときの、周期関数のスパース化のためのスパース化項を用いてもよい。具体的には、目的関数に含まれるＬ(α, β)として、上記（９）式で表わされるsparsityを誘導するための項を用いてもよい。
【０１８０】
また、フーリエ変換を用いて、観測時間周波数成分を求める場合を例に説明したが、これに限定されるものではなく、ウェーブレット変換を用いて、時間周波数展開を行って、観測時間周波数成分を求めるようにしてもよい。
【０１８１】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。
【符号の説明】
【０１８２】
１０   入力部
２０、２２０  演算部
２２   分離信号記憶部
２４   混合信号記憶部
２６   分離信号入力受付部
３０   混合信号入力受付部
２８   分離信号時間周波数変換部
３２   混合信号時間周波数変換部
３４   教師あり音源分離適用部
３６   Ｗｉｎｎｅｒｆｉｌｔｅｒ適用部
３８   ＩＦＦＴ適用部
４０   分離信号描画部
４２   分離信号再生部
５０、２５０  入力受付部
５２   分離信号音源分離適用部
５４   基底スペクトル統合部
５６、２５６  混合信号音源分離適用部
５８   結果出力部
６０   初期値設定部
６２   変数更新部
６４   補助変数更新部
６６   収束判定部
９０   出力部
１００、２００       音源分離装置
２３４ 教師なし音源分離適用部
【図１】
0006535611.tif 000036
【図２】
0006535611.tif 000037
【図３】
0006535611.tif 000038
【図４】
0006535611.tif 000039
【図５】
0006535611.tif 000040
【図６】
0006535611.tif 000041
【図７】
0006535611.tif 000042
【図８】
0006535611.tif 000043
【図９】
0006535611.tif 000044
【図１０】
0006535611.tif 000045
【図１１】
0006535611.tif 000046