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.
