半教師あり NMF による音源分離の実装実験

音源分離したい対象楽器が決まっている場合、その楽器による基底周波数成分を予め学習した上で 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$を決定する最適化問題である。

これらのアルゴリズムの説明や最適化問題における目的関数 (コスト関数) を導出している文献として、以下を参照しながら、実装を試みた。

参考

In [1]:
import numpy as np
np.set_printoptions(formatter={'float': '{: 0.2f}'.format})

Basic NMF

まずは基本の NMF を、意図した通りに動作するよう実装した。近似度の評価にはユークリッド距離を用いた。

In [2]:
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
In [3]:
# 動作確認。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])
original data
---------------
components:
 [[1 0]
 [0 1]
 [1 1]]
activations:
 [[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:
 [[ 0  0  1  0  1  5  0  7  9  6  5  0]
 [ 2  1  0  1  1  2  1  0  0  0  6  0]
 [ 2  1  1  1  2  7  1  7  9  6 11  0]]

decomposed
---------------
H:
 [[ 0.17  4.75]
 [ 2.28  0.00]
 [ 2.47  4.65]]
U:
 [[ 0.84  0.42  0.00  0.42  0.44  0.89  0.42  0.07  0.08  0.06  2.64  0.00]
 [ 0.00  0.00  0.21  0.00  0.20  1.03  0.00  1.47  1.89  1.26  0.96  0.00]]
HU:
 [[ 0.14  0.07  1.01  0.07  1.00  5.02  0.07  6.99  8.99  5.99  5.00  0.00]
 [ 1.91  0.96  0.00  0.96  1.00  2.03  0.96  0.16  0.19  0.13  6.02  0.00]
 [ 2.07  1.04  0.99  1.04  2.00  6.97  1.04  7.01  9.00  6.00  10.99  0.00]]
cost:
 [ 260.52  25.89  24.77  24.13  23.55  22.87  21.88  20.40  18.27  15.55
  12.54  9.63  7.11  5.09  3.56  2.46  1.71  1.20  0.87  0.65  0.51  0.41
  0.35  0.30  0.27  0.24  0.22  0.20  0.19  0.18  0.17  0.16  0.15  0.14
  0.13  0.13  0.12  0.12  0.11  0.11  0.10  0.10  0.09  0.09  0.09  0.08
  0.08  0.08  0.07  0.07]

Semi-Supervised NMF

NMF が実装できたら、次は学習済み基底成分を固定化して NMF を行う Semi-Supervised NMF を実装する。

In [4]:
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
In [5]:
# 動作確認。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])
original data
---------------
components:
 [[1 0]
 [0 1]
 [1 1]]
activations:
 [[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:
 [[ 0  0  1  0  1  5  0  7  9  6  5  0]
 [ 2  1  0  1  1  2  1  0  0  0  6  0]
 [ 2  1  1  1  2  7  1  7  9  6 11  0]]

decomposed
---------------
F:
 [[1]
 [0]
 [1]]
G:
 [[ 0.00  0.00  0.00  0.00  0.05  0.12  0.00  1.08  0.85  0.86  0.20  0.00]]
H:
 [[ 3.37  0.71]
 [ 0.00  1.84]
 [ 3.19  2.60]]
U:
 [[ 0.00  0.00  0.29  0.00  0.16  1.22  0.00  1.72  2.37  1.49  0.67  0.00]
 [ 0.83  0.42  0.03  0.42  0.55  1.13  0.42  0.17  0.22  0.14  3.33  0.00]]
FG + HU:
 [[ 0.59  0.30  1.00  0.30  0.98  5.04  0.30  7.00  9.01  6.01  4.84  0.00]
 [ 1.53  0.77  0.05  0.77  1.02  2.08  0.77  0.31  0.40  0.25  6.12  0.00]
 [ 2.17  1.08  1.00  1.08  2.00  6.95  1.08  6.99  8.98  5.99  11.01  0.00]]
cost:
 [ 212.61  27.14  20.97  19.82  18.41  16.70  14.74  12.61  10.49  8.57
  6.99  5.79  4.90  4.26  3.78  3.41  3.11  2.86  2.64  2.45  2.29  2.14
  2.01  1.89  1.79  1.69  1.61  1.53  1.46  1.40  1.34  1.29  1.24  1.20
  1.15  1.11  1.08  1.04  1.01  0.98  0.95  0.92  0.90  0.87  0.85  0.82
  0.80  0.78  0.76  0.73]

音声データへの適用

ドラム音のみの音声ファイルと、ドラム音+ベース音が合成された音声ファイルがある。

Semi-Supervised NMF を用いてこれらの音源を分離する。

  1. ドラム音のみの波形を短時間フーリエ変換して振幅スペクトルを得る。
  2. ドラム音の振幅スペクトルを基本 NMF によって分解し、基底成分を得る。これを学習済み基底成分と見立てる。
  3. ドラム音+ベース音の波形を短時間フーリエ変換して振幅スペクトルを得る。
  4. ドラム音+ベース音の振幅スペクトルを、2 で得た学習済み基底成分とともに Semi-Supervised NMF によって分解する。
  5. 分解した基底成分とアクティベーション行列から、ドラム音のみの振幅スペクトルおよびドラム音を除いた振幅スペクトルを再合する
  6. ドラム音のみの振幅スペクトルおよびドラム音を除いた振幅スペクトルを逆短時間フーリエ変換して波形データを再合成する。
In [6]:
import librosa
import matplotlib.pyplot as plt
import scipy.io.wavfile as wav
from IPython.display import display, Audio
%matplotlib inline
In [7]:
# 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))
ソース音源: ドラム音+ベース音
ソース音源: ドラム音
ソース音源: ベース音 (後の計算に使用しないが参考のため)