For the FFT sonification I needed to implement my own FFT in order to keep track of which notes are being played while the Fourier transform executes. It turns out that implementing an FFT isn’t too difficult as long as you don’t care how fast it goes (Note: if you do need it to be fast you’ll have a hard time beating fftw which “is an implementation of the discrete Fourier transform (DFT) that adapts to the hardware in order to maximize performance” (cf http://www.fftw.org/fftw-paper-ieee.pdf) ).

The key ingredient to implementing an FFT with the Cooley-Tukey algorithm (which I used for the sonification) is a traversal of the input array in bit-reversed order. Here’s how I did the traversal in Python:

def bit_reverse_traverse(a):
    n = a.shape[0]
    assert(not n&(n-1) ) # assert that n is a power of 2

    if n == 1:
        yield a[0]
    else:
        even_index = arange(n/2)*2
        odd_index = arange(n/2)*2 + 1
        for even in bit_reverse_traverse(a[even_index]):
            yield even
        for odd in bit_reverse_traverse(a[odd_index]):
            yield odd

What I thought was interesting about this code was that I had to use Python’s yield statement, which I had never run into before.

In short, Python’s yield statement is useful when you need to perform recursion efficiently. The difference between a yield statement and a return statement is that yield produces a generator which executes your code on-the-fly rather than storing the entire recursive tree in memory (cf http://stackoverflow.com/a/231855/424631 for more). So, instead of having a global list to write the output of our recursive calls we use a yield statment to make a generator that runs through our list in bit reversed order. Doing something like z=bit_reverse_traverse(a) and then several z.next()s will go through a in the bit-reversed way without eating up all the memory needed for a deep recursive tree. Here’s the code to use the generator:

def get_bit_reversed_list(l):
    n = len(l)

    indexs = arange(n)
    b = []
    for i in bit_reverse_traverse(indexs):
        b.append(l[i])

    return b

For reference, the alternative recursive implementation is below:

def bit_reverse_traverse_no_generator(a):
    n = a.shape[0]
    assert(not n&(n-1))

    if n == 1:
        return a
    else:
        even_indicies = arange(n/2)*2
        odd_indicies = arange(n/2)*2 + 1

        evens = bit_reverse_traverse_no_generator(a[even_indicies])
        odds = bit_reverse_traverse_no_generator(a[odd_indicies])

        return concatenate([evens, odds])

def get_bit_reversed_list_no_generator(l):
    n = len(l)

    indexs = arange(n)
    b = []
    for i in bit_reverse_traverse_no_generator(indexs):
        b.append(l[i])

    return b