scipy.fftpackを使う。

複素FFT正変換/逆変換を行い複素数配列を返す。

import scipy.fftpack as fftpack
x = np.array(x)
Nfq = len(x)/2         #ナイキスト周波数
f = fftpack.fft(x)
power = f[:Nfq]       #ナイキスト周波数以降は折り返しのため有効なものはここまで

• 返り値は複素数
• f[0]は周波数ゼロ成分。f[1:n/2+1]は正の波数、f[n/2+1:]は負の波数。
• nが偶数の時はf[n/2]は正負の波の合計値が入る。ナイキスト周波数以上の係数は決まらないため。
• 正変換では項数Nで割るような正規化はなされない。周波数ゼロは元データの平均のN倍になっている。

NumPyのリファレンスガイドによればNumPyのFFTの正変換、逆変換は以下で定義される。

• 正変換

• 逆変換

SciPyのリファレンスガイドには記述がないがNumPyと同じ。SciPyの関数リファレンスの逆変換の説明文に1/Nのファクターが抜けているがこれは間違いだと思われる。fftwを使った場合と違って正変換して逆変換するとちゃんと元に戻る。

複素フーリエ変換で返り値が実数配列になる。y = fftpack.fft(data)の要素とz = fftpack,rfftの返り値の対応は以下の通り。

• nが偶数のとき

z	z[0]	z[1]	z[2]	z[3]	z[4]	...	z[n-2]	z[n-1]
y	y[0].real	y[1].real	y[1].imag	y[2].real	y[2].imag	...	y[n/2 - 1].imag	y[n/2].real

• nが奇数のとき

z	z[0]	z[1]	z[2]	z[3]	z[4]	...	z[n-2]	z[n-1]
y	y[0].real	y[1].real	y[1].imag	y[2].real	y[2].imag	...	y[n/2].real	y[n/2].imag

(注) nが奇数のときn/2は小数点以下切り捨てで(n-1)/2になる。

画像変換や波数-周波数解析などで周波数ゼロ成分を中心にするように並び替える必要が出てくる。そのための便利なメソッドが用意されている。

fftで返される複素数配列の各要素はそれぞれ以下の周波数に対応している。

• nが奇数

インデックス	0	1	2	...	n/2-1	n/2+1	...	n-2	n-1
周波数	0	1/n	2/n	...	(n-1)/2n	-(n-1)/2n	...	-2/n	-1/n

• nが偶数

インデックス	0	1	2	...	n/2-1	n/2	n/2+1	...	n-2	n-1
周波数	0	1/n	2/n	...	(n-1)/2n	合計	-(n-1)/2n	...	-2/n	-1/n

これを周波数の小さい側から大きい側にゼロ成分が中心となるように並べ替える。

f = fftpack.fft(x)
f = fftpack.fftshift(f)
#多次元配列の場合は入れ替える軸を指定
f = fftpack.fftshift(f, axes=(0,))

#逆にもとに戻すときは
f = fftpack.ifftshift(f)
x = fftpack.ifft(f)

• サンプルが偶数の場合はf[-1]は係数ではないので注意。ナイキスト以上の成分とみなせばフィルターかけるときは問題ない。
• 小さい方から並べるには

#複素FFTの場合
fftpak.fftshift(fftpack.fft(data))
#実FFTの場合はマイナスはない
fftpack.rfft(data)

FFTの係数に対応する周波数を求めるメソッドもある。

#d=でサンプリング間隔を指定。例えば1sで10点ならd = 0.1 /s
frequency = fftpack.fftfreq(x, d=0.1)

• 例

ntime, nlat, nlon = data.shape
#東西波数
wavenum = fftpack.fftsfreq(nlon)*nlon
#周波数 spd１日あたりのサンプリング
freq = fftpack,fftfreq(ntim, d = 1./spd)