今回の目的

今回の目的はタイトルの通り、楽曲中から特定の楽器の音のみを抜き出すことです。

具体的に言うと、まずここにトランペットとピアノで演奏したかえるの歌があります。

もう1つ、トランペットでドレミファソラシドの音を演奏した音源があります。（以降、これを教師音と呼びます）

この2つの音源を使って、かえるの歌からトランペットの音のみを抜き出すことが目的です。

導入

突然ですが、楽器の音色の違いってどこから来るかご存知ですか？試しに、トランペットとピアノでラの音（440Hz）を鳴らした時の波形を出してみます。

f:id:phoro3:20161022153702p:plain

左がピアノの波形、右がトランペットの波形です。トランペットの方がピアノに比べて細かく波打っています。楽器によって波形がだいぶ違うというのが分かります。

この違いをプログラムに処理させるには、何らかの数式や数値で表す必要があります。では、この波形の違いを数式で表すにはどうすればいいでしょうか？「トランペットの方が細かく波打っている」というあいまいな説明では数式として扱えません。

そこで、この2つの波形をフーリエ変換してこの2つの楽器の周波数を見てみましょう。

f:id:phoro3:20161022153701p:plain

この図は、音に含まれている周波数の強さを示したグラフです。横軸が周波数（Hz）、縦軸がその周波数の強さを表します。

例えば、左のピアノの図を見ると、鳴らした音であるラの音（440Hz）にピークが来ているのが分かります。それ以外に、1000Hzの手前（880Hz）も値が大きくなっています。このように、楽器の音は鳴らした音だけでなく、それ以外の周波数も含んでいます。

それを踏まえて、右のトランペットの図を見ると、ピアノとは大きく違いますね。ピアノの時は、440Hzの時の値が一番大きいですが、トランペットでは880Hzの値が一番大きいです。このように、同じ音を鳴らした時でも楽器によって含む周波数が大きく違うことが分かります。

つまり、含んでいる周波数の割合を見れば楽器ごとの違いを判別できる、という予想が立てられます。この考え方をもとに、トランペットの音の抜き出しを考えていきます。

やり方

今回用いるのは、教師ありNMFを用いた音源分離という方法になります。

手法の大まかな流れはこんな感じです。

かえるの歌と教師音をフーリエ変換してスペクトログラムを求める
教師音のスペクトログラムに対してNMFを適用して基底行列を得る
2.で得た基底行列を用いて、かえるの歌のスペクトログラムに対してNMFを適用
3.の結果に逆フーリエ変換をして分離した音源を得る

スペクトログラムやNMFといった聞き慣れない単語が出てきてるので、流れだけだと何をやっているか分かりにくいですね。 1つ1つ説明していきたいと思います。

記事中でもコードを示しながら説明しますが、全体のコードを見たいという方はこちらへどうぞ。

1. かえるの歌と教師音をフーリエ変換してスペクトログラムを求める

かえるの歌からトランペットの音を抜き出すために、まずはもとのかえるの歌をプロットしてみましょう。

f:id:phoro3:20161022115325p:plain

やはり、波形だけを見ても何もヒントは得られなさそうですね。なので、まずはこれをフーリエ変換しちゃいます。ただ、全ての区間に対してフーリエ変換をすると、色々な周波数が混ざってよく分からなくなるので短時間フーリエ変換（STFT）を使います。短時間フーリエ変換は、下図の赤い区間のように音を短く区切りながらフーリエ変換をしていく方法です。

f:id:phoro3:20161022132656p:plain

コードはこんな感じです。

def stft(x, win, step):
    N = len(win)
    M = int(np.ceil(float(len(x) - N + step) / step))

    new_x = np.zeros(N + (M - 1) * step, dtype = "float64")
    new_x[:len(x)] = x
    X = np.zeros([N, M], dtype = "complex64")
    for m in xrange(M):
        start = step * m
        X[:, m] = np.fft.fft(new_x[start:start + N] * win)

    return X

コードでは音声データをNごとに区切ってフーリエ変換をしています。コード中の変数winは、窓関数を表しています。

窓関数については、ここで説明すると本題からどんどん離れていきそうなのでここでは書きません。以下のサイトに丁寧な説明が載ってるので、気になる方は見てみてください。

この短時間フーリエ変換をかえるの歌と教師音の両方に行います。教師音に対して短時間フーリエ変換した結果を図示してみます。

f:id:phoro3:20161022133728p:plain

