この記事は前編の続きです。まだ前編を読んでない方はそちらを先に読んでください。

3.　2.で得た基底行列を用いて、かえるの歌のスペクトログラムに対してNMFを適用

今度は、かえるの歌のスペクトログラムに対してNMFを適用します。ただ、2.で得た $F$ を使いたいので、分解の方法を先ほどとは少し変えます。

かえるの歌のスペクトログラムを $Y$ とした時、これを以下のように分解します。 $Y \simeq FG' + HU$

$F\ (M \times K)$ →2.で得たトランペットの教師音の基底行列

$G'\ (K \times N)$ → $F$ に対応するアクティベーション行列

$H\ (M \times L)$ →トランペット以外の音の基底行列

$U\ (L \times N)$ → $H$ に対応するアクティベーション行列

$F$ は2.で求めたものをそのまま使うので、ここでは $G'$ と $H$ と $U$ を変化させながら $FG'+HU$ を $Y$ に近づけます。 $Y$ を上手く分解できれば、 $FG'$ はトランペットのスペクトログラム、 $HU$ はトランペット以外の音(つまりピアノ）のスペクトログラムとなります。

じゃあ、さっきと同じように目的関数を $D=||Y-(FG'+HU)||_F$ と設定して……と行きたいところですが1つ問題があります。このまま分解しようとすると、 $F$ と $H$ の間に何も制約が無いので上手くいかない可能性があります。

上手くいかない例を1つ挙げると、 $Y$ を分解した時に $G'=O$ （ $O$ は零行列）となってしまう可能性があります。そうなると、 $FG'$ も $O$ になってしまうので、トランペットの音を抜き出したら結果が無音なんていう変な結果になってしまいます。

それを防ぐために、先ほどのNMFの目的関数を少し変えて以下のようにします。

$D=||Y-(FG'+HU)||_F + \mu ||F^TH||_F$

$\mu ||F^TH||_F$ という項が増えました。 $||F^TH||_F$ は直感的には、 $F$ と $H$ がどれくらい似ているかを表します。 $F$ と $H$ が似ているほど値が大きくなり、違う行列であるほど値が小さくなるという感じです。

$Y$ を分解した結果が $G'=O$ となったと仮定します。この時、 $Y$ は $HU$ に分解され、 $||Y-(FG'+HU)||_F$ の値はとても小さくなります。しかし、その場合、 $Y$ は $F$ の周波数パターンを含む音源なので、 $H$ は $F$ と似た値を含む行列になるはずです。なので、 $||F^TH||_F$ の値がとても大きくなってしまい、結果として $D$ の値が大きくなります。以上より、 $\mu ||F^TH||_F$ の項によって $G'=O$ となるような分解は避けられます。

ちなみに、 $\mu$ は罰則項の重みを表すパラメータです。

目的関数に罰則項が増えたので、更新式も変更する必要があります。 $G'$ と $H$ と $U$ の更新式は以下のようになります。導出はβ-divergence規範による罰則条件付き教師あり非負値行列因子分解を用いた目的楽器音抽出を見てください。

$G'_{k,t}=G'_{k,t} \displaystyle \frac{\sum_{\omega} F_{\omega, k} Y_{\omega, t}}{\sum_{\omega} F_{\omega, k} Z_{\omega, t}}$

$H_{\omega, l}=H_{\omega, l} \displaystyle \frac{\sum_t Y_{\omega, t} U_{l,t}}{\sum_t U_{l,t} Z_{\omega, t} + 2 \mu \sum_k F_{\omega, k} \sum_{\omega'} F_{\omega', k} H_{\omega', l}}$

$U_{l,t}=U_{l,t} \displaystyle \frac{\sum_{\omega} H_{\omega, l} Y_{\omega, t}}{\sum_{\omega} H_{\omega, l} Z_{\omega, t}}$

ただし、 $Z_{\omega, t}=\sum_k F_{\omega, k} G'_{k, t} + \sum_l H_{\omega, l} U_{l, t}$

コードはこんな感じです。さっきのNMFの更新式がちょっと難しくなっただけなので大きな違いは無いです。

def PSNMF(Y, F, l, num_iter, mu):
    #D = Y -(FG + HU) + mu * norm(F.T, H)
    M, N = Y.shape
    k = F.shape[1]

    G = np.random.rand(k, N)
    U = np.random.rand(l, N)
    H = np.random.rand(M, l)

    error = []
    eps = np.spacing(1)
    H = normalize(H)

    for i in tqdm(range(num_iter)):
        error.append(euclid_norm(Y, np.dot(F, G) + np.dot(H, U)))

        H *= np.dot(Y, U.T) / (np.dot(np.dot(F, G) + np.dot(H, U), U.T) +  2 * mu * np.dot(F, np.dot(F.T, H)) + eps)
        H = normalize(H)
        U *= np.dot(H.T, Y) / (np.dot(H.T, np.dot(F, G) + np.dot(H, U)) + eps)
        G *= np.dot(F.T, Y) / (np.dot(F.T, np.dot(F, G) + np.dot(H, U)) + eps)

    return (G, H, U, error)


def normalize(X):
    N = X.shape[1]

    for n in range(N):
        X[:, n] /= sum(X[:, n])

    return X

途中でHを正規化してるのは、参考にしたスライドにそう書いてあったからです。ぶっちゃけ入れるべき理由はまだよく分かってないです。

4.　3.の結果に逆フーリエ変換をして分離した音源を得る

2.と3.で得られた $F$ と $G'$ をかけた $FG'$ がトランペット音のみのスペクトログラムです。スペクトログラムは周波数の情報しかないので、これを逆フーリエ変換して音の波形に戻します。要するに、1.でやった作業の逆をやればいいだけです。

def istft(x, win, step):
    N, M = x.shape
    l = (M - 1) * step + N
    X = np.zeros(l, dtype = "float64")
    X_count = np.zeros(l, dtype = "int")

    for m in xrange(M):
        start = m * step
        X[start:start+N] += np.fft.ifft(x[:, m]).real / win
        X_count[start:start+N] += 1

    X /= X_count

    return X

結果

さっそく、上のやり方で抜き出した結果を見てみます。

実際に抜き出した音を聴く前に、まずはNMFで分解したスペクトログラムを確認してみます。つまり、手法の中の $FG'$ と $HU$ ですね。

f:id:phoro3:20161022214156p:plain

左がもとのかえるの歌、真ん中が $FG'$ （トランペット）、右が $HU$ （ピアノ）のスペクトログラムです。この図からだと上手くいってるかどうか分かり辛いですね。

ただ、最初の方の時刻（0～2000のあたり）を見ると人の目でもわかる違いがあります。もとの音源ではトランペットとピアノのスタートが1小節ずれています。なので、最初の方はピアノの音は鳴っていません。それを踏まえて先ほどの図を見ると、右側のピアノのスペクトログラムの最初の方の時刻は周波数の強度が小さいです。これはNMFによってトランペットとピアノが分離てきているからだと期待できます。

では、実際に抜き出したトランペットの音を聴いてみましょう。