2016-02-01 2 views
0

Im пытается переделать функции fft в python. Здесь Iv'e видел аналогичный вопрос 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 все еще немного странно, не я. –

Смежные вопросы