Thursday, February 26, 2015

py28. Fast Fourier Transform in Python

The Discrete Fourier Transform (DFT) is the pair of transforms for the transform and its inverse:


The Fast Fourier Transform (FFT) is the solution of DFT using an algorithm based on symmetry of equations.


We use the fft function from the numpy.fft module to solve for the transform of a 64-length numpy array. It returns a complex numpy array, dtype = 'complex', which is sent to ifft function in the same module. We compare it to the the result from dft and inverse dft. The result are the same, with errors on order of double-precision tolerance.

# ex28.py

from __future__ import print_function, division
import numpy as np
from numpy.fft import fft, ifft
from numpy import pi, exp
def dft(x, inverse = False):
    N = len(x)
    s,d = (1,N) if inverse else (-1,1)
    X = np.zeros(N, dtype='complex')
    for k in range(N):
        for n in range(N):
            X[k] += x[n]*exp(1j*s*2*pi*k*n/N)
    return X/d

x = np.array(range(32)+range(31,-1,-1))
y = fft(x)
y1 = dft(x)
z = ifft(y)
z1 = dft(y1,True)
xR = z.real
x1R = z1.real
error_fft = (x-xR).std()
error_dft = (x-x1R).std()
print("error_fft = %e\nerror_dft = %e" % (error_fft,error_dft))
# error_fft = 1.122558e-15
# error_dft = 1.060977e-13

No comments:

Post a Comment