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.