Recursive FFT in Python Convertible to Verilog/VHDL

This code snippet is a recursive FFT implementation in Python.  This method is used to generate Verilog and VHDL with MyHDL/Python.

Since the FFT is a very popular tool and there are many efficient implementations (see FFTW www.fftw.org) this code is not intended to be a high performance implementation. Rather, this code snippet is more useful as an exercise to gain more understanding of the FFT and for hardware generation.

The code can also be used as a reference for implementing processor specific implementations. Or as an algorithm exploration tool.

Convertible to Verilog and VHDL

This code was used to develop the HDL version here.

Executing the Code

The code can easily be run from the command line.

>>  python rfft.py

The above will execute a simple testbench and verify the recursive FFT by comparing it to the direct implementation of the DFT. If the module is executed as described above, by default, it will run a fairly exhaustive test. Depending on the machine this will take some time to execute.

#
# rfft.py
#
# This file contains a recursive version of the fast-fourier transform and
# support test functions.  This module utilizes the numpy (numpy.scipy.org)
# library.  
#
# References
#   - http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
#   - "A Simple and Efficient FFT Implementation in C++", by Vlodymyr Myrnyy

import numpy
from numpy.fft import fft
from numpy import sin, cos, pi, ones, zeros, arange, r_, sqrt, mean

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def rFFT(x):
    """
    Recursive FFT implementation.

    References
      -- http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
      -- "A Simple and Efficient FFT Implementation in C++"
          by Vlodymyr Myrnyy
    """
    
    n = len(x)

    if (n == 1):
	return x

    w = getTwiddle(n)
    m = n/2;
    X = ones(m, float)*1j
    Y = ones(m, float)*1j
    
    for k in range(m):
        X[k] = x[2*k]
        Y[k] = x[2*k + 1] 

    X = rFFT(X)  
    Y = rFFT(Y) 

    F = ones(n, float)*1j
    for k in range(n):
        i = (k%m)
        F[k] = X[i] + w[k] * Y[i]

    return F

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def getTwiddle(NFFT=8):
    """Generate the twiddle factors"""

    W = r_[[1.0 + 1.0j]*NFFT]

    for k in range(NFFT):
        W[k] = cos(2.0*pi*k/NFFT) - 1.0j*sin(2.0*pi*k/NFFT)

    return W

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def DFT(x, N=8):
    """
    Use the direct definition of DFT for verification
    """
    y = [1.0 + 1.0j]*N
    y = r_[y]
    for n in range(N):
        wsum = 0 + 0j;
	for k in range(N):
	    wsum = wsum + (cos(2*pi*k*n/N) - (1.0j * sin(2*pi*k*n/N)))*x[k]
        
	y[n] = wsum
        
    return y

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def test_rfft(N      = 64,      # FFT order to test
              nStart = 0.2,     # Note aliased signal included
              nStep  = 2.1,     # Samples per period step
              pStep  = pi/4,    # Phase step size
              limErr = 10e-12,  # Error limit to check
              maxErr = 0        # Max difference
              ):
    """
    Use the built in numpy FFT functions and the direct
    implemenation of the DFT to verify the recursive FFT.

    This testbench verifies the different implementations are within
    a certain limit.  Because of the different implemenations the values
    could be slightly off (computer representation calculation error).
    """

    # Use test signal nStart:nStep:N samples per cycle
    for s in arange(nStart, N+nStep, nStep):
        for p in arange(0, pi+pStep, pStep):

            n = arange(N, 0, -1)
            x = cos(2*pi*n/s + p)

            xDFT = DFT(x,N)
            nFFT = fft(x,N)
            xFFT = rFFT(x)

            rmsErrD = sqrt(mean(abs(xDFT - xFFT))**2)
            rmsErrN = sqrt(mean(abs(nFFT - xFFT))**2)

            if rmsErrD > limErr or rmsErrN > limErr:
                print s, p, "Error!", rmsErrD, rmsErrN
                print xDFT
                print nFFT
                print xFFT

            if rmsErrD > maxErr:
                maxErr = rmsErrD
            elif rmsErrN > maxErr:
                maxErr = rmsErrN

    print "N %d maxErr = %f " % (N,maxErr)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If the module is run test a bunch of different size FFTs
if __name__ == '__main__':

    # The following is fairly exhaustive and will take some time
    # to run.
    tv = 2**arange(1,12)
    for nfft in tv:
        test_rfft(N=nfft)