高速フーリエ変換(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のサイン波です。実行すると、
こうなります。前回の離散フーリエ変換の結果と同じです。もう1個、別の波形で試してみます。基本周波数250Hzののこぎり波です。
これも前回と同じ結果が得られました。この程度の例だと高速になったのかほとんどわかりません。リアルタイムで音声を分析したりすると違いがわかるのかもね。
次回は、フーリエ変換の前処理としてよく使われる窓関数の効果を試してみたいと思います。