８．HPSSで音楽を分離してみよう - keep learning blog（キープラーニングブログ）

――私には、完璧ではない音楽のつくり方がわからないんですよ。――

ヴォルフガング・アマデウス・モーツァルト

HPSSで音楽を分離してみよう
- 音声処理の基本中の基本
- 調波打楽器音分離（HPSS）

HPSSで音楽を分離してみよう

私には、完璧ではない仕事のやり方がわからないんですよ。なんて、一度でいいから言ってみたいですね。天才モーツァルトだから許される一言です。

さて、飽き飽きしてるとは思いますが、今回もプログラミング関連です。最近だんだん音声処理の雰囲気が分かってきて少し面白くなってきました。でもさすがに今回でいったん最後にします（書くのが大変だし読む人も退屈だから）。

keep-learning.hatenablog.jp

音声処理の基本中の基本

第６回の記事でいきなりデジタルフィルタの話を始め、ばりばりフィルタの設計とかしてきたにも関わらず、今さら基本中の基本です。教科書なら順番逆だろうという気もしますが、個人的にここで基礎をまとめておきたかったので書きます。

f:id:yuki0718:20190804205310j:plain

デジタルとアナログ

デジタルとアナログ、という言葉をご存知でしょうか。一番身近な例で言うと、時計が分かりやすいと思います。デジタル時計の場合、「9時5分23秒」と表示されていたら、その表示になってから0.999...秒間はずっと同じまま停止し、1秒きっかりで表示が「9時5分24秒」に切り替わります。一方、秒針が滑らかに動くタイプのアナログ時計であれば、「9時5分23秒～9時5分24秒」の間も秒針は停止することなく動き続けます。

まさに、これがコンピュータと現実世界との違いです。コンピュータは、ビットと呼ばれる記憶装置がたくさん並んでいるデジタルな計算機です。ビットは表か裏の状態を記憶できるコインのようなもので、「ちょっと表」とか「そこそこ裏」みたいな中途半端な状態は存在しません。一方、現実の物理現象では、何もかもが連続的で、白黒はっきりつくことは稀です。

サンプリング

音声処理の世界では、現実の連続的な（Consecutive）音の波動をマイクなどで電気信号に変えたあと、コンピュータに入力します。この一連の入力作業の中で、音はコンピュータのビットに記録できるよう、デジタル化される必要があります。

時間に関して言えば、もちろん秒刻みのデジタル時計ほど粗い単位ではなく、１万分の１秒のような極めて短い時間に区切ってデジタル化され（サンプリング）、音の波形をなぞる様に点を打っていきます。音の大きさを表す波の振幅に関しても、通常 $2^{16}(=65536)$ 段階に細かくデジタル化されます（16bit量子化）。

上述した時間の区切り幅をサンプリング間隔Ts（あんまり使わないけど）と呼び、その逆数をサンプリング周波数Fsと呼びます。サンプリング周波数は時間の逆数なので、「1秒間あたりに何個のデジタルな点を打つか」という点の個数を表す量です。以下にFs=15でサンプリングした例を示します。

f:id:yuki0718:20190804171029j:plain — デジタルサンプリングのイメージ

ちなみに、標準規格でCDのサンプリングレートFsは44100（1秒間に4万点以上！）だそうです。つまり、5分間のJ-POPなら、その中には5×60×44100=13230000点ものデータ点が格納されていることになります。ものすごい量です。

デジタルデータと離散フーリエ変換

音の振幅 $x(t)$ は時間領域（Time Domain）の関数なので、フーリエ変換と呼ばれる数学的な操作を行えば、周波数領域（Frequency Domain）に「移す」ことができ、フーリエ逆変換を行えば、時間領域に「戻す」こともできます。

一見これは正しい主張ですが、実は、サンプリングを行ったデジタルデータの場合には、以下のようにちょっとした仮定を置く必要があります*1。

デジタル化され有限の長さ（N個の点）として記録された音の振幅は、時間 $t$ の代わりにサンプリング番号 $n (=0,1,2,...N-1)$ を使って $x_{n}$ と書けます。ここで、 $x_{n}$ のフーリエ変換 $X(\omega)$ が存在するとしたら、 $x_{n}$ は $X(\omega)$ のフーリエ逆変換、

$\displaystyle{x_{n}=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega nT_{s}}=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega n/F_{s}} }$

で書けるはずです（時間 $t=n \times T_{s}$ であることと $T_{s}$ が $F_{s}$ の逆数であることを用いた）。

