2016-02-01 84 views
0

我試圖重新編譯python中的fft函數。我在這裏看到過一個類似的問題Manual fft not giving me same results as fft,但是我很難看到我是否做了同樣的錯誤或不同的錯誤。python manual fft botched

import numpy as np 
import numpy.random as npr 

N=9 ### 10 -1 
MC=10 

###Genrate soem data 
data=complex(1,0)*npr.uniform(size=(N,MC))+complex(0,1)*npr.uniform(size=(N,MC)) 

naive_fft=complex(1,0)*np.zeros((N,MC)) 
for K in range(N): 
    for m in range(N): 
     phase=(2*np.pi*K*m)/float(N+1) 
     naive_fft[K,:]=naive_fft[K,:]+data[m,:]*np.exp(complex(0,1)*phase) 

fft=np.fft.fft(data,axis=0) 
ifft=np.fft.ifft(data,axis=0) 
print('fft') 
print(naive_fft-fft) 
print('ifft') 
print(naive_fft-ifft*(N+1.0)) 

我的成績相較於numpy的FFT我無法重現既不FFT也不IFFT(只有naive_fft[0,:]似乎匹配fft[0,:]值。

回答

2

有幾件事情就更不用說了。首先,在Python我們使用1j來表示虛數單位,而不是complex(0, 1)。如果您想將結果與numpy進行比較,那麼您必須檢查numpy是如何實現fft的。請參閱Numpy FFT docs以瞭解詳細信息。您會發現numpy遵循最常見的fft定義,它使用負指數。此外,在你的階段float(N+1)是完全錯誤的。必須閱讀N。 總而言之,你必須:

# ... 
naive_fft = np.zeros((N,MC), dtype='complex') 
for K in range(N): 
    for m in range(N): 
     phase=(-2*np.pi*K*m)/float(N) 
     naive_fft[K] += data[m] * np.exp(phase*1j) 

xfft = np.fft.fft(data, axis=0) 
# ... 

>>> np.isclose(xfft, naive_fft).all() 
    True 

測試它的逆變換工作類似,但以積極的指數。

+0

謝謝。你救了我。那個範圍(N)給出了0,....,N-1仍然有點奇怪,沒有我。 –