２９．K-shapeで音源分離してみよう - keep learning blog（キープラーニングブログ）

――俯いていたら、虹を見つけることはできないよ。――

K-shapeで音源分離してみよう

K-shapeで音源分離してみよう

喜劇王チャップリン、本当にいいこと言いますね。COVID-19で依然としてストレスフルな日々は続いているものの、顔を上げれば突き抜けるような青空。私は猫背だから姿勢もよくしないといけないです。

最近、妻とアパートの周りを散歩することを日課にしているのですが、ちょっと前まで枯れ木だった木々にも新緑が芽吹き、池のそばの道路では鴨のお母さんがまだ小さい子供たちを引率しています。もうすぐ、私がデンマークに来てから1年が経とうとしています。

つい先週、スペインのバルセロナでICASSP2020*1という音声処理業界で最大の学会が開催される予定でした。しかし、今年はCOVID-19のためにオンラインによるバーチャル開催になりました。

私はもともとICASSPに参加する予定はなかったものの、バーチャル開催の恩恵を受けて誰でも参加費無料になったので、最新の機械学習関連の講演を見てきました（自宅で）。

いくつか面白い講演があり、その中で有望そうなディープラーニングベースの音源分離アルゴリズムを発見しました。自力で実装できたら（たぶんそんなに難しくない）次回に紹介したいと思います。

その前に、今回はずっと長いこと前フリをしつつ記事にしてこなかった、K-shapeアルゴリズムとそれを使った音源分離をご紹介します。

f:id:yuki0718:20200517220727j:plain

時系列クラスタリングとK-shape

SBD（Shape-based distance）とは

K-shapeとは、2015年にコロンビア大学のJohn Paparrizosらが発表した論文k-Shape: Efficient and Accurate Clustering of Time Seriesで初めて提案された、（音声には全く関係ない）新しいクラスタリング手法です。

以前の記事で紹介したとおり、入力した複数のデータ点をk個の近接グループに分割（Hard clustering）する手法としては、K-meansが古くから用いられてきました（今でも現役）。

keep-learning.hatenablog.jp

K-meansは、入力データと重心とのユークリッド距離（L2ノルム）を最小化するアルゴリズムです。例えば、3次元のデータ $\boldsymbol{x}=(1,-2,5)$ があったとき、重心１ $C_{1}=(1,1,1)$ 、重心２ $C_{2}=(2,2,2)$ とのユークリッド距離は以下のように計算できます。

データと重心１の距離 $=\sqrt{(1-1)^{2}+(-2-1)^{2}+(5-1)^{2}}=\sqrt{25}$

データと重心２の距離 $=\sqrt{(1-2)^{2}+(-2-2)^{2}+(5-2)^{2}}=\sqrt{26}$

データに近いのは重心１の方なので、このデータはグループ１に属すると結論付けるわけです。他方、K-shapeはこの距離の測り方をSBD（Shape-based distance）という、時系列データに特化した距離尺度に置き換えたアルゴリズムです。

時系列データとは、上述したデータの「次元」が、「時間（日時や時刻）」に置き換わったものだと考えてください。つまり、4月1日の株価、4月2日の株価、・・・4月30日の株価、のように並んだデータのことです。

なぜSBDを使う必要があるのかを知るために、以下のようなK-meansでのグループ分けがうまくいかない典型例を見てみましょう。

f:id:yuki0718:20200517225726j:plain

ここに図示した時系列データ１と２は、いわゆるsinカーブとcosカーブを離散的にプロットして結んだものです。つまり、時間的にシフトすれば完全に重なる時系列データです。

しかしながら、K-meansのようなL2ノルムで両データの距離を計算すると、各次元（時間）ごとのデータ差分（黒い矢印で示した長さ）を二乗して足し上げるので、両データ間の距離は非常に大きい値になります。

もう一つのK-meansでのグループ分けがうまくいかない典型例は、以下のケースです。

f:id:yuki0718:20200517175221j:plain

ここで図示した時系列データ２は、時系列データ１の振幅を単に0.25倍したものです。つまり、本質的には同じ時間変動を示している時系列データです。にも関わらず、両データ間の距離はかなり大きい値になります。