しかしながら、フーリエ変換は無限に続く関数を想定した数学的操作であるため、N個の有限長さしかない $x_{n}$ には成り立ちません*2。そこで、この有限長さの複製（コピー）が無限に繰り返すという仮定（周期的境界条件 $x_{n+N}=x_{n}$ ）を置くことにより、フーリエ変換を可能にします。

f:id:yuki0718:20190804165236j:plain — 周期的境界条件のイメージ

$\displaystyle{x_{n+N}=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega (n+N)/F_{s}}=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{i\omega n/F_{s}}=x_{n} }$

この周期的境界条件が任意の $n$ について成り立つためには、

$e^{i\omega N/F_{s}}=1$ 　すなわち、 $\displaystyle{\omega=\frac{2{\pi}F_{s}}{N}k}　(kは整数)$

となる必要があります。 $\omega$ は角周波数ですから、周波数 $F$ とは $\omega=2{\pi}F$ という関係にあり、

$\displaystyle{F=\frac{F_{s}}{N}k}　(kは整数)$

というシンプルな条件になります。これは周波数が連続量ではなく飛び飛びの値である（デジタル化されている）ことを意味します。

以上をまとめると、有限なN個のサンプリングデータをフーリエ変換可能にするためには、周期的境界条件を課す必要があり、その結果としてデータは周波数領域でもデジタル化されてしまうことが分かりました。

そして、このデジタル化された周波数点の間隔 $\displaystyle{\Delta F=\frac{F_{s}}{N} }$ が周波数の分解能（最小単位）になり、フーリエ変換は離散フーリエ変換（DFT：Discrete Fourier Transformation）、

$\displaystyle{X_{k}=\sum_{n=0}^{N-1} x_{n}e^{-i\omega nT_{s}}=\sum_{n=0}^{N-1} x_{n}e^{-i2{\pi}kn/N} }$

によって実行されます。

短時間フーリエ変換（STFT）

現実に音声を分析しようと思ったとき、デジタル記録された音の波動を周波数の重ね合わせに分解（Decompose）できる離散フーリエ変換DFTは、たいへん強力な武器です。そのため、音声処理は常にフーリエ変換のコンピュータ処理能力の向上と共に発展してきました。

しかしながら、分析対象とする音声データというのは、数十秒から数分という長大なデータであり、その間も時々刻々と周波数の分布が変わります。そのため、音声データ全体にわたってフーリエ変換しても、ただの平均値が出てくるだけで、あまり意味のある周波数分布は採れません*3。

そこで音声認識の分野で用いられているのが、短時間フーリエ変換（STFT：Short Time Fourier Transformation）という手法です。

STFTでは、まず長大な音声データから幅N個分のサンプリング点（時間にしてN×Ts）を含む単位で波形を切り取り、それを離散フーリエ変換します。切り取られた単位を時間窓（Time Window）と呼び、切り取り幅N（N×Ts）を窓幅（Window Width）と呼びます。

続いて、時間窓を時間軸方向に少しずつ（普通は窓幅の半分ずつ）ずらしながら、それぞれの時間窓で切り取られた音声データに対して離散フーリエ変換を実行します。以上がSTFTの行われる手順です*4。

このSTFTの出力結果は、当然ながら周波数Fの関数であると同時に、窓が切り取った時刻（時間）の関数でもあり、グラフに表示するにはサーモグラフィのような二次元マッピングが必要です。一般的に横軸を時間、縦軸を周波数とし、STFTの振幅Aを $20log_{10}(A)$ （デシベル単位）に変換して色で表現したグラフを、スペクトログラム（Spectrogram）と呼びます。スペクトログラムは、音声処理の分野ではポピュラーなデータの可視化手法です。

f:id:yuki0718:20190728035947p:plain — 第６回の記事で紹介した音楽のスペクトログラム

調波打楽器音分離（HPSS）

第６回の記事では既存の音楽にコンサートホールのような反響を加えて加工し、第７回の記事では無からギター風の音色を生成しました。そして、音の加工と音の生成ができたので、次は音楽の中から特定の成分を分離するタスクを実装してみようと思い、簡単にできそうな例を探しました。

そして偶然発見したのが、2009年頃に東大工学部の守谷／亀岡研究室というグループが研究されていた「調波打楽器音分離（HPSS：Harmonic/Percussive Sound Separation）」というアルゴリズムです。私はこちらの論文を参考にしました。

最も基本的なHPSSの理論

論文の中で説明されているHPSSのアイディアは、以下のとおりです。

（１）音楽には大きく分けて「コード楽器音」と「歌声」と「打楽器音」の３種類が含まれている。コード楽器音は旋律を奏でるギターやピアノなどの音、打楽器音はビートを刻むドラムやベースなどの音である。

