高速フーリエ変換（FFT） - 人工知能に関する断創録

Pythonで音声信号処理（2011/05/14）

今回は、高速フーリエ変換（FFT）を試してみます。FFTとは~~Final Fantasy Tactics~~ Fast Fourier Transformの略でその名の通り、前回の離散フーリエ変換（DFT）を大幅に高速化したしたアルゴリズムです。一般にフーリエ変換といったらFFTが使われるようです。DFTは自分で公式に忠実に実装してみましたが、FFTはPythonのnumpyやscipyに実装があるのでそれを使ってみます。numpyの実装はnumpy.fft.fftでscipyの実装はscipy.fftpack.fftです。使い方はほとんど同じですが、この記事によるとscipyの実装の方が高速とのこと。scipy版には他にもいろいろ関数があります。おいおい使っていきたいと思います。

#coding:utf-8
import wave
import numpy as np
import scipy.fftpack
from pylab import *

if __name__ == "__main__" :
    wf = wave.open("data/sine.wav" , "r" )
    fs = wf.getframerate()  # サンプリング周波数
    x = wf.readframes(wf.getnframes())
    x = frombuffer(x, dtype= "int16") / 32768.0  # -1 - +1に正規化
    wf.close()

    start = 0  # サンプリングする開始位置
    N = 256    # FFTのサンプル数

    X = np.fft.fft(x[start:start+N])  # FFT
#    X = scipy.fftpack.fft(x[start:start+N])         # scipy版

    freqList = np.fft.fftfreq(N, d=1.0/fs)  # 周波数軸の値を計算
#    freqList = scipy.fftpack.fftfreq(N, d=1.0/ fs)  # scipy版

    amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X]  # 振幅スペクトル
    phaseSpectrum = [np.arctan2(int(c.imag), int(c.real)) for c in X]    # 位相スペクトル

    # 波形を描画
    subplot(311)  # 3行1列のグラフの1番目の位置にプロット
    plot(range(start, start+N), x[start:start+N])
    axis([start, start+N, -1.0, 1.0])
    xlabel("time [sample]")
    ylabel("amplitude")

    # 振幅スペクトルを描画
    subplot(312)
    plot(freqList, amplitudeSpectrum, marker= 'o', linestyle='-')
    axis([0, fs/2, 0, 50])
    xlabel("frequency [Hz]")
    ylabel("amplitude spectrum")

    # 位相スペクトルを描画
    subplot(313)
    plot(freqList, phaseSpectrum, marker= 'o', linestyle='-')
    axis([0, fs/2, -np.pi, np.pi])
    xlabel("frequency [Hz]")
    ylabel("phase spectrum")

    show()

前回と同じくFFTの結果から振幅スペクトルと位相スペクトルを求めています。横軸の周波数は、fftfreq()を用いると簡単に計算できます（前回みたいに自分で計算しても簡単ですけど）。第2引数にはサンプリング周波数の逆数、サンプリング周期を与えます。sine.wavは前回使った基本周波数250Hzのサイン波です。実行すると、

f:id:aidiary:20110618121656p:plain