# importing cmath makes exp and other math functions handle # complex numbers from cmath import * def fft(h): """compute the discrete Fourier transform of the sequence h. Assumes that len(h) is a power of two. """ N = len(h) # the Fourier transform of a single value is itself if N == 1: return h # recursively compute the FFT of the even and odd values He = fft(h[0:N:2]) Ho = fft(h[1:N:2]) # merge the half-FFTs i = complex(0,1) W = exp(2*pi*i/N) ws = [pow(W,k) for k in range(N)] H = [e + w*o for w, e, o in zip(ws, He+He, Ho+Ho)] return H # make a signal with two sine components, f=6 and f=12 N = 128 t = [1.0*n/N for n in range(N)] h = [sin(2*pi*6*tn) + sin(2*pi*12*tn) for tn in t] # compute the Fourier transform H = fft(h) # print the spectral density function sdf = [Hn * Hn.conjugate() for Hn in H] for n in range(N/2 + 1): print '%d %.3f' % (n, sdf[n].real)