音源分離したい対象楽器が決まっている場合、その楽器による基底周波数成分を予め学習した上で NMF を実行できれば、分離性能が高まると考えられる。 このアイデアで NMF を拡張したものが半教師ありNMF (Semi-Supervised NMF : SSNMF) である。
合成音源の振幅スペクトル行列を$Y$と表すと、通常の NMF は周波数基底行列$H$とその周波数のアクティベーション行列$U$の積で$Y$を近似を求める。
$Y \approx HU$
$Y$と$HU$を近似していくことは、$Y$と$HU$の距離を最小化する$H$および$U$を求める最適化問題である。
これに対し Semi-Supervised NMF は予め学習した基底周波数成分$F$を用いて、次の近似を求める。
$Y \approx FG + HU$
ここで$F$は所与の値なので、$Y$と$FG + HU$の距離を最小化する周波数基底行列$H$およびアクティベーション行列$U$に加え、教師基底におけるアクティベーション行列$G$を決定する最適化問題である。
これらのアルゴリズムの説明や最適化問題における目的関数 (コスト関数) を導出している文献として、以下を参照しながら、実装を試みた。
import numpy as np
np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
まずは基本の NMF を、意図した通りに動作するよう実装した。近似度の評価にはユークリッド距離を用いた。
def NMF(Y, R=3, n_iter=50, init_H=[], init_U=[], verbose=False):
"""
decompose non-negative matrix to components and activation with NMF
Y ≈ HU
Y ∈ R (m, n)
H ∈ R (m, k)
HU ∈ R (k, n)
parameters
----
Y: target matrix to decompose
R: number of bases to decompose
n_iter: number for executing objective function to optimize
init_H: initial value of H matrix. default value is random matrix
init_U: initial value of U matrix. default value is random matrix
return
----
Array of:
0: matrix of H
1: matrix of U
2: array of cost transition
"""
eps = np.spacing(1)
# size of input spectrogram
M = Y.shape[0]
N = Y.shape[1]
# initialization
if len(init_U):
U = init_U
R = init_U.shape[0]
else:
U = np.random.rand(R,N);
if len(init_H):
H = init_H;
R = init_H.shape[1]
else:
H = np.random.rand(M,R)
# array to save the value of the euclid divergence
cost = np.zeros(n_iter)
# computation of Lambda (estimate of Y)
Lambda = np.dot(H, U)
# iterative computation
for i in range(n_iter):
# compute euclid divergence
cost[i] = euclid_divergence(Y, Lambda)
# update H
H *= np.dot(Y, U.T) / (np.dot(np.dot(H, U), U.T) + eps)
# update U
U *= np.dot(H.T, Y) / (np.dot(np.dot(H.T, H), U) + eps)
# recomputation of Lambda
Lambda = np.dot(H, U)
return [H, U, cost]
def euclid_divergence(Y, Yh):
d = 1 / 2 * (Y ** 2 + Yh ** 2 - 2 * Y * Yh).sum()
return d
# 動作確認。Y と HU がある程度近似的に見えれば OK
np.random.seed(1)
comps = np.array(((1,0), (0,1), (1,1)))
activs = np.array(((0,0,1,0,1,5,0,7,9,6,5,0), (2,1,0,1,1,2,1,0,0,0,6,0)))
Y = np.dot(comps, activs)
print('original data\n---------------')
print('components:\n', comps)
print('activations:\n', activs)
print('Y:\n', Y)
computed = NMF(Y, R=2)
print('\ndecomposed\n---------------')
print('H:\n', computed[0])
print('U:\n', computed[1])
print('HU:\n', np.dot(computed[0], computed[1]))
print('cost:\n', computed[2])
NMF が実装できたら、次は学習済み基底成分を固定化して NMF を行う Semi-Supervised NMF を実装する。
def SSNMF(Y, R=3, n_iter=50, F=[], init_G=[], init_H=[], init_U=[], verbose=False):
"""
decompose non-negative matrix to components and activation with Semi-Supervised NMF
Y ≈ FG + HU
Y ∈ R (m, n)
F ∈ R (m, x)
G ∈ R (x, n)
H ∈ R (m, k)
U ∈ R (k, n)
parameters
----
Y: target matrix to decompose
R: number of bases to decompose
n_iter: number for executing objective function to optimize
F: matrix as supervised base components
init_W: initial value of W matrix. default value is random matrix
init_H: initial value of W matrix. default value is random matrix
return
----
Array of:
0: matrix of F
1: matrix of G
2: matrix of H
3: matrix of U
4: array of cost transition
"""
eps = np.spacing(1)
# size of input spectrogram
M = Y.shape[0];
N = Y.shape[1];
X = F.shape[1]
# initialization
if len(init_G):
G = init_G
X = init_G.shape[1]
else:
G = np.random.rand(X, N)
if len(init_U):
U = init_U
R = init_U.shape[0]
else:
U = np.random.rand(R, N)
if len(init_H):
H = init_H;
R = init_H.shape[1]
else:
H = np.random.rand(M, R)
# array to save the value of the euclid divergence
cost = np.zeros(n_iter)
# computation of Lambda (estimate of Y)
Lambda = np.dot(F, G) + np.dot(H, U)
# iterative computation
for it in range(n_iter):
# compute euclid divergence
cost[it] = euclid_divergence(Y, Lambda + eps)
# update of H
H *= (np.dot(Y, U.T) + eps) / (np.dot(np.dot(H, U) + np.dot(F, G), U.T) + eps)
# update of U
U *= (np.dot(H.T, Y) + eps) / (np.dot(H.T, np.dot(H, U) + np.dot(F, G)) + eps)
# update of G
G *= (np.dot(F.T, Y) + eps)[np.arange(G.shape[0])] / (np.dot(F.T, np.dot(H, U) + np.dot(F, G)) + eps)
# recomputation of Lambda (estimate of V)
Lambda = np.dot(H, U) + np.dot(F, G)
return [F, G, H, U, cost]
def euclid_divergence(V, Vh):
d = 1 / 2 * (V ** 2 + Vh ** 2 - 2 * V * Vh).sum()
return d
# 動作確認。Y と FG + HU がある程度近似的に見えれば OK
# このパラメータでは近似度がまだ低く、コスト関数の値をみてもまだまだ収束していない感がある
# SSNMF のパラメータ引数 n_iter を大きくすると計算量が増えるが、近似度が高くなる。
np.random.seed(1)
comps = np.array(((1,0), (0,1), (1,1)))
activs = np.array(((0,0,1,0,1,5,0,7,9,6,5,0), (2,1,0,1,1,2,1,0,0,0,6,0)))
Y = np.dot(comps, activs)
print('original data\n---------------')
print('components:\n', comps)
print('activations:\n', activs)
print('Y:\n', Y)
computed = SSNMF(Y, F=np.array(((1,0,1),)).T, R=2)
print('\ndecomposed\n---------------')
print('F:\n', computed[0])
print('G:\n', computed[1])
print('H:\n', computed[2])
print('U:\n', computed[3])
print('FG + HU:\n', np.dot(computed[0], computed[1]) + np.dot(computed[2], computed[3]))
print('cost:\n', computed[4])
ドラム音のみの音声ファイルと、ドラム音+ベース音が合成された音声ファイルがある。
Semi-Supervised NMF を用いてこれらの音源を分離する。
import librosa
import matplotlib.pyplot as plt
import scipy.io.wavfile as wav
from IPython.display import display, Audio
%matplotlib inline
# load wav
y_db, sr_db = librosa.load('data/drums+bass.wav')
y_d, sr_d = librosa.load('data/drums.wav')
y_b, sr_b = librosa.load('data/bass.wav')
plt.subplot(311)
plt.title('mixed')
plt.plot(y_db)
plt.subplot(312)
plt.title('drums only')
plt.plot(y_d)
plt.subplot(313)
plt.title('bass only')
plt.plot(y_b)
plt.tight_layout()
print('ソース音源: ドラム音+ベース音')
display(Audio(y_db, rate=sr_db))
print('ソース音源: ドラム音')
display(Audio(y_d, rate=sr_d))
print('ソース音源: ベース音 (後の計算に使用しないが参考のため)')
display(Audio(y_b, rate=sr_b))
# STFT
S_db = librosa.stft(y_db)
S_d = librosa.stft(y_d)
S_b = librosa.stft(y_b)
# learn drums as base
# ドラム音学習で分解する基底数
R_d = 50
# ドラム音学習で行う反復計算回数
n_iter=300
nmf_d = NMF(np.abs(S_d), R=R_d, n_iter=n_iter)
base = nmf_d[0]
# decompose by SSNMF
# Semi-Supervised NMF で分解する基底数
R = 20
# Semi-Supervised NMF で行う反復計算回数
n_iter = 50
%time ssnmf_db_d = SSNMF(np.abs(S_db), F=base, R=R, n_iter=n_iter)
# confirm each {F / G / H / U} 's size and divergence
for i in range(5):
print(i, np.shape(ssnmf_db_d[i]))
print('cost', ssnmf_db_d[4])
# remix with separated matrix by SSNMF
remixed_only_d = librosa.istft((np.dot(ssnmf_db_d[0], ssnmf_db_d[1]) * (np.cos(np.angle(S_db) + 1j * np.sin(np.angle(S_db))))))
remixed_without_d = librosa.istft(np.dot(ssnmf_db_d[2], ssnmf_db_d[3]) * (np.cos(np.angle(S_db) + 1j * np.sin(np.angle(S_db)))))
print('ドラム音のみの音源から学習した基底を用いて、合成音源からドラム音のみを再合成したもの')
display(Audio(remixed_only_d, rate=sr_db))
print('ドラム音のみの音源から学習した基底を用いて、合成音源からドラム音を除いて再合成したもの')
display(Audio(remixed_without_d, rate=sr_db))
plt.subplot(211)
plt.title('remixed only drums')
plt.plot(remixed_only_d)
plt.subplot(212)
plt.title('remixed without drums')
plt.plot(remixed_without_d)
plt.tight_layout()
NMF は反復的な計算が必要であり精度を高めるためにはそれなりの実行時間が必要になる。 リアルタイム適用を考慮すると、反復計算の回数や、基底数などを調整する必要がある。
基底数については、少ない基底数でもとの振幅スペクトルを復元できれば、分離性能および計算量のどちらに対しても理想的であるが、実際のところ、基底数が少なすぎると復元される音源の品質が悪くなるため、適切な値を見つけるためには試行錯誤が必要である。
学習する基底数と Semi-Supervised NMF で求める基底数の比率は分離結果に大きく影響を与える。学習する基底数のほうが大きいと分離されにくく、逆に小さいと分散されやすい傾向があるように感じられる。