# Troubles with implementing a Cooley-Tukey style FFT in python

For education purposes I decided to try to implement a Cooley-Tukey style FFT algorithm in python, I must say that I am in no way a professional python dev nor a signal processing expert, so ugly code and approximate comprehension of the subject matter are to be somewhat expected.

Anyway, I followed the wikipedia article on the topic and re-demonstrated myself the necessary equations on paper from the original DFT equation.

I tried to implement this in Python, using both the equations and the provided pseudo-code. I then injected a sine wave into my FFT code and instead of getting the frequency of my sine wave on either side of the shannon frequency, I do get the expected frequency plus another unexpected frequency that takes up more energy out of my signal the more I increase my sample frequency.

For the sake of troubleshooting, I did only split my samples into even and odd indexes once.

Here's the sine wave generation code:

fs = 500           #sampling frequency
Ts = 1/fs          #sampling period
samples = 128      #number of samples
t = np.arange(0, samples * Ts, Ts)

f1 = 100           #signal frequency

signal = np.array([np.sin(2 * np.pi * f1 * t[i]) for i in range(0,samples)])


And here's my FFT implementation:

(In a nutshell, I start my code by making a few check ups and proceed to split my sample array of size N into two other arrays of size N/2 filled respectively with the even and odd indexes of the original sample array. I then apply a DFT on those and compute the twiddle factor array of size N/2. Eventually, I sum the DFT'd even array with the DFT'd odd array multiplied by the twiddle factor and do the exact same thing with a substraction in order to get the full-scale FFT.)

N = len(signal)

#pre-calculation check-ups
if N < 2:
print("Selected sample is too small (less than 2 points)")
else:
powerOfTwo = np.log2(N)
if powerOfTwo.is_integer() == False:
print("Selected sample isn't a power of 2")

#start of the actual algorithm
else:
FFT1 = []
FFT2 = []
FFT = []
evenArr = []
evenArrFFT = []
oddArr = []
oddArrFFT = []
TF = []

#sorts samples by parity
for i in np.arange(0, N):
if(i % 2 == 0):
evenArr = np.append(evenArr, signal[i])
else:
oddArr = np.append(oddArr, signal[i])

#applies DFT to the even and odd parts of the original sample array individually
evenArrFFT = DFT(evenArr, "blackman")
oddArrFFT = DFT(oddArr, "blackman")

#initializes the twiddle factor
for k in np.arange(0, int(N/2)):
TF = np.append(TF, np.exp(-1j * 2 * np.pi * k/N))

FFT1 = evenArrFFT + TF * oddArrFFT  #second half of the FFT
FFT2 = evenArrFFT - TF * oddArrFFT  #first half of the FFT
FFT = np.append(FFT2, FFT1)         #complete FFT


This is the spectrum obtained with my own DFT implementation (with a blackman window): And this is the spectrum obtained with Numpy's fft.fft (with a blackman window): However this is what I'm getting when I'm plotting the output of my FFT implementation (once again, you guessed it, with a blackman window): I did try to check if my split into even and odd indexes worked, if my DFT implementation worked, if I didn't make mistakes in initializing the twiddle factor, but I cannot seem to find the root of the issue and I'm starting to feel like that this issue is more likely to originate from a lack of understanding of signal processing rather than a programming error.

I am fully aware that more efficient implementations of the FFT exist on the internet, the point of this exercise being to code it myself in a way that I understand. Any help would be appreciated.

• What exactly is the difference between the first and last picture? Apr 18 at 12:54
• @Hilmar The first one is the spectrum computed by my DFT algorithm and displays the expected frequency of my signal. The last one is the spectrum computed by my FFT algorithm, which uses the aforementionned DFT algorithm, and displays my signal's frequency as well as that parasitic frequency my post is about. Apr 18 at 13:01
• Are you just trying to implement a single stage of the FFT using another DFT implementation for the sub-transforms? If yes, don't use a window. Apr 18 at 14:08