numpyでFFTを使う

1次元FFT

numpy.fft.fftを使う。
※FFTの結果の格納の順番に注意
最初に周波数プラスのものを昇順に、次に周波数マイナスのものを昇順に、という順番で格納されている。なのでそのままプロットしても結果を把握しづらい。

格納順への対応方法

1. numpy.fft.fftfreqで上記の格納順に対応する周波数の配列を作成
2. numpy.fft.fftshiftでFFTの結果を周波数マイナスのもの〜周波数プラスのものの順に並び替える

以下numpy.fft.fftの使用例

1次元FFT使用例

import numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt


y = []
N = 10000
step = 0.001
# 周波数の配列を作成
freq = fft.fftfreq(N, d=step)
# 周波数、振幅その1
f1 = 2
amp1 =1
# 周波数、振幅その2
f2 = 5
amp2 = 3

for i in range(N):
    x = i*step
    # 周波数、振幅その1、その2の波の重ね合わせ
    y_i = amp1*np.sin(f1*2*np.pi*x) + amp2*np.sin(f2*2*np.pi*x)
    y.append(y_i)

fy = fft.fft(y)
# 絶対値を規格化してプロット
plt.scatter(freq, np.abs(fy)/np.linalg.norm(fy))
plt.xlim([0, 7])
plt.grid()
plt.show()

• 結果
  それぞれの周波数のピークが表れている。絶対値の比もそれぞれの振幅の比と一致する。

2次元FFT

numpy.fft.fft2を使う。
2次元の場合、x、y方向両方とも上記のように周波数プラスのもの〜周波数マイナスのものの順で格納されている。
numpy.fft.fftshiftを使用すればx、y方向両方とも周波数マイナス〜プラスの順に並べ替えてくれる。

2次元FFT使用例

mport numpy as np
import numpy.fft as fft
import matplotlib.pyplot as plt


z_ij = []
N = 1000
mid = N/2
step = 0.001
freq = fft.fftfreq(N, d=step)
freq = freq.tolist()
freq.sort()
freq_x = []
freq_y = []
for i in range(N):
    freq_x.append(freq)
for i in range(N):
    freq_y.append([freq[i]]*N)

# x方向、y方向周波数
fx1 = 4
fy1 = 6

for i in range(N):
    z_j = []
    for j in range(N):
        x = j*step
        y = i*step
        z = 2*np.sin(fx1*2*np.pi*x + fy1*2*np.pi*y)
        z_j.append(z)
    z_ij.append(z_j)

fz_ij = fft.fft2(z_ij)
# fftshiftでFFTの結果を並び替え
shift_fz_ij = fft.fftshift(fz_ij)

plt.contourf(freq_x, freq_y, np.log(np.abs(shift_fz_ij)), cmap='coolwarm')
delta = 10
plt.xlim([0, delta])
plt.ylim([0, delta])
plt.grid()
plt.show()

• 結果
  x方向周波数4、y方向周波数6の場所にピークがある