以上のように、K-meansで使っているような、三平方の定理で単純計算されたユークリッド距離では、時系列データ同士の距離を適切に測れないことが分かります。距離が適切でなければ当然、重心位置やグループ分けも適切ではありません。

前者のような違いで距離が変化しないことをシフト不変性（Time-shift invariance）、後者のような違いで距離が変化しないことをスケール不変性（Scale invariance）と呼び、SBDはまさにシフト不変かつスケール不変な距離尺度と言えます。

SBDの計算方法

端的に言えば、SBDは、データ１とシフト済みデータ２の相互相関関数（Cross correlation function）を、データ１とデータ２の自己相関関数（Autocorrelation function）で割り算（規格化）した量に関係しています。すなわち、

$\displaystyle SBD(x, y)=1-{\rm max}_{s} \biggl( \frac{R_{s}(x, y)}{\sqrt{R_{0}(x, x)R_{0}(y, y)}} \biggr)$

ただし、 $\displaystyle R_{s}(a, b) \stackrel{\mathrm{def}}{\equiv} \sum_{n=1}^{N} a_{n}b_{n-s}$

maxと書いてあるのは、片方の時系列データを横にs=1点分シフトしたときの相互相関、s=2点分シフトしたときの相互相関、・・・をそれぞれのsについて計算して、相互相関の最大値を採用するからです。

例えば、上述した１つ目の典型例では、橙色のグラフを右方向に6点分くらいシフトすれば青色のグラフと概ね重なりそうです。この最も相関が高いシフト状態の相関値を採用することで、シフト不変な距離になるわけです。

また、両データのシフトなし自己相関関数を分母にしているので、データのスケールが規格化され、２つ目の典型例のようなスケールが違うデータ同士でも大きな差が出にくくなります。

ところで、信号処理の分野では、シフト不変性とスケール不変性を考慮した動的時間伸縮法（DTW：Dynamic time warping）という優れた古典的手法があるのですが、K-shapeはその代替手法として提案されています（論文中でもDTWと同等の性能を目標にしている）。

実は、SBDを使う最大のメリットは、その計算コストの低さです。つまり、DTWと同じ問題意識を持ち同レベルのパフォーマンスを出せるのに、計算時間はSBDの方がDTWよりはるかに短い点です*2。

どうしてそんなにSBDの計算が速いのかというと、相関関数は現代のPCなら離散フーリエ変換（FFT）で超高速に計算できるからです*3。すなわち、

$\displaystyle SBD(x, y)=1-{\rm max}_{s} \biggl( \frac{FFT^{-1}_{s}(FFT(x)FFT^{*}(y))}{\sqrt{FFT^{-1}_{0}(FFT(x)FFT^{*}(x))*FFT^{-1}_{0}(FFT(y)FFT^{*}(y))}} \biggr)$

コンピュータではこの式で計算ができます。このとき、シフト量 $s$ が逆離散フーリエ変換の係数番号に対応している点に注目してください。

PythonのNumpyやScipyライブラリでは、逆離散フーリエ変換が各係数番号をインデックスとする配列データとして出力されます。その配列の中から最大値を抜き出すだけなので、SBDの実装は非常に簡単です。

重心の計算方法

距離尺度をSBDに置き換えると、データと重心との距離をSBDで測り、SBDが最も小さい重心に所属するようにグループ分けできます。しかし、問題はSBDにおける重心 $\boldsymbol{c}_{j}$ の計算方法をどうするかということです。

実はここが原著論文のキモでもありまして、数学的な導出は飛ばして結果だけ論文から抜き出しますと、SBDベースの重心は以下のような最大化問題に帰着します。

$\displaystyle \newcommand{\argmax}{\mathop{\rm arg~max}\limits} \boldsymbol{c}_{j}=\begin{equation} \argmax_{\boldsymbol{c}_{j}} \biggl( \frac{\boldsymbol{c}^{T}_{j}\boldsymbol{M}\boldsymbol{c}_{j}}{\boldsymbol{c}^{T}_{j}\boldsymbol{c}_{j}} \biggr) \end{equation}$

