1
$\begingroup$

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):

DFT

And this is the spectrum obtained with Numpy's fft.fft (with a blackman window):

np.fft.fft

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):

my fft implementation

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.

$\endgroup$
3
  • $\begingroup$ What exactly is the difference between the first and last picture? $\endgroup$
    – Hilmar
    Apr 18, 2023 at 12:54
  • $\begingroup$ @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. $\endgroup$
    – tampler
    Apr 18, 2023 at 13:01
  • $\begingroup$ 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. $\endgroup$
    – Hilmar
    Apr 18, 2023 at 14:08

1 Answer 1

2
$\begingroup$

The use of the window is definitely incorrect here. If you want to window, apply the window BEFORE doing the FFT and don't use any windows in the internal DFTs.

$\endgroup$
3
  • $\begingroup$ It may be incorrect and I'll correct that, however it does not change the core of the issue. After trying it, removing the blackman window does increase the amplitude of secondary lobes a little bit as expected but the unexpected frequency still remains. $\endgroup$
    – tampler
    Apr 19, 2023 at 16:09
  • $\begingroup$ @tampler Then edit your question to omit any windowing discussions because that's completely separate. We prefer focused questions, and this'd be closed as an XY problem on other networks. $\endgroup$ Apr 22, 2023 at 9:15
  • $\begingroup$ @tampler: can you take out all the windowing and post the results? By visual inspection the rest of your code looks good (other than the last line where the order of the cascade seems reversed). $\endgroup$
    – Hilmar
    Apr 23, 2023 at 12:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.