PythonにおけるGriffinlim波形復元のシンプルな実装

はじめに

今回は、Pythonを用いてスペクトログラムから波形情報を復元する手法について簡単に説明していきます。
注意点として、この記事ではGriffinlimアルゴリズムの仕組みに関する説明は述べません。Python初心者の方向けに、librosaのライブラリ関数の一つであるgriffinlim関数をどのように使えばいいのかを解説する記事になります。

ところで、librosaの関数の仕様を知りたいのですから、当然librosaの公式リファレンスを見るのが最善の手法です。
このリファレンスを読んで理解できる方はこれ以降読む必要はありませんし本来はそうすべきでしょう。
ですが、特に日本人の初学者の場合「慣れない内容を英語で説明されてもよくわからない」というのが本音だと思います。そのような方のためにここに書かれている内容を重要なポイントだけかみ砕いて説明するのが本記事の趣旨になります。

書かれている内容が不正確である可能性がありますが、予めご了承ください。

前提

前提として、rfftを用いた短時間フーリエ変換によってスペクトログラムを得る手法は理解しているものとします。
ここを理解できていない方は、短時間フーリエ変換について説明したほかの記事を参照してください。

目標

スペクトログラムを表すnumpy配列をもとにして、音声波形を復元する。

サンプル実行環境

（多少環境が違っても動くとは思います）
python 3.6.8
librosa==0.9.2
numpy==1.19.5
matplotlib==3.3.4

事前準備

ひとまず、短時間フーリエ変換によってスペクトログラムを得ます。

コード１

import numpy as np
import librosa

SR = 16000

x, _ = librosa.load('input.wav', sr=SR)
size_frame = 2048 #短時間フーリエ変換のフレームサイズ
size_shift = 160 #短時間フーリエ変換のフレームシフト量
hamming_window = np.hamming(size_frame)
spectrogram = []
for i in np.arange(0, len(x)-size_frame, size_shift):
    x_frame=x[int(i):int(i)+size_frame] #短時間フレームの切り出し
    fft_spec = np.fft.rfft(x_frame * hamming_window) #ハミング窓の適用
    fft_log_abs_spec = np.abs(fft_spec) #絶対値に直す
    spectrogram.append(fft_log_abs_spec)
spectrogram=np.array(spectrogram)


print(spectrogram.shape)

このとき、最後の行の出力は、（時間方向のサイズ,周波数方向のサイズ）となっているはずです。
具体的にはおおよそ以下の式で求められます。（小数点の切り上げ等の関係で若干ずれます）

時間方向のサイズ= \frac{len(x)-\mathrm{size\_frame}}{\mathrm{size\_shift}}

\\

周波数方向のサイズ= \frac{\mathrm{size\_frame}}{2}+1

周波数方向のサイズがsize_frameの半分になっているのは、np.rfftが入力配列の半分+1のサイズのnumpy配列を返すという仕様になっているためです。（np.fftの場合は入力配列と同じサイズのnumpy配列を返します）

griffinlim関数の仕様

重要な引数

引数１つ目　スペクトログラムS