ただし、 $\displaystyle \boldsymbol{M} \stackrel{\mathrm{def}}{\equiv} \biggl( \boldsymbol{I}-\frac{1}{p} \biggr) ^{T} \boldsymbol{X}^{T}\boldsymbol{X} \biggl( \boldsymbol{I}-\frac{1}{p} \biggr)$

$\boldsymbol{X}$ はそのグループに所属するデータを各行に格納した行列であり、 $\boldsymbol{I}$ は単位行列、 $p$ はデータの次元数（時系列データの時間点数）です。

これは線形代数学で最大固有値問題として知られる形式です。つまり、行列 $\boldsymbol{M}$ に関する固有値問題を解き、固有値が最大になるときの固有ベクトル $\boldsymbol{c}_{j}$ が求める重心になります。

ちなみに、行列の固有値問題を解く際には、NumpyやScipyライブラリの「linalg」という線形代数用クラスに便利な関数が用意されているので、こちらもPythonでの実装は簡単です。

Pythonで実装

毎度自己流で恐縮ですが、以上の計算式に基づいて、PythonでK-shapeのアルゴリズムをコードに書き起こしてみました。

import sys
import time
import matplotlib.pyplot as plt
import numpy as np
from scipy import fftpack as fp
from scipy import linalg

#Function for getting a shape-based distance (SBD)
def get_SBD(x, y):
    
    #Define FFT-size based on the length of input
    p = int(x.shape[0])
    FFTlen = int(2**np.ceil(np.log2(2*p-1)))
    
    #Compute the normalized cross-correlation function (NCC)
    CC = fp.ifft(fp.fft(x, FFTlen)*fp.fft(y, FFTlen).conjugate()).real
    
    #Reorder the IFFT result
    CC = np.concatenate((CC[-(p-1):], CC[:p]), axis=0)
    
    #To avoid zero division
    denom = linalg.norm(x) * linalg.norm(y)
    if denom < 1e-10:
        denom = numpy.inf
    NCC = CC / denom
    
    #Search for the argument to maximize NCC
    ndx = np.argmax(NCC, axis=0)
    dist = 1 - NCC[ndx]
    #Get the shift parameter (s=0 if no shift)
    s = ndx - p + 1
    
    #Insert zero padding based on the shift parameter s
    if s > 0:
        y_shift = np.concatenate((np.zeros(s), y[0:-s]), axis=0)
    elif s == 0:
        y_shift = np.copy(y)
    else:
        y_shift = np.concatenate((y[-s:], np.zeros(-s)), axis=0)
    
    return dist, y_shift

#Function for updating k-shape centroid
def shape_extraction(X, v):
    
    #Define the length of input
    N = int(X.shape[0])
    p = int(X.shape[1])
    
    #Construct the phase shifted signal
    Y = np.zeros((N, p))
    for i in range(N):
        #Call my function for getting the SBD between centeroid and data
        _, Y[i, :] = get_SBD(v, X[i, :])
    
    #Construct the matrix M for Rayleigh quotient
    S = Y.T @ Y
    Q = np.eye(p) - np.ones((p, p)) / p
    M = Q.T @ (S @ Q)
    
    #Get the eigenvector corresponding to the maximum eigenvalue
    eigen_val, eigen_vec = linalg.eig(M)
    ndx = np.argmax(eigen_val, axis=0)
    new_v = eigen_vec[:, ndx].real
    
    #The ill-posed problem has both +v and -v as solution
    MSE_plus = np.sum((Y - new_v)**2)
    MSE_minus = np.sum((Y + new_v)**2)
    if MSE_minus < MSE_plus:
        new_v = -1*new_v
    
    return new_v

#Function for checking empty clusters
def check_empty(label, num_clu):
    
    #Get unique label (which must include all number 0~num_clu-1)
    label = np.unique(label)
    
    #Search empty clusters
    emp_ind = []
    for i in range(num_clu):
        if i not in label:
            emp_ind.append(i)
    
    #Output the indices corresponding to the empty clusters
    return emp_ind

