読者です 読者をやめる 読者になる 読者になる

音楽プログラミングの超入門(仮)

Python / 音楽情報処理 初心者が、初心者にも分かるような記事を書きたい。

Pythonで短時間フーリエ変換(STFT)と逆変換

※ iSTFT の実装が間違っていました。というよりも、スペクトログラムをイジらない場合しか滑らかな再合成ができない不適切な実装となっていたので、まともな(ロバストな)実装に変えました。(今のところ)一番需要がありそうな記事で残念過ぎるミスを犯してすいません・・・。
(2014, 2/7)

周波数(スペクトル)の時間変化を見る

少し前に音響信号を周波数スペクトルに変換するフーリエ変換を扱いました。

Pythonでフーリエ変換(と逆変換) - 音楽プログラミングの超入門(仮)

 ここで、例えば市販曲などの3分程度の音響信号の周波数を分析することを考えます。この音響信号全体をフーリエ変換して得られたスペクトルはあまり意味がありません。このような長い音響信号では、ある時間ではドの音が鳴っていたり、別の時間ではラの音がなっていたりしますが、それらが全部ごちゃまぜで解析されてしまうからですね。各時刻のスペクトルを計算して、その時間変化を見たいですね。
 これを行う処理を短時間フーリエ変換(STFT)と言います。アルゴリズムは簡単で、信号を前から順番に 1.決められた長さを切り出して、その部分をフーリエ変換2.少しシフトして同じ処理を行う、ということを繰り返していきます(下図)。
f:id:yukara_13:20131117183055p:plain
毎フレームに窓関数をかけるのを忘れないようにしましょう。ちなみにスペクトル時間変化を表示した図をスペクトログラムと言います。

周波数分解能と時間分解能

 アルゴリズムは分かりましたが、毎フレームどのくらい切り出せばいいのでしょうか?切り出す幅を大きくすればするほど、多くの周波数を含むので周波数分解能は高くなりますが、時間的に広い範囲をとることになるので細かい時間変化などをとらえにくくなり時間分解能が低くなります。切り出す幅を小さくすると、この逆のことが起こりますね。このように周波数分解能と時間分解能にはトレードオフの関係があり、不確定性原理と呼ばれています。
 実際に使う幅は、サンプリング周波数と目的に応じて適宜決めることになりますが、44.1[kHz]のときは 1024 サンプル、16[kHz] のときは 512 サンプルくらいが一般的だと思います。シフトする幅は(個人的に)切り出す幅の半分とすることが多いです。

実装しよう(逆変換も)

それでは短時間フーリエ変換を実装してみましょう。ここで、フーリエ変換は可逆変換でしたので、短時間フーリエ変換の処理を逆に行なっていけば、逆変換することができます。

# -*- coding: utf-8 -*-
# ==================================
#
#    Short Time Fourier Trasform
#
# ==================================
from scipy import ceil, complex64, float64, hamming, zeros
from scipy.fftpack import fft# , ifft
from scipy import ifft # こっちじゃないとエラー出るときあった気がする
from scipy.io.wavfile import read

from matplotlib import pylab as pl

# ======
#  STFT
# ======
"""
x : 入力信号(モノラル)
win : 窓関数
step : シフト幅
"""
def stft(x, win, step):
    l = len(x) # 入力信号の長さ
    N = len(win) # 窓幅、つまり切り出す幅
    M = int(ceil(float(l - N + step) / step)) # スペクトログラムの時間フレーム数
    
    new_x = zeros(N + ((M - 1) * step), dtype = float64)
    new_x[: l] = x # 信号をいい感じの長さにする
    
    X = zeros([M, N], dtype = complex64) # スペクトログラムの初期化(複素数型)
    for m in xrange(M):
        start = step * m
        X[m, :] = fft(new_x[start : start + N] * win)
    return X

# =======
#  iSTFT
# =======
def istft(X, win, step):
    M, N = X.shape
    assert (len(win) == N), "FFT length and window length are different."

    l = (M - 1) * step + N
    x = zeros(l, dtype = float64)
    wsum = zeros(l, dtype = float64)
    for m in xrange(M):
        start = step * m
        ### 滑らかな接続
        x[start : start + N] = x[start : start + N] + ifft(X[m, :]).real * win
        wsum[start : start + N] += win ** 2 
    pos = (wsum != 0)
    x_pre = x.copy()
    ### 窓分のスケール合わせ
    x[pos] /= wsum[pos]
    return x


if __name__ == "__main__":
    wavfile = "./aiueo.wav"
    fs, data = read(wavfile)

    fftLen = 512 # とりあえず
    win = hamming(fftLen) # ハミング窓
    step = fftLen / 4

    ### STFT
    spectrogram = stft(data, win, step)

    ### iSTFT
    resyn_data = istft(spectrogram, win, step)

    ### Plot
    fig = pl.figure()
    fig.add_subplot(311)
    pl.plot(data)
    pl.xlim([0, len(data)])
    pl.title("Input signal", fontsize = 20)
    fig.add_subplot(312)
    pl.imshow(abs(spectrogram[:, : fftLen / 2 + 1].T), aspect = "auto", origin = "lower")
    pl.title("Spectrogram", fontsize = 20)
    fig.add_subplot(313)
    pl.plot(resyn_data)
    pl.xlim([0, len(resyn_data)])
    pl.title("Resynthesized signal", fontsize = 20)
    pl.show()

こんな感じです。あんまり難しくないですね。STFTをどこから始めるかとか、個人によって細かい流儀がある気がしますが、自分は(確か)何回 STFT と iSTFT を繰り返しても長さが変わらないようにしていると思います。(多分)
 これを実行した結果が下図です。ちなみに、「あいうえお」と言っている音声ファイルを入力しました。それっぽいスペクトログラムが出ていますね。逆変換で元の信号が得られているのが分かると思います。
f:id:yukara_13:20131117205253p:plain

まとめ

短時間フーリエ変換ができるようになれば、サウンドプログラミング中級者であるといっても過言ではありません。(もしかしたら)