ディジタル音声処理特論レポート

課題３「NMFによる音源分離」

学籍番号：G22TX0XX，　氏名：○○ ○○

　6.3節ではマイクロホンアレイによるビームフォーミングを行いました。特定の方向の音を収集する方法と説明しましたが、これは複数の音が混合した音源から目的の音源を分離する方法とも言い換えられます。このようなタスクは音源分離と名付けられ、様々な状況における様々な分離方式が日々提案される一大研究分野となっています。本節では音源分離タスクの概観について触れたのち、ビームフォーミングでは対応できないような状況でも分離が可能な手法の一つ、非負値行列因子分解を実装することで音源分離に対する理解をもう少し深めていきましょう。

6.x.1 音源分離について

まず、音源分離について簡単に解説しておきましょう。音源分離とは先ほど述べた通り、混合音源から各音源信号を分離する技術のことです。分離対象の音源には様々なものがあり、例えばZoom等のオンライン会議システムにおいて会議に参加していない人の音声が入らないように話者ごとの音声を推定する話者分離、通話や録音をした際に入ってしまう雑音を除去する雑音除去、各楽器の音が混合した楽曲から目的の楽器の音だけを取り出す音楽信号分離など、音源分離はあらゆる音に対して適用される技術なのです。2023年現在、このような技術はLINE,Discord,Zoomといった通話機能を搭載したソフトウェアの多くに搭載されています。また音楽の文脈でも活用が進んでおり、代表的なDJソフトの１つseratoや代表的なミキシング/マスタリングVSTプラグインのizotopeシリーズ等に音源分離技術は使われています。

音源分離技術の概観をつかむために、様々な条件によって分類してみましょう[2]。まず、観測信号のチャネル数による違いがあります。

• 単一チャネル信号

単一チャネル信号は要するにマイク1本で録音したモノラル信号で、空間の情報が全く得られないので最も分離が難しい条件になります。しかし、手軽な録音方法でもあり、分離対象がこの形で観測されることもとても多いです。

• 劣決定条件

音源数＞マイク数 の場合のことを指します。よく見られる形としては、2chステレオの音楽信号があります。劣決定条件の場合でも空間に関する情報が得られるため、これをもちいることが多くあります。

• 優決定条件

音源数≦マイク数 の場合を指します。マイクロホンアレイで録音した場合はこれにあたることもあります。得られる空間に関する情報がとても多く、高性能な音源分離が可能です。

また、事前情報や前提条件の面からも各方式を分類できます。

• ブラインド音源分離(Blind Source Separation/BSS)

一切のヒントが与えられないまま音源分離を行います。単一チャネルでは何らかの仮定の元でクラスタリングを行います。劣決定条件であれば時間周波数マスキング、優決定条件であれば独立成分分析(ICA)で分離を行います。

• 空間情報

マイクアレイがあれば空間の情報は与えられます。よって、6.3で行ったビームフォーマはここに分類されます。

• 音色情報

音色に関するサンプルが与えられている場合（目的音だけがなっている時間のある楽曲や、大量に目的音と同じ楽器のサンプルを集めたMUSDB18のようなデータセットが与えられている場合）、教師ありの非負値行列因子分解(Non-Negative Matrix Factrization/NMF)や深層学習(Deep Neural Network/DNN)を利用して音源モデルから分離を行えます。

• その他種々のヒント

楽譜の情報を与えられる場合や、演奏している動画を同時に与えた場合などがこれに該当します。現代ではDNNで音源分離を行う際に追加情報として同時にモデルに入力するという方法で取り入れられる研究が多くあります。

6.3で行ったビームフォーマは優決定条件かつ空間情報が与えられている場合に音源分離が可能ですが、単一チャネルで空間情報が一切与えられていない状況では刃がたちません。状況によって選択すべき方式は全く異なってくるのです。単一チャネルにおいて音源分離を行えるシステムの１つが教師ありNMFです。今回はこれを実装して先ほどとは全く違った音源分離の世界を体感してみましょう。