#Function for getting KShape clustering
def get_KShape(X, num_clu, max_iter, num_init):
    
    #Define the length of input
    N = int(X.shape[0])  #The number of data
    p = int(X.shape[1])  #The length of temporal axis
    
    #Repeat for each trial (initialization)
    minloss = np.inf
    for init in range(num_init):
        
        #Initialize label, centroid, loss as raondom numbers
        label = np.round((num_clu-1) * np.random.rand(N))
        center = np.random.rand(num_clu, p)
        loss = np.inf
        
        #Normalize the centroid
        center = center - np.average(center, axis=1)[:, np.newaxis]
        center = center / np.std(center, axis=1)[:, np.newaxis]
        
        #Copy the label temporarily
        new_label = np.copy(label)
        new_center = np.copy(center)
        
        #Repeat for each iteration
        for rep in range(max_iter):
            
            #Reset loss value
            new_loss = 0
            
            ### Refinement step (update center) ###
            #Repeat for each cluster
            for j in range(num_clu):
                
                #Construct data matrix for the j-th cluster
                clu_X = []
                for i in range(N):
                    #If the i-th data belongs to the j-th cluster
                    if label[i] == j:
                        clu_X.append(X[i, :])
                clu_X = np.array(clu_X)
                
                #Call my function for updating centroid
                new_center[j,:] = shape_extraction(clu_X, center[j,:])
                
                #Normalize the centroid
                new_center = new_center - np.average(new_center, axis=1)[:, np.newaxis]
                new_center = new_center / np.std(new_center, axis=1)[:, np.newaxis]
            
            ### Assignment step (update label) ###
            #Repeat for each data
            for i in range(N):
                
                #Define the minimum distance
                mindist = np.inf
                
                #Repeat for each cluster
                for j in range(num_clu):
                    
                    #Call my function for getting the shape based distance
                    dist, _ = get_SBD(new_center[j,:], X[i, :])
                    
                    #Assign the label corresponding to the minimum distance
                    if dist < mindist:
                        #Update minimum distance
                        mindist = dist
                        new_label[i] = j
                
                #Get summation of the SBD
                new_loss = new_loss + mindist
            
            ### Error handling (avoid empty clusters) ###
            #Call my function for checking empty clusters
            emp_ind = check_empty(new_label, num_clu)
            if len(emp_ind) > 0:
                for ind in emp_ind:
                    #Assign the same index of data as the one of cluster
                    new_label[ind] = ind
            
            #Get out of the loop if loss and label unchange
            if loss - new_loss < 1e-6 and (new_label == label).all():
                #print("The iteration stopped at {}".format(rep+1))
                break
            
            #Update parameters
            label = np.copy(new_label)
            center = np.copy(new_center)
            loss = np.copy(new_loss)
            #print("Loss value: {:.3f}".format(new_loss))
        
        #Output the result corresponding to minimum loss
        if loss < minloss:
            out_label = np.copy(label).astype(np.int16)
            out_center = np.copy(center)
            minloss = loss
    
    #Output the label vector and centroid matrix
    return out_label, out_center, minloss

アルゴリズムの全体構造としてはK-meansとほぼ同じです。最も近い重心を探して適切なグループ番号 $Label_{i}=1,2,\cdots,k$ を割り当てるAssignment stepと、グループの所属データに基づいて重心の再計算を行うRefinement stepとを交互に繰り返して収束させています。

K-meansと異なるのは、重心を計算する関数「shape_extraction(X, v)」が最大固有値問題になっている点と、グループ番号を重心に割り当てるときの距離尺度が関数「get_SBD(x, y)」で計算されている点です*4。

シフト不変性とスケール不変性を備えた尺度なので、時系列データ同士をグループ分けする際にはかなり強力な手法です*5。

音源分離にK-shapeを応用してみる

前節までK-shapeの基礎理論とPythonでの実装をご紹介しましたが、これは音に関わらず時系列データのクラスタリング一般に適用可能な手法です。

とはいえ、ご存知のとおり音声というのは典型的な時系列データです。今回はK-shapeをうまく使って、世界初*6の音源分離アルゴリズムを考案してみました。

アイデアはシンプルです。以前の記事で、2009年にM. SpiertzとV. Gnannが発表した「Source-Filter Based Clustering For Monaural Blind Source Separation」の実装をご紹介しました。

