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

すこしふしぎ.

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

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

python SoundProcessing

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

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

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

FFTアルゴリズム

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

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

調べた感じ,FFTには大きく2種類のアルゴリズムがあるそうです. 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(N2)なので,N/2サイズのDFTは計算量が1/4になります. k = 2l, k=2l+1の2回分計算することを考えると,全体の計算量は1/2になる計算です.

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

pythonFFT

とまぁ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

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

まとめ