（２）音楽のスペクトログラムの振幅 $A(t, F)$ （時間tと周波数Fの関数）は、上述した３種類の重ね合わせで表現される。そして一般的に、コード楽器音は同じ音が比較的長い時間継続するため時間方向に滑らかな分布をしていることが多く、逆に打楽器音は周波数方向に滑らかな分布をしていることが多い。歌声はどちらかというと時間方向に滑らかだが、はっきりとどちら側とは言えない。

（３）振幅 $A(t, F)$ は、時間方向に滑らかなHarmonic成分 $H(t, F)$ と周波数方向に滑らかなPercussive成分 $P(t, F)$ の足し算で構成されると仮定する。つまり、 $A(t, F)=H(t, F)+P(t, F)$ とする。

（４）コンピュータで複素数を扱うのは大変。そこで、数学的に単純化するため、ちょっと大胆に $H(t, F)$ と $P(t, F)$ の位相は常に揃っていると仮定する。すると、この問題は振幅の絶対値（実数） $|A(t, F)|=|H(t, F)|+|P(t, F)|$ に関する問題として扱える。事前に $A(t, F)$ の位相を何かの変数に記憶しておき、H,Pの絶対値が得られたら複素数に復元すればよい。

（５）「時間／周波数方向に滑らか」とは、スペクトログラムの二次元マップにおいて「横隣り／縦隣りにある色の差が小さい」ことを意味する。これを数学的に表現するには、以下のコスト関数Jを最小にするような $|H_{n, k}|$ と $|P_{n, k}|$ を見つければよい（時間と周波数はデジタル化されているので整数の番号nとkに置き換えた）。

$\displaystyle{J(H, P) \stackrel{\mathrm{def}}{\equiv} \sum_{n=0}^{N-1}\sum_{k=0}^{K-1} {w_{H}(|H_{n+1, k}|^{0.5}-|H_{n, k}|^{0.5} )^{2}+w_{P}(|P_{n, k+1}|^{0.5}-|P_{n, k}|^{0.5} )^{2} } }$

本当は、人間の耳の感度に合わせると振幅絶対値（|H|や|P|）の0.6乗くらいが適切だが、ルートの方が計算が楽なので0.5乗にした（らしい）。

（６） $|A_{n,k}|=|H_{n,k}|+|P_{n,k}|$ という拘束条件付きで、コスト関数 $J(H, P)$ を最小化するようなHとPは、ラグランジュの未定乗数法を用いて解析的に結果を導ける。計算は論文に書いてあるので省略するが、結果は以下のようにHとPを選べばよいことになる。

$\displaystyle{|H_{n, k}|^{0.5}=\frac{w_{H}(|H_{n+1, k}|^{0.5}+|H_{n-1, k}|^{0.5})}{\sqrt{w_{H}^{2}(|H_{n+1, k}|^{0.5}+|H_{n-1, k}|^{0.5})^{2}+w_{P}^{2}(|P_{n, k+1}|^{0.5}+|P_{n, k-1}|^{0.5})^{2} } }|A_{n,k}|^{0.5} }$

$\displaystyle{|P_{n, k}|^{0.5}=\frac{w_{P}(|P_{n, k+1}|^{0.5}+|P_{n, k-1}|^{0.5})}{\sqrt{w_{H}^{2}(|H_{n+1, k}|^{0.5}+|H_{n-1, k}|^{0.5})^{2}+w_{P}^{2}(|P_{n, k+1}|^{0.5}+|P_{n, k-1}|^{0.5})^{2} } }|A_{n,k}|^{0.5} }$

（７）ラグランジュの未定乗数法で解析的に導かれた結果は、ある点（n, k）における値を縦隣りと横隣りの点(n+1, k)、(n-1, k)、(n, k+1)、(n, k-1)における値で表した、漸化式の形をしている。よって、HとPの初期条件（例えば両方Aの1/2にする等）を適当に決めて右辺を計算し、その結果で再び右辺を計算する、というイテレーションをコンピュータで繰り返せば、逐次計算で近似的な解が得られる。

長くなりましたが、やっている計算は大学の初等数学なので、論文を見れば容易に数式は追えると思います（10年ぶりに研究を始めた私程度のレベルでも大丈夫でした）。

PythonでHPSSを実装

例のごとく、いきなりPythonのソースコードを載せます。数式はややこしいですが、コンピュータに処理させる内容は極めてシンプルなので、コード自体はさほど複雑ではないです。

import soundfile as sf
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sg

