pythonによるFFTの実装
pythonでFFT(高速フーリエ変換)を実装しようと思っています
コードはご覧の通りです
(FFT_sort.py)
import numpy as np
def sort(N):
flag = ~(N & (N - 1))
if flag != -1:
return None
result = np.zeros(N, dtype=np.int64)
result[0] = 0
result[1] = N / 2
result[N / 2] = 1
result[N / 2 + 1] = N / 2 + 1
T = n = N / 2
oldN = n
oldCat = 0
cat = 1
count = 1
judge = 0
next = 2
linspace = np.array([i for i in range(N)])
while(n != 2):
n /= 2
cat *= 2
judge += cat
normal = 2 ** (count + 1))
for i in range(cat / 2):
j = 1 + i * 2
index = j * N / normal
result[index] = next
result[index + 1] = next + T
result[index + T] = next + 1
result[index + T + 1] = next + T + 1
next += 2
count += 1
oldN = n
return result
(FFT.py)
import FFT_sort as fs
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
N = input()
x = np.linspace(0, 2 * np.pi, N)
y = np.sin(x)
data = fs.sort(N)
if data is None:
print("error")
quit()
else:
print(data)
result = np.zeros(N, dtype=np.complex)
result[:] = y[:]
num = 1
while(num != N):
count = 0
T = num
F = N / (num * 2)
print("F:", F)
for i in range(N):
judge = i / num
if judge % 2 == 0:
k = i + T
else:
k = i - T
iF = i * F % N
w = np.exp(float(i) * iF * -1j * 2 * pi / float(N))
print("iF:", iF)
print("w:", w)
n = data[k]
result[i] += w * result[n]
count += 1
num *= 2
print(result)
answer = np.fft.fft(y)
print(answer)
コードの入力はFFTのサンプリングの数で
アウトプットは
私のFFTの結果とnumpyでやった正規のFFTの結果です
ここで質問ですが
私のコードでやったFFTと
本家numpyのFFTの結果が全然違います
[ 6.23785596 -4.42077404j 5.0713722 -3.1314097 j
-3.8474327 -1.99237623j 2.05949337 -0.61485772j
5.84009223 -4.84151225j -2.02845113 -0.82589246j
-1.36232081 -8.6730753 j 1.26764366 -7.27397125j
-0.01871131 -4.59613684j 1.0188491 -4.96565877j
3.58757676-10.1782615 j 2.92161011 -6.39504849j
-2.59079252 -3.55000541j 2.78417111-13.26650134j
-6.0297906 -1.17150871j -4.92477154 -4.98328949j]
[-2.49800181e-16+0. j 1.49800460e+00-7.53097769j
-2.88537029e-01+0.69659001j -2.36488255e-01+0.35392969j
-2.22614343e-01+0.22261434j -2.16932373e-01+0.14494958j
-2.14217111e-01+0.08873163j -2.12937210e-01+0.04235584j
-2.12556562e-01+0. j -2.12937210e-01-0.04235584j
-2.14217111e-01-0.08873163j -2.16932373e-01-0.14494958j
-2.22614343e-01-0.22261434j -2.36488255e-01-0.35392969j
-2.88537029e-01-0.69659001j 1.49800460e+00+7.53097769j]
このような結果になるのは何故でしょうか
コードのどの部分に問題があるのでしょうか?
上の括弧がサンプル点が8つと仮定した時の結果が私のFFTの結果で
下の括弧が本家のnp.fft.fftのアルゴリズムの結果です
ご覧の通り、結果は全然違います
自分では何回も検査しましたが、なかなかみつかりません
バタフライ演算の方向もあっていますし
wの値のとり方に問題がありそうですが
自力では見つかりません
理由を教えて貰えないでしょうか?