6.x.2 非負値行列因子分解

非負値行列因子分解は、画像処理の文脈で最初に登場した技術です[4]。当時は顔画像から顔パーツを抽出するのが目的でした。画像のように、実世界で得られる物理的量は非負値であることが多いです。よって、そのようなデータを非負制約付きの任意基底数で低ランク近似するのが非負値行列因子分解の骨子となります。同じように、振幅スペクトログラムに対して適用すれば頻出スペクトルという形で特徴が得られるという発想から、NMFは音響信号処理へ盛んに応用されてきました。

要するにNMFが何をしたいのか、式で示すとこのようになります。

$$ \bf X \approx HU $$

行列XをHとUの積で近似する、これだけです。行列Hは基底行列または辞書行列と呼ばれ、Xに頻出するスペクトルが表れます。つまり、Xが音声のスペクトログラムであれば音色の辞書がHとなります。そして、Hに集められたスペクトルの1本1本がどのタイミングでどのような強さで発現しているか、の情報がUに現れ、これは活性化行列またはアクティベーション行列といいます。

行列XをHとUの積で近似するというタスクを解けるようにもう少し考えると、次のような式で表されることがわかります。

$$ arg \min_{\bf H,U} \mathcal{D}(\bf X|HU \rm ) $$

つまり、何らかの距離関数Dを設定して、これをコストとして変数H,Uについての最小化を行う形になります。距離関数は任意に決められますが、以下の3つが特に有名です。(式を簡単にするため、yのxからの距離の形で示す。)[5]

• 二乗誤差 $$ \mathcal{D}_{EU}(y|x) = (y-x)^2 $$

• Iダイバージェンス $$ \mathcal{D}_{KL}(y|x) = y\log{\frac{y}{x}}-y+x $$

• 板倉斎藤擬距離 $$ \mathcal{D}_{KL}(y|x) = \frac{y}{x}-\log{\frac{y}{x}}-1 $$

あとはこれらを最小化するH,Uを求めればよいのだが、残念ながら非負制約付きの非線形最適化問題であることから、解析的に解くことはできません[5]。そこで、補助関数法というものが考案されました。目的関数の上限となる関数をJensenの不等式から導き出すことが出来るため、こちらを更新して降下させていくことで目的関数を間接的に降下させていくことができるのです。

導出は非常に大変なので専門書や[5]を読んでいただくとして、更新式は結果的に以下のようになります。（Iダイバージェンスの場合）

$$ \bf H \leftarrow H \odot \frac{\frac{XU}{HU}}{U} $$$$ \bf U \leftarrow U \odot \frac{\frac{XH}{HU}}{H} $$

（$ \odot $はアダマール積（行列の要素ごとの積）)

実際にはこの式の通りに更新していくだけでNMFは実装できてしまいます。

NMFを実装して、実際に入力スペクトログラムを行列HとUの積で近似してみましょう。

# NMF
def nmf(S=None,k=None,iter=500,H=None):
    '''NMFを行う関数の定義
      引数 S: 入力スペクトログラム
          k:  基底数
          iter: 繰り返し回数
          H: 基底の初期値
    '''
    k=sum(k)
    l,t = S.shape
    U = np.random.rand(k,t)
    print(S.shape)
    print(U.shape)
    print(H.shape)
    for i in range(iter):
        a = np.tile(np.sum(U,1),(l,1))
        b = H * (np.dot(U,(S/np.dot(H.T,U)).T))

        H = b / a.T + 1e-22

        d = np.tile(np.sum(H,1),(t,1))
        e = U * (np.dot(H,S/np.dot(H.T,U)))

        U = e / d.T + 1e-22
    return H, U

音声を読み込みます。

# 音声の読み込み
fs, wave_data = scipy.io.wavfile.read ('sample/myAIUEO.wav')