横軸が時間、縦軸が周波数（音の高さ）、色が周波数の強度を表しています（赤いほど強度が強い）。教師音はドレミファソラシドとだんだん音が上がっていく音源でした。なので、図を見ると、時間が経つにつれて周波数が上がっていっています。（あまりキレイには上がってませんが・・・（笑））

この図をスペクトログラムと呼びます。スペクトログラムを見ることで、各時間ごとの周波数の強さが分かります。

2. 教師音のスペクトログラムに対してNMFを適用して基底行列を得る

NMFとは

まずは、NMFの説明をします。 NMFを知っている方はこの部分は飛ばして大丈夫です。

NMFはNon-negative Matrix Factorizationの略で、ある行列 $X\ (M \times N)$ を2つの非負の行列 $F\ (M \times K)$ と $G\ (K \times N)$ の積に分解する手法です。

つまり、 ${ \displaystyle X \simeq FG }$ となるような $F$ と $G$ を求めるのが目的です。

アルゴリズムは割と単純です。 $X$ と $FG$ がなるべく同じ行列になればいいので目的関数 $D$ を以下のように設定します。（ $||\cdot||_F$ はフロベニウスノルムを表します）

${\displaystyle D(X|FG) = ||X-FG||_F \ = \sum_{\omega , t} ( Y_{\omega,t} - \sum_{k}H_{\omega, k}U_{k, t})^{2} }$

あとは目的関数 $D$ を $F$ と $G$ で微分すると、以下の更新式を導出できます。

$F_{\omega,k} = F_{\omega,k} \frac{\sum_{t} X_{\omega,t} G_{k,t}}{\sum_{t} G_{k, t}\sum_{k'}F_{\omega, k'}G_{k', t}}, \hspace{3mm} G_{k,t} = G_{k,t}\frac{\sum_{\omega} X_{\omega,t} F_{\omega,k}}{\sum_{\omega} F_{\omega, k}\sum_{k'} F_{\omega, k'} G_{k', t}}$

この更新式に基づいて、 $F$ と $G$ を繰り返し更新していけば、 $FG$ が $X$ に近づいていきます。実は導出がちょっと面倒くさいんですが、これも説明し始めると長いので他の方にお任せします（笑）。以下のサイトが分かりやすかったです。

NMFアルゴリズムの導出（ユークリッド距離版） - LESS IS MORE

数式通り実装するだけですが、一応コードも示しておきます。

def NMF(X, k, num_iter):
    M, N = X.shape
    H = np.random.rand(M, k)
    U = np.random.rand(k, N)


    error = []
    eps = np.spacing(1)
    for i in tqdm(range(num_iter)):
        Y = np.dot(H, U)
        error.append(euclid_norm(X, Y))

        H *= np.dot(X, U.T) / (np.dot(H, np.dot(U, U.T)) + eps)
        U *= np.dot(H.T, X) / (np.dot(H.T, np.dot(H, U)) + eps)

    return (H, U, error)


def euclid_norm(X, Y):
    d = (X ** 2 + Y ** 2 - 2 * X * Y).sum()
    return d

NMFを教師音のスペクトログラムに適用

NMFについて頑張って説明しましたが、ここからが本題です（笑）

このNMFをスペクトログラムに適用するとどうなるか？を表したのが次の図です。

f:id:phoro3:20161022164853p:plain

(図は直交化及び距離最大化則条件を用いた教師あり非負値行列因子分解による音楽信号分離の6ページより引用)

図のようにスペクトログラムにNMFを適用すると、 $F$ は頻出する周波数のパターン、 $G$ は周波数の時間的な変化を表します。つまり、今回の教師音のように単一の楽器のスペクトログラムにNMFを適用すると、 $F$ はその楽器の周波数のパターンとなります。導入の部分で言った「含んでいる周波数の割合」がこの $F$ にあたります。

これを使えば、楽器の音を上手く抜き出せそうです。

ちなみに、 $F$ は基底行列、 $G$ はアクティベーション行列と呼びます。

後編に続きます。

SHIROBAKO大好き人間のブログ

SHIROBAKOが好きなエンジニアによる技術ブログ

楽曲から特定の楽器のみを抜き出す　前編

今回の目的

導入

やり方

1. かえるの歌と教師音をフーリエ変換してスペクトログラムを求める

2. 教師音のスペクトログラムに対してNMFを適用して基底行列を得る

NMFとは

NMFを教師音のスペクトログラムに適用