高速フーリエ変換（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]
|
けっこう良さそう！！