音を振幅スペクトログラムの形で扱うことが特徴であるため、短時間フーリエ変換、逆短時間フーリエ変換が必要になります。ここまで音を時間周波数領域で扱うプログラムは自分の手で実装してきましたが、それは本節の本質的な部分ではないので、一回既存のライブラリを使ってみましょう。scipy等の音声を扱えるライブラリには基本的に短時間フーリエ変換を行ってスペクトログラムを簡単に算出できる関数があり、逆短時間フーリエ変換もしっかり実装されているので今回はこれを使ってみましょう。ちなみに、scipy.signalに含まれる関数を使うには、scipyではなくscipy.signalを読み込まなければならないので気を付けましょう。

# ライブラリ読み込み
import scipy.signal

# 短時間フーリエ変換
f, t, csp_data = scipy.signal.stft(wave_data, # 時系列データ
                                   fs = fs, # サンプリング周波数
                                   window = 'hann', # 窓関数
                                   nperseg = 512, # 各セグメントの長さ(FFT点数)
                                   noverlap = 512//2, # オーバーラップ
                                   )
# 振幅スペクトログラム
psp_data = np.abs(csp_data) ** 2

# スペクトログラムの描画
plt.figure()
plt.subplot(2, 1, 1)
plt.title('Waveform')
plot_wave([], wave_data)
plt.subplot(2, 1, 2)
plt.pcolormesh(t, f, 20*np.log10(psp_data), vmin=0, vmax=20*np.log10(psp_data).max()) # 振幅はdBで表して見やすくしてある
plt.title('Magnitude Spectrogram')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [s]');

NMFに掛けてHとUを求めます。

H, U = nmf(S=psp_data, k=10, iter=500, H=np.random.rand(10, 257))

(257, 334)
(10, 334)
(10, 257)

# NMFの可視化
fig, ax = plt.subplots(2, 2,
                       gridspec_kw={
                           'width_ratios': [1, 4],
                           'height_ratios': [1, 4]},
                       figsize=(12,8))
fig.suptitle("Decomposing Amplitude Spectrograms with NMF", fontsize=16)
ax[0][0].axis("off")

ax[1][0].pcolormesh(range(10), f, 20*np.log10(H.T), vmin=0, vmax=20*np.log10(H).max())
ax[1][0].set_title('Basis H')
ax[1][0].set_ylabel('Frequency [Hz]')
ax[1][0].set_xlabel('k');

ax[0][1].pcolormesh(t, range(10), 20*np.log10(U), vmin=0, vmax=20*np.log10(U).max())
ax[0][1].set_title('Activation U')
ax[0][1].set_ylabel('k')

ax[1][1].pcolormesh(t, f, 20*np.log10(psp_data), vmin=0, vmax=20*np.log10(psp_data).max())
ax[1][1].set_xlabel('Time [s]');

実際に基底Hにスペクトログラムに出てきているスペクトルが集まっていますね。そして、活性化行列Uの値が立っているところでそのスペクトルがスペクトログラムに表れているのもこの図からよくわかると思います。実際にHとUの積がスペクトログラムになっており、HとUにはそれぞれ特徴量が表れているのです。

それでは、原理をつかんだところで実際の音源分離を行ってみましょう。

6.x.3 半教師付きNMFによる音源分離

NMFを音源分離につなげるにはいくつか方法があります。

H内の基底をクラスタリングして、目的音のスペクトルだけを集めて対応するUを用いて再構成を行うと、目的の音が取り出せます。この方式は空間情報も教師情報も一切必要ない完全なブラインド音源分離です。

しかし、いくつか挑戦例はあるのですがクラスタリングに相当な精度が必要になってしまうため、高い性能を出すことは難しいのです。

そこで、今回は半教師付きNMF(semi-supervised NMF/SSNMF)を実装してみましょう。 半教師付き、というのはどういう意味かというと、目的の音源の辞書行列(=基底)だけ必要になるという意味になります。手順としては、学習ステージ→分離ステージのように行います。

学習ステージでは、目的の音源のサンプルをNMFに掛けて音色(=スペクトル)の辞書を学習させます。ここで得た辞書をTとします。分離ステージでは、学習ステージで得たTを用いて、他の音源の基底を表現するHを追加してNMFに掛けます。

