Updated on 

高速フーリエ変換(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

  1. 入力配列を偶数行、奇数行に分割→x_evenx_odd
  2. 偶数行、奇数行のバタフライ演算を再帰的に計算→F_evenF_odd
  3. 出力配列を用意→F
  4. 回転因子の配列を用意→W
  5. 出力配列の前半分をF_even + W*F_oddとする
  6. 出力配列の後半分をF_even - W*F_odd とする

↓実際のコード(参考サイト

FFT
1
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

# 再帰的にFFTを計算
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とそんなに変化なし. 回転因子の の符号を逆にして,最後に をかけるだけです.

IFFT
1
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

# 再帰的にFFTを計算
X_even = IFFT(F[0::2], 0)
X_odd = IFFT(F[1::2], 0)

# 回転因子(IFFTなので 2*pi*i)
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]

けっこう良さそう!!


このページはHexoStellarを使用して作成されました。