#■Harmonic-Percussive_Sound_Separation(HPSS)の更新用関数を定義■
def HPSS(A):
    
    #論文の更新式がルートで表記されているので平方根をとる
    rootA = np.sqrt(A)
    
    #HarmonicとPercussiveの和は振幅なので初期条件は均等にしておく
    rootH = np.sqrt(0.5) * rootA
    rootP = np.sqrt(0.5) * rootA
    rootH_new = rootH
    rootP_new = rootP
    
    #行数と列数を取得
    row = A.shape[0]
    column = A.shape[1]
    
    #更新時の重みパラメータとりあえず1にする
    Wh = 1
    Wp = 1
    
    #更新回数は5～50回くらいがいいらしい
    epoch = 5
    
    #更新を単純にepoch回繰り返す
    for x in range(epoch):
        
        #スペクトログラム上を走査して論文の更新式どおりに値を更新
        #周波数軸を走査(行番号kの繰り返し)
        for k in range(1, row-1):
            
            #時間軸を走査(列番号tの繰り返し)
            for t in range(1, column-1):
                #rootHの更新式
                temp1 = Wh * (rootH[k, t+1] + rootH[k, t-1]) * rootA[k, t]
                temp2 = (Wh**2) * (rootH[k, t+1] + rootH[k, t-1])**2 + (Wp**2) * (rootP[k+1, t] + rootP[k-1, t])**2
                rootH_new[k, t] = temp1 / np.sqrt(temp2)
                
                #rootPの更新式
                temp1 = Wp * (rootP[k+1, t] + rootP[k-1, t]) * rootA[k, t]
                temp2 = (Wh**2) * (rootH[k, t+1] + rootH[k, t-1])**2 + (Wp**2) * (rootP[k+1, t] + rootP[k-1, t])**2
                rootP_new[k, t] = temp1 / np.sqrt(temp2)
            
        #値をスペクトログラム全体で一気に更新
        rootH = rootH_new
        rootP = rootP_new
        rootA = np.sqrt(rootH**2 + rootP**2)
    
    return [rootH**2, rootP**2]

#■メイン処理■
if __name__ == "__main__":
    
    #ファイルパス
    fname = './song.wav'
    
    #波形データwavdataとサンプリング周波数Fsを取得
    data, Fs = sf.read(fname)
    
    #ステレオ音源の場合は左右の音を足して2で割ってモノラルに変換
    if(data.shape[1] == 1):
        wavdata = data
    else:
        wavdata = 0.5 * data[:, 1] + 0.5 * data[:, 0]
    
    #曲が長すぎる場合は1分～1分30秒の区間でトリミング
    #wavdata = wavdata[Fs*60 : Fs*90]
    
    #窓関数の種類はハニング窓,窓の大きさは30msec,窓の走査幅は15mecで設定
    window_size = 0.03*Fs #30msec
    window_hop = 0.015*Fs #15msec
    overlap = window_size - window_hop
    
    #短時間フーリエ変換(STFT)を実行
    F, T, dft = sg.stft(wavdata, fs=Fs, window='hann', nperseg=window_size, noverlap=overlap)
    
    #STFTの振幅と位相を取得
    A = np.abs(dft)
    arg = np.angle(dft)
    
    #Original soundのスペクトログラムを表示
    #plt.figure(figsize=(16, 5))
    #plt.title('Spectrogram of the song', fontsize=12)
    #plt.xlabel('Time [sec]')
    #plt.ylabel('Frequency [kHz]')
    #plt.pcolormesh(T, F/1000, 20*np.log10(A), cmap = 'rainbow')
    #plt.colorbar(orientation='horizontal').set_label('Amplitude [dB]')
    #plt.show()
    
    #自己定義関数を呼び出してHarmonicの振幅HとPercussionの振幅Pを得る
    [H, P] = HPSS(A)
    
    #論文ではSTFT時の位相が保存されるという仮定を置いている
    Hdft = H * np.exp(1j*arg)
    Pdft = P * np.exp(1j*arg)
    Adft = A * np.exp(1j*arg)
    
    #Harmonic partのスペクトログラムを表示
    #plt.figure(figsize=(16, 5))
    #plt.title('Spectrogram of the harmonic part', fontsize=12)
    #plt.xlabel('Time [sec]')
    #plt.ylabel('Frequency [kHz]')
    #plt.pcolormesh(T, F/1000, 20*np.log10(H), cmap = 'rainbow')
    #plt.colorbar(orientation='horizontal').set_label('Amplitude [dB]')
    #plt.show()
    
    #Percussion Partのスペクトログラムを表示
    #plt.figure(figsize=(16, 5))
    #plt.title('Spectrogram of the harmonic part', fontsize=12)
    #plt.xlabel('Time [sec]')
    #plt.ylabel('Frequency [kHz]')
    #plt.pcolormesh(T, F/1000, 20*np.log10(P), cmap = 'rainbow')
    #plt.colorbar(orientation='horizontal').set_label('Amplitude [dB]')
    #plt.show()
    
    #逆STFTして分離後の音声データを復元
    _, Adata = sg.istft(Adft, fs=Fs, window='hann', nperseg=window_size, noverlap=overlap)
    _, Hdata = sg.istft(Hdft, fs=Fs, window='hann', nperseg=window_size, noverlap=overlap)
    _, Pdata = sg.istft(Pdft, fs=Fs, window='hann', nperseg=window_size, noverlap=overlap)
    
    #IPythonライブラリを使ってjupyter notebook上に再生バーを表示
    print('Original Sound')
    ipd.display(ipd.Audio(data=Adata, rate=Fs))
    sf.write('original.wav', Adata, Fs)
    
    print('Harmonic Part')
    ipd.display(ipd.Audio(data=Hdata, rate=Fs))
    sf.write('harmonic.wav', Hdata, Fs)
    
    print('Percussive Part')
    ipd.display(ipd.Audio(data=Pdata, rate=Fs))
    sf.write('percussive.wav', Pdata, Fs)