keep-learning.hatenablog.jp

これは、NMFで分解された基底行列 $H_{k,m}$ の各列ベクトルが、楽器の音色ごとに共通の特徴を持っていると仮定し、MFCCとK-meansで基底番号 $m$ をグループ分けしたのでした。私のアイデアは、NMFで分解されたもう片方の行列 $U_{m,n}$ を利用するものです。

NMFの記事でご説明したとおり、アクティベーション行列 $U_{m,n}$ は、各行ベクトルが音量の時間変化（時系列データ）を表しています。そして、同じ楽器に由来するなら音を出すタイミングもほぼ同じはずなので、 $U_{m,n}$ の各行ベクトルはK-shapeでグループ分けできるのではないでしょうか。

この仮説に基づいてPythonコードを実装*7し、テストしてみた結果がこちらです。今回も、SiSEC2011のDevelopment datasetで提供されている曲「Que_pena_tanto_faz」と「Roads」を分離してみます。

まずは音圧が同一（0dB）となるように「Que_pena_tanto_faz」のギターパートとボーカルパートを混合したものがこちらです。

これを私の提案手法（Activation K-shape）で分離したところ、こうなりました。

予想以上にきれいに分離できています。少なくとも以前記事にしたSpiertzらの手法（MFCC+K-means）より分離性能（SDRスコア）は高いことが分かりました。

続いて、二曲目の「Roads」でギターパートとドラムスパートを混合したものがこちらです。

これを私の提案手法で分離したところ、こうなりました。

こちらはSpiertzらの手法ではほとんど分離できなかった（のでブログにも載せなかった）のですが、提案手法ではまともに分離できています。

シンプルなアイディアではあるものの、なかなか有望そうな音源分離アルゴリズムです。私が情報系の学部生だったら卒論のネタくらいにはなったかもしれません。もしご興味があったら実装してみてください。

今回はここまでです。ちょっと長かったですが、K-shapeは時系列データにも使いやすいですし、何より計算コストが低いので今後もいろいろ使っていきたいアルゴリズムです。

なお、2016年くらいにFateme Fahimanさんという人がFuzzy c-shapeなる拡張版アルゴリズムを提案していたので今回ついでに実装してみました。しかし、こちらは収束性に致命的な欠陥があるような気がするので、さほど有望ではなさそうです（勘違いだったらごめんなさいFatemeさん）。

次回は記念すべき30回目の記事です。別にアクセス数を稼ぐことを目的に書いてるブログではないので*8、30回更新できるかどうかが1年間の（ひそかな）目標でした。実際にやってみると定期的にブログを書くのってとても大変でした。でもいい経験でした。

稚文をお読みいただきありがとうございました。

*1:International Conference on Acoustics, Speech, and Signal Processing

*2:実際、私も自分で実装して試してみたのですが、5倍～10倍くらい計算時間が違いました。

*3:ウィナー・ヒンチンの定理（Wiener–Khinchin theorem）という有名な定理により、自己相関関数はフーリエ変換の絶対値二乗の逆フーリエ変換、相互相関関数はフーリエ変換同士の積（片方は複素共役）の逆フーリエ変換で計算できることが示されています。これを使わず、定義どおりにFor文などで愚直に実装するとめちゃくちゃ時間がかかります。

*4:私が以前の記事に掲載したように、簡単なK-meansの方を自力で実装しておくと全体構造に間違いがなくなるのでおすすめです。

*5:なお、他の論文を読んでいると、DTWに比べて次元（時間点数）が大きい時系列データになるとパフォーマンスが落ちることが報告されています。あまりにも長尺なデータにはDTWの方が向いているかもしれません。

*6:私がすぐに思いつくくらいだから、簡単すぎて誰もやらないだけかもしれないですけど。

*7:ちょっとした工夫ですが、アクティベーション行列に時間方向の滑らかさ（Smoothness）を与えるため、通常のNMFにT. Virtanenらの考案したTemporal continuity costを導入しています。

*8:もし人に見せるためのブログなら、もっと一文一文を短くしないと読んでもらえません。自分が他人としてこのブログに偶然アクセスしたら、一目見てすぐに戻るボタンをクリックする気がします。