【python】FFTで周波数解析【サウンドプログラミング】

こんにちは，1000chです．前回はギター，リコーダーの音を離散フーリエ変換によって周波数軸でのデータに変換しました．

しかし実行してみるとわかるのですが，前回のようなDFTの定義そのままに書いたプログラムでは，処理に時間がかかってしまいます．とてもじゃないですがリアルタイムのデータ解析なんてこたぁできません．

じゃー計算方法を工夫してなんとか素早く処理を行おう，という話になります．そこで考案されたのが，高速フーリエ変換（Fast Fourier Transform，FFT）という手法です．ということで，今回はFFTのアルゴリズムとpythonでの使い方についてみていきます．（ぶっちゃけただfftで周波数解析するだけであれば，アルゴリズムはブラックボックスでもいいんですけどねー）

FFTのアルゴリズム

信号処理の話を調べたとき，やたらと出てくるキーワードの一つが"FFT"だと思います．勘違いしやすいのですが，FFTはフーリエ変換の特殊例という訳ではなく，「離散フーリエ変換を高速に解くアルゴリズム」であることは頭に置いておいた方が良いでしょう．なので，言ってしまえばFFTにより得られる出力は，DFTにより得られる出力と何ら変わりありません．あくまで同じ計算を，いかに効率的に解くか，という話がFFTです．具体的に言うと，DFTがO(N²)，FFTがO(NlogN)だそうです．

では，FFTのアルゴリズムをみていきましょう．参考はこちら．

調べた感じ，FFTには大きく２種類のアルゴリズムがあるそうです． Cooley-Tukey型とPrime-Factor型だそうで．ここではCooley-Tukey型の基本的な考え方について触れます．

まず，DFTの式をおさらいします．

${ \displaystyle X(k) = \sum_{n=0}^{N-1} x(n) {w_n}^k \hspace{3.5em} (0 \leq k \leq N-1) }$

ただし，

${ \displaystyle w_n = \exp( \frac{-j2\pi}{N}) }$

こんな感じでした．

N を偶数と仮定した上で，X(k)について，kの偶奇で場合分けしてみます．

（はてなのLatex書き方に慣れない慣れないので，おいおい数式化します）

k = 2l(l = 0~N/2)の場合

f:id:ism1000ch:20140527153848j:plain

k=2l+1の場合 (l = 0 ~ (N-1)/2)

f:id:ism1000ch:20140527155117j:plain

さて，これらをまとめると，X(k)を求めることができます．

f:id:ism1000ch:20140527155836j:plain

計算量を考えましょう．

DFTはO(N²)なので，N/2サイズのDFTは計算量が1/4になります． k = 2l, k=2l+1の2回分計算することを考えると，全体の計算量は1/2になる計算です．

さて，今DFTの計算をkの偶奇で２つのDFTに分解することで，計算量が1/2になることがわかりました．分割されたDFTに対し，この操作を繰り返すと，計算量は1/4,1/8,..とどんどん減っていくことになります．分解するための条件はN = 偶数であることなので，N = 2 ^ n と表現できるとき，DFTはn回分分割できることになります．分割された結果としては，N個のスカラ積となります．そのオーダーを計算すると，O(N * log(N))となるわけですね．

pythonでFFT

とまぁFFTのアルゴリズムがわかったところで，実際にfftを使ってみましょう．

numpyのfftモジュールを使うととても簡単です．

import numpy as np
freq_data = np.fft.fft(time_data) #time_data は時間軸上のデータ，サイズは2 ** n

これだけです． time_dataのサイズを2のn乗にしなくても動作しているのですが...どういうことなんでしょうね...

まあいいか．

テストコードとしてこんなんかきました．

#coding:utf-8

import numpy as np
from wave_loader import *
import matplotlib.pyplot as plt

fs = 8000.0
d = 1.0 / fs
size = 256

wave1 = wave_load("guitar_A4.wav")
wave2 = wave_load("recorder_A4.wav")
dt1 = np.fft.fft(wave1[10000:10000 + size])
dt2 = np.fft.fft(wave2[10000:10000 + size])
frq = np.fft.fftfreq(size,d)
print(len(frq))

plt.subplot(2,1,1)
plt.title("FFT - Guitar A4")
plt.plot(frq,abs(dt1))
plt.axis([0,fs/2,0,max(abs(dt1))])
plt.subplot(2,1,2)
plt.title("FFT - recorder_A4")
plt.plot(frq,abs(dt2))
plt.axis([0,fs/2,0,max(abs(dt2))])

plt.show()

実行結果

f:id:ism1000ch:20140527163055p:plain

前回とほぼ一緒ですが，きちんと分析されていることがわかりますね．

まとめ

FFTのアルゴリズムを勉強
numpyのfftモジュール使ってみた

すこしふしぎ．

VR/HI系院生による技術ブログ．まったりいきましょ．（友人ズとブログリレー中．さぼったら焼肉おごらなきゃいけない）

【python】FFTで周波数解析【サウンドプログラミング】

FFTのアルゴリズム

pythonでFFT

まとめ