FFT App Implementation (WIP)

[bholt]

A simple 1D FFT on a massive global array, leaving the result in the original array.

Implementation details

• Cooley-Tukey algorithm
  1. Reinterpret 1D array as 2D matrix: N1xN2
  2. Perform N1 DFTs of size N2
  3. Multiply by twiddle factors
  4. Perform N2 DFTs of size N1
• Thoughts
  • Maybe N1 or N2 is sizeof(complex<double>)/block_size?
  • Use FFTW when decomposition is small enough for elements to fit in a single block

[bholt]

FFT from RosettaCode:

from cmath import exp, pi

def fft(x):
    N = len(x)
    if N <= 1: return x
    even = fft(x[0::2])
    odd =  fft(x[1::2])
    return [even[k] + exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] + \
           [even[k] - exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]

print fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])

[bholt]

Iterative version, roughly interpretted from Introduction to Parallel Computing:

def iterative_fft(X, Y, n):
  r = log2(n)
  omega = exp(2j*pi/n)
  # bitrev(i,r): flips the lowest 'r' bits of 'i'

  for i in range(n):
    R[i] = X[i]
  for m in range(r): # outer loop
    for i in range(n):
      S[i] = R[i]
    for i in range(n): # inner loop
      rm = (r-1) - m
      j = i & ~(1 << rm) # i with m'th bit set to 0
      k = i | (1 << rm) # i with m'th bit set to 1
      R[i] = S[j] + S[k] * omega**( (bitrev(i,r) << rm) & ((1<<r)-1) )
  for i in range(n):
    Y[i] = R[i]