scipy.fftpack †

高速フーリエ変換に関するモジュールです．
実際にはfftwなどの下位ライブラリへのラッパーとなっています．
システムの標準ディレクトリにfftwがあればそれを用いるようです（未確認）．

↑

はじめに †

以下の例を実行するには

from scipy import *
from scipy.fftpack import *

しておいて下さい。

↑

1次元FFT †

↑

複素数のFFT †

fftpack.fftを用います．

>>> x = sin(2*pi*arange(64.0)/64.0)
>>> y = fft(x)

これでxをフーリエ変換した結果がyに格納されます。 ここではxは実数の配列ですが、自動的にこれは複素数の実部と解釈されます．

FFTを使うときはデータの並びがライブラリによって異なって混乱しますが、 ドキュメントによればfftpackでは以下のような格納方法が採用されているようです。

   fftで返される配列は
     [y(0),y(1),..,y(n/2-1),y(-n/2),...,y(-1)]        (nが偶数のとき)
     [y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  (nが奇数のとき)
   ただし、
     y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n)
     j = 0..n-1
   である．ここで、y(-j) = y(n-j)であることに注意.
↑

実数のFFT †

fftpack.rfftを用います．

>>> y = rfft(x)

実数の場合は以下のようなデータ並びになっています．

   rfftで返される実数配列は
     [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))]              (nが偶数のとき)
     [y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))]   (nが奇数のとき)
   ただし、
     y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n)
     j = 0..n-1
   である．ここで、 y(-j) = y(n-j)であることに注意．
↑

周波数配列作成 †

fftfreqという周波数配列をfftが返す配列と同じ並びで作成してくれる関数があります．

>>> n, dt = 8, 1.0
>>> fftfreq(n, dt)
array([ 0.   ,  0.125,  0.25 ,  0.375, -0.5  , -0.375, -0.25 , -0.125])

となります．nはデータ数，dtはサンプリング周期です．作成された配列がfftが返す 配列と同じ並びになっているのが分かると思います．
dtはオプションでデフォルトでは1と解釈されます．

また，rfftfreqというrfft用の関数もあります．

>>> rfftfrq(n, dt)
array([ 0.   ,  0.125,  0.125,  0.25 ,  0.25 ,  0.375,  0.375,  0.5  ])

となります．

↑

並べ替え †

fftで返される配列をfftshiftという関数を用いて扱いやすい形に並べ替えることがで きます．
実際にどのように並べ替えられるかを以下で見てみみます．

>>> n, dt = 8, 1.0
>>> freq = fftfreq(n, dt)
>>> fftshift(freq)
array([-0.5  , -0.375, -0.25 , -0.125,  0.   ,  0.125,  0.25 ,  0.375])

つまり周波数が小さい方から大きくなるように並び替えてくれます．

↑

規格化 †

fftで返される配列はそれぞれのモードの実際の振幅をn/2倍したものになっています(n はデータ数)．従って実際の振幅を求めるにはfftの結果を2/n倍します．

>>> y = fft(x) * 2.0 / x.size

また逆変換はデフォルトで，

ifft(fft(x)) == x

が数値誤差の範囲内で成り立つようになっています．従って上の規格化をした場合は

>>> x = ifft(x)* 0.5

が規格化された逆変換となります．

↑

パワースペクトル †

複素数のfftが返す配列を仮定すると

>>> y = fft(x) * 2.0 / x.size # 規格化したFFT
>>> abs(y[0:n/2+1])**2

がパワースペクトルを表す。ここでnはデータ数。 ここではFFTを規格化しているので、パワースペクトルはそれぞれのモードの 振幅の2乗となる。

↑

多次元FFT †

ちょっとデータ並びとかがややこしい．