(57)【要約】【課題】各構成音が混合した観測信号から、各構成音を精度よく分離することができる。【解決手段】パラメータ推定部３６が、各構成音が混合された観測信号の時系列データを入力として、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、基底スペクトル及びアクティベーションパラメータの積との誤差の大きさを表す規準を小さくするように、各構成音の信号と、基底スペクトルと、アクティベーションパラメータとを推定する。このとき、パラメータ更新部４４が、補助関数を小さくするように、各構成音の信号と、基底スペクトルと、アクティベーションパラメータとを更新する。各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新する。【選択図】図１
2018072664.tif 000002
【特許請求の範囲】
【請求項１】
各構成音が混合された観測信号の時系列データを入力として、予め学習された各構成音及び各周波数における振幅スペクトルを表す基底スペクトルに基づいて、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、前記基底スペクトル及び各構成音の各時刻における音量を表すアクティベーションパラメータの積との誤差の大きさを表す規準を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを推定するパラメータ推定部
を含む信号解析装置であって、
前記パラメータ推定部は、
前記規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新するパラメータ更新部であって、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新するパラメータ更新部と、
予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させる収束判定部と、
を含む信号解析装置。
【請求項２】
前記パラメータ更新部は、以下の式に従って、各構成音の信号を更新する請求項１記載の信号解析装置。
【数１】
2018072664.tif 000135
2018072664.tif 000136
ただし、ｓlは、構成音ｌの信号を表し、
2018072664.tif 000137
は、前記補助関数の偏微分を表し、γはステップサイズを表し、Ｌは、構成音の数を表し、ｙは、前記観測信号を表す。
【請求項３】
前記規準は、以下の式で表わされる請求項１又は２記載の信号解析装置。
【数２】
2018072664.tif 000138
2018072664.tif 000139
ただし、
2018072664.tif 000140
であり、βl,k,mは、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みを表し、ψk,mは、時刻ｍ、周波数ｋの振幅スペクトルを測るための基底関数を表し、（・）Hは、ベクトルの複素共役転置であり、Ｈk,lは、各周波数ｋ及び各構成音ｌの基底スペクトルを表し、Ｕl,mは、各構成音ｌの各時刻ｍにおけるアクティベーションパラメータを表し、λは、重みパラメータであり、ｐは定数である。
【請求項４】
パラメータ推定部が、各構成音が混合された観測信号の時系列データを入力として、予め学習された各構成音及び各周波数における振幅スペクトルを表す基底スペクトルに基づいて、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、前記基底スペクトル及び各構成音の各時刻における音量を表すアクティベーションパラメータの積との誤差の大きさを表す規準を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを推定する
信号解析方法であって、
前記パラメータ推定部が推定することでは、
パラメータ更新部が、前記規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新するパラメータ更新部であって、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新し、
収束判定部が、予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させることを含む信号解析方法。
【請求項５】
前記パラメータ更新部が更新することでは、以下の式に従って、各構成音の信号を更新する請求項４記載の信号解析方法。
【数３】
2018072664.tif 000141
2018072664.tif 000142
ただし、ｓlは、構成音ｌの信号を表し、
2018072664.tif 000143
は、前記補助関数の偏微分を表し、γはステップサイズを表し、Ｌは、構成音の数を表し、ｙは、前記観測信号を表す。
【請求項６】
前記規準は、以下の式で表わされる請求項４又は５記載の信号解析方法。
【数４】
2018072664.tif 000144
2018072664.tif 000145
ただし、
2018072664.tif 000146
であり、βl,k,mは、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みを表し、ψk,mは、時刻ｍ、周波数ｋの振幅スペクトルを測るための基底関数を表し、（・）Hは、ベクトルの複素共役転置であり、Ｈk,lは、各周波数ｋ及び各構成音ｌの基底スペクトルを表し、Ｕl,mは、各構成音ｌの各時刻ｍにおけるアクティベーションパラメータを表し、λは、重みパラメータであり、ｐは定数である。
【請求項７】
請求項１〜請求項３の何れか１項に記載の信号解析装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号解析装置、方法、及びプログラムに係り、特に、パラメータを推定する信号解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
近年、非負値行列因子分解(Non-negative Matrix Factorization; NMF)に基づく音源分離手法が注目されている。観測信号の振幅（あるいはパワー）スペクトログラムを非負値行列と見なし、これを二つの非負値行列の積で近似する。これは観測スペクトログラムをランク1のスペクトログラムの和で近似していると考えられる。このアプローチではスペクトルの加法性を仮定しているが、振幅スペクトルやパワースペクトルは実際には加法的ではない。この問題を解決するため複素NMFと呼ばれる枠組みが提案されている。複素NMFではNMFにおける各スペクトルテンプレートに位相スペクトルを時変パラメータとして付加したモデルを用いている。同様の動機からParry et al.とFevotte et al.はそれぞれ、観測信号の短時間フーリエ変換(Short-Time Fourier Transform; STFT)から得られた複素スペクトログラムに対し各構成音のパワースペクトログラムがランク１の構造をもち、位相スペクトログラムが一様に分布した潜在変数とみなした生成モデルを提案している。各構成音の複素スペクトログラムの各周波数成分が平均0の複素正規分布に従うとみなすと、モデルパラメータの最尤推定は観測パワースペクトログラムを板倉斎藤擬距離規準でフィッティングしたNMFとなることが示されている。これは板倉斎藤NMFと呼ばれている。同様に複素正規分布の代わりに複素コーシー分布を生成モデルとして用いたものも提案されている。
【０００３】
短時間フーリエ変換(Short-Time Fourier Transform; STFT)やウェーブレット変換などの時間周波数変換によって得られる複素スペクトログラムは変換元の信号の冗長表現であるため、複素数の要素をもつ任意の二次元配列が時間領域の信号に対応する複素スペクトログラムであるためには、各要素間で満たすべき制約条件が存在する。例えばSTFTではある幅ごとに短時間フレームの信号を切り出し、それぞれに対してフーリエ変換を行ったものを連結することで複素スペクトログラムが得られる。従って、各フレームの複素スペクトルを逆フーリエ変換した際にオーバーラップした区間内の波形が無矛盾となる条件を満たす必要がある。複素NMFや板倉斎藤NMFが位相を考慮したNMFとして提案されているが、位相スペクトログラムの各要素が独立なパラメータ（もしくは潜在変数）として扱われており、各時間周波数成分の位相が満たすべき制約条件が考慮されていなかった。このため、推定された複素スペクトログラムが必ずしも時間領域信号に対応したものにならないという問題があった。
【０００４】
この課題を解決することを目的とし、発明者らは以前NMFと同様のコンセプトに基づく信号分解を時間領域で実現する手法「時間領域低ランクスペクトログラム近似法(Time-domain Spectrogram Factorization; TSF)」と呼ぶ手法を提案している（非特許文献２）。この枠組の特徴は、分解成分の振幅スペクトログラムができるだけランク1構造に近くなるように観測信号を時間領域で分解することが可能である点にある。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】亀岡弘和, 小野順貴, 柏野邦夫, 嵯峨山茂樹, "複素NMF: 新しいスパース信号分解表現と基底系学習アルゴリズム," 日本音響学会2008 年秋季研究発表会講演論文集, 2-8-13, pp. 657-660, Sep. 2008.
【非特許文献２】Hirokazu Kameoka, "Multi-resolution signal decomposition with time-domain spectrogram factorization," in Proc. 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP2015), pp. 86{90, Apr. 2015.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかし、この最適化問題を解決するために以前提案したパラメータ推定アルゴリズムには大きな行列の逆行列計算が含まれており、計算コストに解決すべき課題が残されていた。
【０００７】
本発明では、上記事情を鑑みて成されたものであり、計算コストを抑えて、各構成音が混合した観測信号から、各構成音を精度よく分離することができる信号解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、本発明に係る信号解析装置は、各構成音が混合された観測信号の時系列データを入力として、予め学習された各構成音及び各周波数における振幅スペクトルを表す基底スペクトルに基づいて、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、前記基底スペクトル及び各構成音の各時刻における音量を表すアクティベーションパラメータの積との誤差の大きさを表す規準を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを推定するパラメータ推定部を含む信号解析装置であって、前記パラメータ推定部は、前記規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新するパラメータ更新部であって、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新するパラメータ更新部と、予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させる収束判定部と、を含んで構成されている。
【０００９】
本発明に係る信号解析方法は、パラメータ推定部が、各構成音が混合された観測信号の時系列データを入力として、予め学習された各構成音及び各周波数における振幅スペクトルを表す基底スペクトルに基づいて、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、前記基底スペクトル及び各構成音の各時刻における音量を表すアクティベーションパラメータの積との誤差の大きさを表す規準を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを推定する信号解析方法であって、前記パラメータ推定部が推定することでは、パラメータ更新部が、前記規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新するパラメータ更新部であって、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新し、収束判定部が、予め定められた収束条件を満たすまで、前記パラメータ更新部による更新を繰り返させることを含む。
【００１０】
また、本発明のプログラムは、コンピュータを、上記の信号解析装置を構成する各部として機能させるためのプログラムである。
【発明の効果】
【００１１】
以上説明したように、本発明の信号解析装置、方法、及びプログラムによれば、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、前記基底スペクトル及び前記アクティベーションパラメータの積との誤差の大きさを表す規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新することにより、計算コストを抑えて、各構成音が混合した観測信号から、各構成音を精度よく分離することができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の実施の形態に係る信号解析装置の機能的構成を示すブロック図である。
【図２】本発明の実施の形態に係る信号解析装置における学習処理ルーチンを示すフローチャート図である。
【図３】本発明の実施の形態に係る信号解析装置におけるパラメータ推定処理ルーチンを示すフローチャート図である。
【図４】実験結果を示す図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１４】
＜問題の定式化＞
＜目的関数＞
観測信号を
2018072664.tif 000003
とする。TSFではＬ個の信号成分
2018072664.tif 000004
が
【００１５】
【数１】
2018072664.tif 000005
【００１６】
となる条件の下で、各成分の振幅スペクトログラムができるだけランク１構造に近くなるように
2018072664.tif 000006
を求める手法である。ここで、
2018072664.tif 000007
を時刻ｍ、周波数ｋの時間周波数成分を測るための基底関数とする。
【００１７】
2018072664.tif 000008
の振幅スペクトログラムは、
2018072664.tif 000009
を用いて
2018072664.tif 000010
と表される。ただし、(・)Hはベクトルの複素共役転置である。TSFは、目的関数
【００１８】
【数２】
2018072664.tif 000011
2018072664.tif 000012
【００１９】
を最小化する最適化問題として定式化される。ただし、
2018072664.tif 000013
が推定すべきパラメータであり、βl,k,mはβl,k,m＞０、Σlβl,k,m=1を満たすものとする。この問題はβl,k,m の逆数を重みとした重み付き二乗誤差最小化問題として捉えることができる。この重みの重要性については後述する。
【００２０】
2018072664.tif 000014
の一つ目の項は振幅スペクトログラム
2018072664.tif 000015
が完全にランク１構造を有したとき0になる。複素NMFと同様、このモデルでは要素ごとに値を相殺することを許容するため
2018072664.tif 000016
をスパースに誘導するための正則化項が必要になる。
【００２１】
2018072664.tif 000017
の二つ目の項はこのために導入した
2018072664.tif 000018
ノルムであり、λ＞０は正則化項の強さを決める重みパラメータである。
【００２２】
2018072664.tif 000019
が有界であり、0＜ｐ＜２のとき
2018072664.tif 000020
ノルムは
2018072664.tif 000021
をスパースに誘導できる。
【００２３】
2018072664.tif 000022
を有界にするため、
【００２４】
【数３】
2018072664.tif 000023
【００２５】
という制約を設ける。
【００２６】
＜補助関数法＞
パラメータ
2018072664.tif 000024
の目的関数を
2018072664.tif 000025
とするとき
【００２７】
【数４】
2018072664.tif 000026
【００２８】
を満たす
2018072664.tif 000027
を
2018072664.tif 000028
の補助関数、
2018072664.tif 000029
を補助変数と呼ぶ。式(4)を満たすような補助関数を設計できれば、パラメータ
2018072664.tif 000030
と補助変数
2018072664.tif 000031
を
【００２９】
【数５】
2018072664.tif 000032
【００３０】
により反復的に更新することで目的関数
2018072664.tif 000033
を局所最小化する
2018072664.tif 000034
を得ることができる。従来のNMFの反復アルゴリズムは上述の原理に基づいて導かれる。
【００３１】
＜Ｊの補助関数＞
目的関数
2018072664.tif 000035
は
2018072664.tif 000036
が微分不可能な点を含んでいるため、目的関数を局所最小化する
2018072664.tif 000037
の更新式は解析的に求まらない。しかし次節に記載する不等式(10)(11) を用いることで
2018072664.tif 000038
の補助関数
【００３２】
【数６】
2018072664.tif 000039
【００３３】
が得られる。ここで
2018072664.tif 000040
である。補助変数
2018072664.tif 000041
についての更新式は次の式で得られる。
【００３４】
【数７】
2018072664.tif 000042
2018072664.tif 000043
【００３５】
＜不等式＞
ここでは補助関数を設計するために用いた不等式をまとめる。
【００３６】
（Lemma 1.）任意の複素数ｚ、｜ｃ｜=1を満たす複素数ｃについて次式が成り立つ。
【００３７】
【数８】
2018072664.tif 000044
【００３８】
等号はｃ＝ｚ／｜ｚ｜のとき成立する。
【００３９】
（Lemma 2.）0＜ｐ＜2のとき、任意の実数（または複素数）ｘについて次式が成り立つ。
【００４０】
【数９】
2018072664.tif 000045
【００４１】
等号はｖ＝ｘのとき成立する。
【００４２】
＜射影勾配法を用いた最適化＞
補助関数が最小となる
2018072664.tif 000046
は、ラグランジェの未定乗数法を用いることで
【００４３】
【数１０】
2018072664.tif 000047
【００４４】
が得られる。ただし
【００４５】
【数１１】
2018072664.tif 000048
2018072664.tif 000049
2018072664.tif 000050
【００４６】
であることが示される。式(12)が
2018072664.tif 000051
の更新式として用いられる場合、大きな行列の逆行列を計算する必要があるため計算コストが課題となる。そこで、上記非特許文献２では
2018072664.tif 000052
において
2018072664.tif 000053
が単位行列となることに着目し、逆行列計算が不要な方法が提案されているが、βを固定することが信号分解性能に限界を与えることが実験的に明らかになっている。そこで本発明の実施の形態では
2018072664.tif 000054
の更新に関して射影勾配法を用いることで、βを自由に更新できるようにしつつ逆行列計算を避けることが可能な効率的なアルゴリズムを提案する。
【００４７】
＜ｓlに関する更新式＞
ここでは
2018072664.tif 000055
に関する更新のみ考えるため、
2018072664.tif 000056
の他の変数は固定して考える。そのため、
2018072664.tif 000057
を
2018072664.tif 000058
の関数とし、
2018072664.tif 000059
とする。このとき
2018072664.tif 000060
の偏微分は次の式で得られる。
【００４８】
【数１２】
2018072664.tif 000061
【００４９】
式(16)の各項は効率的に計算可能であり、特に
2018072664.tif 000062
がSTFTの基底関数の状況について考える。このとき、
2018072664.tif 000063
は
2018072664.tif 000064
をSTFTした複素スペクトログラムの(k，m)成分である。そのため式(13)の
2018072664.tif 000065
は逆STFTを計算することに対応しており、式(16)の初項は
2018072664.tif 000066
の逆STFTを計算することで得られる。同様に式(14)の
2018072664.tif 000067
は
2018072664.tif 000068
の逆STFTを計算することで得られる。以上のことから
2018072664.tif 000069
の更新は勾配法を用いて効率的に計算可能であり、
【００５０】
【数１３】
2018072664.tif 000070
【００５１】
という更新を行えばよい。ただし
2018072664.tif 000071
はステップサイズである。ここで
2018072664.tif 000072
は勾配
2018072664.tif 000073
のリプシッツ定数であり、行列
2018072664.tif 000074
の最大固有値から得られる。これを効率的に計算するため、次の式で得られるγをステップサイズとして用いる。
【００５２】
【数１４】
2018072664.tif 000075
2018072664.tif 000076
【００５３】
ただし
2018072664.tif 000077
はSTFTの窓関数、αmはオーバーラップの割合を示すホップサイズである。
【００５４】
次に
【００５５】
【数１５】
2018072664.tif 000078
2018072664.tif 000079
【００５６】
を定義する。ここで
2018072664.tif 000080
はＮ×Ｎの単位行列である。これを用いて線形制約
2018072664.tif 000081
は
2018072664.tif 000082
と表記される。また、
【００５７】
【数１６】
2018072664.tif 000083
【００５８】
を定義する。このときアフィン空間
2018072664.tif 000084
への射影は
【００５９】
【数１７】
2018072664.tif 000085
【００６０】
となる。さらに
【００６１】
【数１８】
2018072664.tif 000086
2018072664.tif 000087
【００６２】
であることを用いると
【００６３】
【数１９】
2018072664.tif 000088
【００６４】
が得られる。更新式(17)、(25) を交互に繰り返すことにより、
2018072664.tif 000089
を
2018072664.tif 000090
について最小化する更新が得られた。実際はρの上界が
2018072664.tif 000091
であるため式(18)、(19)で得られるγは
2018072664.tif 000092
を満たす。従って射影勾配法による更新式は、線形制約の下で
2018072664.tif 000093
を最小化する値に収束することが保証される。
【００６５】
＜提案手法全体の更新式＞
2018072664.tif 000094
に関して
2018072664.tif 000095
の偏微分を0と置くことで更新式
【００６６】
【数２０】
2018072664.tif 000096
【００６７】
が得られる。また
2018072664.tif 000097
についてはラグランジェの未定乗数法を用いることで
【００６８】
【数２１】
2018072664.tif 000098
2018072664.tif 000099
【００６９】
が得られる。
【００７０】
以上より、全体のアルゴリズムは下記のようにまとめることができる。
【００７１】
Step1)
2018072664.tif 000100
の初期値を設定する。
【００７２】
Step2)  式（８）、（９）に従って補助変数
2018072664.tif 000101
を更新する。
【００７３】
Step3)  式（１６）、（１８）、（１９）、（２５）に従って、構成音の信号
2018072664.tif 000102
と重みγとを更新する。
【００７４】
具体的には、繰り返し回数ｔ＝１から、以下の式に従って、繰り返し回数ｔにおける各構成音ｌの信号
2018072664.tif 000103
を繰り返し更新する。
【００７５】
【数２２】
2018072664.tif 000104
【００７６】
Step4)  式（２６）〜（２８）に従って、モデルパラメータ
2018072664.tif 000105
を更新する。
【００７７】
＜複素ＮＭＦとの関係＞
実は、上述した補助関数
2018072664.tif 000106
は複素NMFで得られる補助関数と似たものであることがわかる。
【００７８】
複素NMFでは観測複素スペクトログラムＹk,mを次のモデルを用いて近似することが目的である。
【００７９】
【数２３】
2018072664.tif 000107
【００８０】
ここでφl,k,mはl番目の音源の位相スペクトログラムである。先行研究ではＨk,l， Ｕl,m，φl,k,mを変数とする目的関数
【００８１】
【数２４】
2018072664.tif 000108
【００８２】
を最小化する問題として扱われている。このとき
2018072664.tif 000109
である。また
【００８３】
【数２５】
2018072664.tif 000110
【００８４】
が
2018072664.tif 000111
の補助関数になっている。ただし
2018072664.tif 000112
が補助変数であり、
2018072664.tif 000113
である。またβl,k,mはΣlβl,k,m = 1 を満たす正の数であり、Ｘl,k,mはＹk,m=ΣlＸl,k,mを満たす。このときＸl,k,mが推定されたl番目の音源の複素スペクトログラムの推定値である。
【００８５】
式(7) と式(31)を比べると，複素NMFの補助関数におけるＸl,k,mと
2018072664.tif 000114
がTSFの補助関数における
2018072664.tif 000115
とｃl,k,m に対応していることがわかる。複素NMFの一つの問題はＸl,k,mが複素スペクトログラムが満たすべき条件について考慮できていないという点である。つまり複素NMF アルゴリズムでは広い空間で解を探索することになってしまっている。一方、提案手法ではｌ番目の複素スペクトログラム
2018072664.tif 000116
を手がかりに、適切な解空間の中で時間領域信号
2018072664.tif 000117
を推定することができる。
【００８６】
＜本発明の実施の形態に係る信号解析装置の構成＞
次に、本発明の実施の形態に係る信号解析装置の構成について説明する。図１に示すように、本発明の実施の形態に係る信号解析装置１００は、ＣＰＵと、ＲＡＭと、後述する学習処理ルーチン及びパラメータ推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することができる。この信号解析装置１００は、機能的には図１に示すように入力部１０と、演算部２０と、出力部９０と、を含んで構成されている。
【００８７】
入力部１０は、各構成音の構成音信号の時系列データを受け付ける。また、入力部１０は、複数の構成音が混じっている音響信号（以後、観測信号）の時系列データを受け付ける。
【００８８】
演算部２０は、時間周波数展開部２４と、基底スペクトル学習部２８と、基底スペクトル記憶部３０と、パラメータ推定部３６と、を含んで構成されている。
【００８９】
時間周波数展開部２４は、各構成音信号の時系列データに基づいて、各時刻における各周波数のスペクトルを表す振幅スペクトログラム又はパワースペクトログラムを計算する。なお、本実施の形態においては、短時間フーリエ変換やウェーブレット変換などの時間周波数展開を行う。
【００９０】
基底スペクトル学習部２８は、時間周波数展開部２４によって計算された、各構成音信号の各時刻における各周波数のスペクトルに基づいて、従来技術であるＮＭＦを用いて、各構成音ｌ及び各周波数kにおける振幅スペクトルを表す基底スペクトルＨk,lを推定する。
【００９１】
基底スペクトル記憶部３０は、基底スペクトル学習部２８によって推定された、各構成音ｌ及び各周波数kにおける振幅スペクトルを表す基底スペクトルＨk,lを記憶している。
【００９２】
パラメータ推定部３６は、観測信号の時系列データ、及び基底スペクトル記憶部３０に記憶された各構成音及び各周波数における基底スペクトルに基づいて、各構成音
2018072664.tif 000118
の信号から求められる各周波数の各時刻における振幅スペクトルと、基底スペクトルＨk,l及びアクティベーションパラメータＵl,mの積との誤差の大きさを表す上記（２）式の規準を小さくするように、各構成音ｌの信号
2018072664.tif 000119
と、各構成音及び各周波数の基底スペクトルＨk,lと、各構成音の各時刻におけるアクティベーションパラメータＵl,mとを推定する。
【００９３】
具体的には、パラメータ推定部３６は、初期値設定部４０、補助変数更新部４２、パラメータ更新部４４、及び収束判定部４６を備えている。
【００９４】
初期値設定部４０は、各構成音信号の基底スペクトルＨの初期値として、基底スペクトル記憶部３０に記憶された各構成音信号の各周波数における基底スペクトルＨk,lを設定する。また、初期値設定部４０は、各構成音信号のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000120
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに初期値を設定する。
【００９５】
補助変数更新部４２は、初期値である、又は前回更新した、各構成音ｌの信号
2018072664.tif 000121
と、各構成音信号のアクティベーションパラメータＵl,mとに基づいて、上記（８）式〜（９）式に従って、各構成音ｌ、各周波数k、及び各時刻ｍに対するｃl,k,m、各構成音ｌ及び各時刻ｍに対するＶl,mを更新する。
【００９６】
パラメータ更新部４４は、観測信号の時系列データｙと、補助変数更新部４２によって更新された各構成音ｌ、各周波数k、及び各時刻ｍに対するｃl,k,m、各構成音ｌ及び各時刻ｍに対するＶl,mと、初期値である、又は前回更新した、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000122
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに基づいて、上記（７）式に示す補助関数を小さくするように、各構成音ｌの信号
2018072664.tif 000123
と、各構成音及び各周波数の基底スペクトルＨk,lと、各構成音の各時刻におけるアクティベーションパラメータＵl,mと、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mを更新する。
【００９７】
具体的には、パラメータ更新部４４は、観測信号の時系列データｙと、補助変数更新部４２によって更新された各構成音ｌ、各周波数k、及び各時刻ｍに対するｃl,k,m、各構成音ｌ及び各時刻ｍに対するＶl,mと、初期値である、又は前回更新した、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000124
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに基づいて、上記（１６）式、（１８）〜（１９）式、（２５）式に従って、各構成音ｌの信号
2018072664.tif 000125
と、ステップサイズγとを更新する。ここで、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、以下の式に従って、各構成音の信号を繰り返し更新する。
【００９８】
【数２６】
2018072664.tif 000126
【００９９】
ただし、Ｌは、構成音の数を表し、ｔは繰り返し回数を表す。
【０１００】
また、パラメータ更新部４４は、初期値である、又は前回更新した、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000127
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに基づいて、上記（２６）式〜（２８）式に従って、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとを更新する。
【０１０１】
収束判定部４６は、収束条件を満たすか否かを判定し、収束条件を満たすまで、補助変数更新部４２における更新処理と、パラメータ更新部４４における更新処理とを繰り返させる。
【０１０２】
収束条件としては、例えば、繰り返し回数が、上限回数に到達したことを用いることができる。あるいは、収束条件として、上記（２）式の規準の値と前回の規準の値との差分が、予め定められた閾値以下であることを用いることができる。
【０１０３】
出力部９０は、パラメータ推定部３６において取得した各構成音ｌの信号
2018072664.tif 000128
を出力する。
【０１０４】
＜本発明の実施の形態に係る信号解析装置の作用＞
次に、本発明の実施の形態に係る信号解析装置１００の作用について説明する。まず、入力部１０において各構成音の構成音信号の時系列データを受け付けると、信号解析装置１００は、図２に示す学習処理ルーチンを実行する。
【０１０５】
まず、ステップＳ１００では、入力部１０において受け付けた各構成音の構成音信号の時系列データに基づいて、各構成音信号の各時刻における各周波数のスペクトルを計算する。
【０１０６】
次に、ステップＳ１０４では、上記ステップＳ１００において取得した各構成音信号の各時刻における各周波数のスペクトルに基づいて、従来技術であるＮＭＦにより、各構成音ｌ及び各周波数kにおける振幅スペクトルを表す基底スペクトルＨk,lを推定し、基底スペクトル記憶部３０に格納して、学習処理ルーチンを終了する。
【０１０７】
次に、入力部１０において、各構成音が混在した観測信号の時系列データを受け付けると、信号解析装置１００は、図３に示すパラメータ推定処理ルーチンを実行する。
【０１０８】
まず、ステップＳ１２０では、に基づいて、観測スペクトログラムＹを計算する。
【０１０９】
ステップＳ１２２では、各構成音ｌの基底スペクトルＨk,lの初期値として、基底スペクトル記憶部３０に記憶された各構成音及び各周波数における基底スペクトルを設定する。また、各構成音のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000129
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに初期値を設定する。
【０１１０】
ステップＳ１２４では、初期値である、又は後述するステップＳ１２６、Ｓ１２８で前回更新した、各構成音ｌの信号
2018072664.tif 000130
と、各構成音信号のアクティベーションパラメータＵl,mとに基づいて、上記（８）式〜（９）式に従って、各構成音ｌ、各周波数k、及び各時刻ｍに対するｃl,k,m、各構成音ｌ及び各時刻ｍに対するＶl,mを更新する。
【０１１１】
次に、ステップＳ１２６では、入力部１０において受け付けた観測信号の時系列データｙと、上記ステップＳ１２４で更新された各構成音ｌ、各周波数k、及び各時刻ｍに対するｃl,k,m、各構成音ｌ及び各時刻ｍに対するＶl,mと、初期値である、又はステップＳ１２６、Ｓ１２８で前回更新した、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各構成音ｌの信号
2018072664.tif 000131
と、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとに基づいて、上記（１６）式、（１８）〜（１９）式、（２５）式に従って、各構成音ｌの信号
2018072664.tif 000132
と、ステップサイズγとを更新する。ここで、各構成音の信号を更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、以下の式に従って、各構成音の信号を繰り返し更新する。
【０１１２】
【数２７】
2018072664.tif 000133
【０１１３】
ステップＳ１２８では、初期値である、又は前回更新した、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mと、上記ステップＳ１２６で更新した各構成音ｌの信号
2018072664.tif 000134
と、に基づいて、上記（２６）式〜（２８）式に従って、各構成音信号の各周波数における基底スペクトルＨk,lと、各構成音信号のアクティベーションパラメータＵl,mと、各時刻ｍ、各周波数ｋ、及び各構成音ｌの重みβl,k,mとを更新する。
【０１１４】
次に、ステップＳ１３０では、収束条件を満たすか否かを判定する。収束条件を満たした場合には、ステップＳ１３２へ移行し、収束条件を満たしていない場合には、ステップＳ１２４へ移行し、ステップＳ１２４〜ステップＳ１２８の処理を繰り返す。
【０１１５】
ステップＳ１３２では、上記ステップＳ１２６で最終的に更新された各構成音の信号を、出力部９０から出力して、パラメータ推定処理ルーチンを終了する。
【０１１６】
＜実験例＞
本実施の形態の手法の有効性を確認するため、SiSEC 2013 Database の"professionally produced music recordings"（5音源からなるポピュラー音楽の音響信号）を用いて3分割交差検証による教師あり音源分離実験を行い，既存のTSF(β1,k,m= ・・・ =βL,k,m)、複素NMF(CNMF-EUC)、Iダイバージェンス規準NMF(NMF-KL)、二乗誤差規準NMF(NMF-EUC)、板倉斎藤擬距離規準NMF(NMF-IS)による各音源のSignal to Distortion Ratio (SDR) 改善値を比較した。図４にその結果を示す。提案手法が他手法に比べて高い改善値を得られていることが確認できた。
【０１１７】
以上説明したように、本発明の実施の形態に係る信号解析装置によれば、時間領域スペクトログラム分解法のアルゴリズムにおいて、各構成音の信号から求められる各周波数の各時刻における振幅スペクトルと、基底スペクトル及びアクティベーションパラメータの積との誤差の大きさを表す規準の上界関数である補助関数を小さくするように、各構成音の信号と、各構成音及び各周波数の基底スペクトルと、各構成音の各時刻におけるアクティベーションパラメータとを更新する際、各構成音の信号の和が前記観測信号と等しくなる線形制約を満たす空間へ各構成音の信号を射影するように、各構成音の信号を更新することにより、計算コストを抑えて、各構成音が混合した観測信号から、各構成音を精度よく分離することができる。
【０１１８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１１９】
例えば、各構成音信号の基底スペクトルを学習する処理と、観測信号からパラメータを推定するパラメータ推定とを別々の装置で行うように構成してもよい。
【０１２０】
また、更新するパラメータの順番には任意性があるため、上記の実施の形態の順番に限定されない。
【０１２１】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。
【符号の説明】
【０１２２】
１０ 入力部
２０ 演算部
２４ 時間周波数展開部
２８ 基底スペクトル学習部
３０ 基底スペクトル記憶部
３６ パラメータ推定部
４０ 初期値設定部
４２ 補助変数更新部
４４ パラメータ更新部
４６ 収束判定部
９０ 出力部
１００ 信号解析装置
【図１】
2018072664.tif 000147
【図２】
2018072664.tif 000148
【図３】
2018072664.tif 000149
【図４】
2018072664.tif 000150