つまり、分離時に解くNMFでは、学習ステージで得られた辞書をそのまま使うため、以下の式で表されます。

$$ arg \min_{\bf G,H,U} \mathcal{D}(\bf X_{\rm mix}|TG+HU \rm ) $$

ここで、学習した基底T以外を更新しているところがポイントです。これによって、最終的に

$$ \bf Y_{目的} \approx TG $$

で分離された音源を推定することが出来るのです。

このように実装した半教師付きNMFが次のセルのものになります。上の式におけるTを更新しない以外は先ほどのNMFと同じものになります。

# 半教師付きNMF
def ssnmf(S=None,k=None,iter=500,H=None):
    '''SSNMFを行う関数の定義
      引数 S: 入力スペクトログラム
          k:  基底数
          iter: 繰り返し回数
          H: 基底の初期値
    '''
    k=sum(k)
    l,t = S.shape
    U = np.random.rand(k,t)
    print(S.shape)
    print(U.shape)
    print(H.shape)
    for i in range(iter):
        a = np.tile(np.sum(U,1),(l,1))
        b = H * (np.dot(U,(S/np.dot(H.T,U)).T))

        H_ = b / a.T + 1e-22

        H[k//2:, :] = H_[k//2:, :]

        d = np.tile(np.sum(H,1),(t,1))
        e = U * (np.dot(H,S/np.dot(H.T,U)))

        U = e / d.T + 1e-22
    return H, U

それでは、実際に分離を行ってみましょう。

今回のタスク設定は、sample.wavにおける前半の音を既知のものとして、sample.wavの後半＋sample2.wav(ドラム音)の混合音源からsample.wavの後半の音を分離するというものにします。sample.wavの前半の音を教師音源として話者特有の基底を学習し、未知であるsample.wavの後半の音を分離します。

まずはsample.wavを分けましょう。

# ターゲット音源の読み込み
fs, smp_data = scipy.io.wavfile.read('sample/sample.wav')
# 分割
f_data = smp_data[:89000]
s_data = smp_data[89000:]
plt.subplot(2,1,1)
plt.title('f_data')
plot_wave([], f_data)
plt.subplot(2,1,2)
plt.title('s_data')
plot_wave([], s_data)

それでは、混合音源を作成してみましょう。

ドラム音源であるsample2.wavを読み込みます。

# ドラム音源の読み込み
fs, drum_data = scipy.io.wavfile.read('sample/sample2.wav')
drum_data = drum_data[:len(s_data)] # 混合する音源と同じ長さにそろえる
# 可視化
plt.figure()
plt.subplot(2, 1, 1)
plt.title("sample2.wav");
plot_wave([], drum_data);
plt.subplot(2, 1, 2)
plt.specgram(drum_data, NFFT=512, Fs=fs, noverlap=512//2);
plt.ylabel('Frequency [Hz]');
plt.xlabel('Time [s]');

/usr/local/lib/python3.10/dist-packages/matplotlib/axes/_axes.py:7773: RuntimeWarning: divide by zero encountered in log10
  Z = 10. * np.log10(spec)

# ドラム音源を聞いてみよう
audio = Audio(drum_data, rate=fs)
audio

# 混合音源の作成
mixed_data = s_data[:len(drum_data)] + drum_data
# 可視化
plt.figure()
plt.subplot(2, 1, 1)
plt.title("mixed wave");
plot_wave([], mixed_data);
plt.subplot(2, 1, 2)
plt.specgram(mixed_data, NFFT=512, Fs=fs, noverlap=512//2);
plt.ylabel('Frequency [Hz]');
plt.xlabel('Time [s]');

# 混合音源を聞いてみよう
audio = Audio(mixed_data, rate=fs)
audio

# 目的音の振幅スペクトログラムを得る
f, t, f_csp = scipy.signal.stft(f_data, # 時系列データ
                                   fs = fs, # サンプリング周波数
                                   window = 'hann', # 窓関数
                                   nperseg = 512, # 各セグメントの長さ(FFT点数)
                                   noverlap = 512//2, # オーバーラップ
                                   )
f_psp = np.abs(f_csp) ** 2

学習ステージになります。サンプル音源から基底を学習しましょう。

# 基底の学習
k_smp = 10
H_smp, U_smp = nmf(f_psp, k=k_smp, iter=1000, H=np.random.rand(k_smp, len(f_psp)))

(257, 349)
(10, 349)
(10, 257)

それでは、混合音源のスペクトログラムに対して半教師付きNMFを実行して分離を行ってみましょう。

# 混合音源のスペクトログラム
# 目的音の振幅スペクトログラムを得る
f, t, mix_csp = scipy.signal.stft(mixed_data, # 時系列データ
                                   fs = fs, # サンプリング周波数
                                   window = 'hann', # 窓関数
                                   nperseg = 512, # 各セグメントの長さ(FFT点数)
                                   noverlap = 512//2, # オーバーラップ
                                   )
mix_psp = np.abs(mix_csp) ** 2
# ssnmf用基底の作成
init_H = np.concatenate([H_smp, np.random.rand(k_smp, 512//2+1)], axis=0)
# 半教師付きNMFの実行
H_mix, U_mix = ssnmf(mix_psp, k=k_smp*2, iter=1000, H=init_H)

(257, 623)
(20, 623)
(20, 257)

仕上げに目的音源の辞書と活性化行列を掛け合わせて分離音の再構成を行います。

# 分離音の再構成
s_psp = np.dot(H_mix[:10].T, U_mix[:10])
plt.pcolormesh(t, f, 20*np.log10(s_psp), vmin=0, vmax=20*np.log10(s_psp).max())

<matplotlib.collections.QuadMesh at 0x7e10253a2200>

音源分離では、分離前の混合音源の位相を使って分離後の音源を複素スペクトログラムに戻すことが多いです。

# 混合音源の位相から複素スペクトログラムに戻す
mix_asp = np.angle(mix_csp)
s_csp = s_psp * np.exp(1j*mix_asp)
# 逆短時間フーリエ変換に掛ける
sep_t, sep_data = scipy.signal.istft(s_csp, fs=fs, noverlap=512//2, nfft=512)
plt.figure()
plt.title("separated audio")
plot_wave([], sep_data);
plt.figure()
plt.title("oracle audio")
plot_wave([], s_data)

それでは、分離結果を聞いてみましょう！

audio = Audio(sep_data, rate=fs)
audio

そこまで完璧ではないですが、混合音源と比べるとsample2.wavのドラム音がほとんど消えていることが分かるでしょうか。単一チャネルで何も空間情報が得られない状況でも、音色の特徴が分かるだけでここまで分離できるのがNMFの魅力です。現在はこのような問題に対して、データが多くない場合にはNMFで、データがたくさんある場合にはDNNを使うことによってかなりの精度で分離することが出来るようになってきています。興味があれば調べてみてはいかがでしょうか。

参考文献：

[1]"Pythonで学ぶ音源分離",

[2]北村大地, "音源分離における音響モデリング," 日本音響学会 サマーセミナー 招待講演, September 11th, 2017.

[3]

亀岡弘和. (2011). チュートリアル：非負値行列因子分解. 研究報告音楽情報科学（MUS）, 2011(4), 1–1. https://cir.nii.ac.jp/crid/1572543026893302912.bib?lang=ja

[4]Lee, D., Seung, H. Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791 (1999). https://doi.org/10.1038/44565

[5]亀岡 弘和, 非負値行列因子分解の音響信号処理への応用(<小特集>近年の音響信号処理における数理科学の進展), 日本音響学会誌, 2012, 68 巻, 11 号, p. 559-565, 公開日 2017/06/02, Online ISSN 2432-2040, Print ISSN 0369-4232, https://doi.org/10.20697/jasj.68.11_559, https://www.jstage.jst.go.jp/article/jasj/68/11/68_KJ00008328910/_article/-char/ja