高速フーリエ変換(FFT)- 実装編
FFTのPython3
での実装を目指します。
参考
回転因子(配列)
入力の配列サイズ に対して、 サイズの回転因子を用意すれば良い。
サイズ の回転因子の配列を とすると、
であるので、
1 2 3
| import numpy as np
W_N = np.exp(-2j * np.pi * np.arange(N // 2) / N)
|
FFT
関数は下のようなフローとなる。
FFT(入力配列x
)
- 入力配列を偶数行、奇数行に分割→
x_even
、x_odd
- 偶数行、奇数行のバタフライ演算を再帰的に計算→
F_even
、F_odd
- 出力配列を用意→
F
- 回転因子の配列を用意→
W
- 出力配列の前半分を
F_even + W*F_odd
とする - 出力配列の後半分を
F_even - W*F_odd
とする
↓実際のコード(参考サイト)
FFT1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
| import numpy as np
def FFT(X): N = len(X) n = N // 2
if N == 1: return X F_even = FFT(X[0::2]) F_odd = FFT(X[1::2])
W = np.exp(-2j * np.pi * np.arange(n) / N)
F = np.zeros(N, dtype='complex') F[:n] = F_even + W * F_odd F[n:] = F_even - W * F_odd
return F
|
IFFT
逆高速フーリエ変換(高速フーリエ逆変換?)
↑これを実装すればいいので,コード自体はFFTとそんなに変化なし.
回転因子の の符号を逆にして,最後に をかけるだけです.
IFFT1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
| import numpy as np
def IFFT(F, origin=1): N = len(F) n = N // 2
if N == 1: return F X_even = IFFT(F[0::2], 0) X_odd = IFFT(F[1::2], 0)
W = np.exp(2j * np.pi * np.arange(n) / N)
X = np.zeros(N, dtype='complex') X[:n] = X_even + W * X_odd X[n:] = X_even - W * X_odd
if origin: X /= N
return X
|
テスト
numpy.fft
モジュールを使って検算をしてみます.
テストコード
テスト1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
| print("# 入力配列") X = np.array([1, 4, 3, 2, 0, 8, 4, 7]) print(X)
print("\n# FFTの実行結果") F = FFT(X) print(F)
print("\n# IFFTの実行結果") X_ = IFFT(F) print(X_)
print("\n# numpyでのFFT") F_numpy = fft.fft(X) print(F_numpy)
print("\n# numpyでのIFFT") X_numpy = fft.ifft(F_numpy) print(X_numpy)
|
結果
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
| # 入力配列 [1 4 3 2 0 8 4 7]
# FFTの実行結果 [ 29. +0.j 1.70710678+7.36396103j -6. -3.j 0.29289322+5.36396103j -13. +0.j 0.29289322-5.36396103j -6. +3.j 1.70710678-7.36396103j]
# IFFTの実行結果 [1.+0.0000000e+00j 4.-4.5924255e-17j 3.+3.0616170e-17j 2.+4.5924255e-17j 0.+0.0000000e+00j 8.-4.5924255e-17j 4.-3.0616170e-17j 7.+4.5924255e-17j]
# numpyでのFFT [ 29. +0.j 1.70710678+7.36396103j -6. -3.j 0.29289322+5.36396103j -13. +0.j 0.29289322-5.36396103j -6. +3.j 1.70710678-7.36396103j]
# numpyでのIFFT [1.+0.j 4.+0.j 3.+0.j 2.+0.j 0.+0.j 8.+0.j 4.+0.j 7.+0.j]
|
けっこう良さそう!!