librosaの公式リファレンスを見てみると、S:np.ndarray [shape=(…, n_fft // 2 + 1, t), non-negative]　とあります。

これは、shape=(…, n_fft // 2 + 1, t)でかつすべての要素が非負（non-negative）であるようなnumpy配列を入力として与えてください、という意味です。librosaのstft（短時間フーリエ変換）関数の出力のshapeもこれと同じになっています。
ですから、短時間フーリエ変換にstftを使ったならば、その結果をそのままgriffinlimへの入力としてもよい、ということです。

ここで、n_fftはフーリエ変換の式におけるNの値、すなわちnp.fftの出力サイズと同じ意味を表しています。ということはすなわち、n_fft // 2 + 1とはnp.rfftの出力サイズ、さらにはコード１で示したspectrogramにおける「周波数方向のサイズ」と一致していることが分かります。ちなみにtとは時間方向のサイズのことです。

従ってこのshape=(…, n_fft // 2 + 1, t)は、コード１で出力した（時間方向のサイズ,周波数方向のサイズ）と逆になっており、spectrogramをlibrosa.griffinlim関数の引数に与えたければ 転置をする必要がある 、ということが分かります。

引数２つ目　反復回数n_iter

Griffinlimのアルゴリズムは、はじめランダムな位相を仮定し、これを反復計算によって徐々に更新していくことによって正しい波形を復元するアルゴリズムです。ここで重要なのはこれは「反復計算」であるということです。
griffinlim関数では、この計算の反復回数を指定することができます。何も指定しなければ反復回数は32回です。当然、反復回数が少ないほど計算が速く終わりますが、あまり少なくしすぎると精度が低下するため、適当な値を見つける必要があります。特にこだわりがなければ指定しなくても問題ありません。

引数３つ目　フレームシフトhop_length

hop_lengthには、短時間フーリエ変換時のフレームシフト（切り出すフレームをどれだけの幅ずつずらしていくか）を指定します。
指定しなければn_fft//4がデフォルトで指定されますが、短時間フーリエ変換時のsize_frameとsize_shiftの関係が4:1になっていない場合にはエラーが出たり結果がおかしくなったりするため、指定しておくのが無難です。（基本的にはsize_frameがsize_shiftの4倍よりも大きいときにエラーとなります）

引数４つ目　フレームサイズwin_length

win_lengthには、短時間フーリエ変換時のフレームサイズを指定します。指定しなければn_fftがデフォルトで指定されます。前述の通り、np.rfftを用いてスペクトログラムを得ている場合、フレームサイズとn_fftは等しくなっているため何も指定しなくても問題ありません。

その他の引数（これらの引数については無視しても特に問題ありません）

引数５つ目　ウィンドウ関数window

windowには、使用するwindow関数を指定することができます。例えば"hamming"と指定すればハミング窓を指定できます。
正直なところ指定しなくても十分問題なく動作するので、初学者が気にすることではないように思います。
指定できる内容を詳しく知りたい方は公式リファレンスを参照してください。

引数６つ目　中央揃えフラグcenter

短時間フーリエ変換時のフレームが、中央揃え(True)であるか、左揃え(False)であるかをbool値で指定することができます。
デフォルトはTrueです。

引数７つ目　型dtype

信号の型を指定することができます。基本的には入力信号に自動で合わせられるのでそれほど気にする必要はありません。

引数８つ目　長さlength

出力の長さを指定することができます。入力よりも短い長さを指定すれば、指定された長さだけ音声が切り出されて出力されます。逆に長い時間を指定すれば、無音区間によって埋められます。必要がなければ特に指定する必要はありません。

引数９つ目　パディングモードpad_mode

中央揃えフラグがTrueの場合、信号の開始と終了の両端部分をどのように処理するのかを指定できます。必要がなければ特に気にする必要はありません。

引数１０個目　更新速度momentum

Griffinlimアルゴリズムの反復計算における更新速度を指定できます。小さい値を指定するとなかなか収束しなくなり、逆に1を超える値を指定すると発散して処理がうまくいかなくなる可能性があります。分からなければ指定しなくて問題ありません。

引数１１個目　初期推定値init

位相の初期推定値を指定します。デフォルトは'random'でランダムに初期化されます。与えているスペクトログラムSが、すでにある程度位相推定が進んだものであり、位相推定を途中から再開したい場合には'None'を指定することができます。分からなければ指定しなくて問題ありません。

引数１２個目　乱数設定random_state

位相の初期推定値のランダムな初期化において、どのような乱数を用いるかを指定します。デフォルトは'None'ですが、整数値でシードを指定したり、numpyのRandomStateを指定したりすることができます。特に指定する必要はありませんが、処理の結果が初期値依存性によって毎回変化してしまっては困るような場合には、適宜シードを指定することで対処することが可能です。

返り値

返り値１つ目　波形y

復元された波形を表すnumpy配列です。

サンプルコード

コード２

#Griffinlimアルゴリズム

y = librosa.griffinlim(spectrogram.T,hop_length=size_shift,n_iter=50)

#グラフ表示

fig = plt.figure()
ax11 = fig.add_subplot(211)
ax12 = fig.add_subplot(212)
ax11.plot(x)
ax12.plot(y)
plt.show()

結果

まとめ

今回はlibrosaのgriffinlim関数の使い方についてまとめました。使い方が意外と簡単であることが分かっていただけたかと思います。
リファレンスの解釈（特に引数５つ目以降）については筆者が使いこなせていない部分もあるので、間違いなどに気付かれた方はコメント等いただけると助かります。