実行結果

今回は音分離のタスクなので、どうしても元となる音楽（望ましくは歌付き）が必要です。とはいえ自分で作詞・作曲するようなスキルもないため、「魔王魂」さんというフリー楽曲を提供しているウェブサイト*5から音楽を拝借しました。

maoudamashii.jokersounds.com

それでは、以下に元の曲、分離されたHarmonic成分（コード楽器音＋歌声）、分離されたPercussive成分（打楽器音）の順に結果を掲載します。

1曲目（テクノポップっぽい曲）

【Harmonic（Epoch=25, Window Width=30msec）】

【Percussive（Epoch=25, Window Width=30msec）】

2曲目（ジャズっぽい曲）

【Harmonic（Epoch=25, Window Width=30msec）】

【Percussive（Epoch=25, Window Width=30msec）】

いかがでしょうか。個人的には、極めてシンプルなアルゴリズムの割に効果的な分離ができている印象を受けました。1曲目はシンセサイザーの旋律がHarmonicできれいに拾えており、2曲目は比較的ドラムのビート音が強いのでPercussiveによく分離できています。

ちなみに、論文をよく読むと後半部分に「STFTのWindow Widthを30msecではなく200msec程度まで長くすると、今度は歌声がPercussive側に分離されて、Harmonicはコード楽器音の旋律だけになる」旨の記載があったので、ちょっと1曲目でWindow Widthを変えてみました。

【Harmonic（Epoch=25, Window Width=200msec）】

【Percussive（Epoch=25, Window Width=200msec）】

確かに論文で言っているとおりの結果になりました。Harmonicが旋律だけになって寂しい感じになりましたが、確かにPercussiveの方に歌声が分離されています。なんかSONYのWalkmanにこんな機能（カラオケモード？）あったなあ、と思いました。

今日のところは以上です。数式多めの記事でTex表記するのは意外と大変だったものの、記事にまとめたらすごく頭の中が整理されました。次回は研究から離れておもしろデンマーク生活的な記事が書けたらと思っています。

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

*1:FFT（Fast Fourier Transformation）などでデジタルデータをフーリエ変換するとき、この仮定が前提として存在することを忘れてはいけません。

*2:疑い深い方なら、「デジタル化されていることによっても相互変換の厳密性を欠くんじゃないか」と思うかもしれません。しかし実は、積分がシグマの無限和に代わり、ディラックのデルタ関数がクロネッカーのデルタに代わるだけで、概ね数学的な厳密性は欠きません。この「概ね」というのはサンプリング定理（フーリエ変換でデジタルデータからアナログ波形を再現できるのは周波数Fs/2まで）と関係しているのですが、説明が面倒なので省略します。

*3:例えば、「あいうえお」と喋った音声データ全体をフーリエ変換しても、「あ」～「お」までの周波数分布の平均が採れるだけで、あまり意味がありません。むしろ、分析者は「あ」～「お」の周波数分布を個別に取り出したいと考えるでしょう。

*4:もちろん、窓幅（Window Width）と窓のずらし幅（Window Shift）を指定すれば、計算はコンピュータが勝手にやってくれます。

*5:アマチュアとは思えないほどクオリティが高いにも関わらず、注意書きを付ければ楽曲の転載を許可してくださっている、たいへん素晴らしいウェブサイトです。