コンピュータを研究に使うための私的メモ。Python、Fortran、Ubuntu、etc...


FFT

scipy.fftpackを使う。

fft/ifft

複素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を使った場合と違って正変換して逆変換するとちゃんと元に戻る。

rfft/irfft

複素フーリエ変換で返り値が実数配列になる。y = fftpack.fft(data)の要素とz = fftpack,rfftの返り値の対応は以下の通り。
  • nが偶数のとき
zz[0]z[1]z[2]z[3]z[4]...z[n-2]z[n-1]
yy[0].realy[1].realy[1].imagy[2].realy[2].imag...y[n/2 - 1].imagy[n/2].real
  • nが奇数のとき
zz[0]z[1]z[2]z[3]z[4]...z[n-2]z[n-1]
yy[0].realy[1].realy[1].imagy[2].realy[2].imag...y[n/2].realy[n/2].imag
(注) nが奇数のときn/2は小数点以下切り捨てで(n-1)/2になる。

FFTユーティリティ

並べ替え

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

fftで返される複素数配列の各要素はそれぞれ以下の周波数に対応している。
  • nが奇数
インデックス012...n/2-1n/2+1...n-2n-1
周波数01/n2/n...(n-1)/2n-(n-1)/2n...-2/n-1/n
  • nが偶数
インデックス012...n/2-1n/2n/2+1...n-2n-1
周波数01/n2/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
#周波数 spd1日あたりのサンプリング
freq = fftpack,fftfreq(ntim, d = 1./spd)
×

この広告は60日間更新がないwikiに表示されております。

情報

管理人のみ編集できます