Skip to content

Instantly share code, notes, and snippets.

@lukicdarkoo
Created December 19, 2015 12:24
Show Gist options
  • Star 6 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save lukicdarkoo/1ab6a9a7b24025cb428a to your computer and use it in GitHub Desktop.
Save lukicdarkoo/1ab6a9a7b24025cb428a to your computer and use it in GitHub Desktop.
Basic implementation of Cooley-Tukey FFT algorithm in Python
'''
Basic implementation of Cooley-Tukey FFT algorithm in Python
Reference:
https://en.wikipedia.org/wiki/Fast_Fourier_transform
'''
__author__ = 'Darko Lukic'
__email__ = 'lukicdarkoo@gmail.com'
import numpy as np
import matplotlib.pyplot as plt
import random
SAMPLE_RATE = 8192
N = 128 # Windowing
def fft(x):
X = list()
for k in xrange(0, N):
window = 1 # np.sin(np.pi * (k+0.5)/N)**2
X.append(np.complex(x[k] * window, 0))
fft_rec(X)
return X
def fft_rec(X):
N = len(X)
if N <= 1:
return
even = np.array(X[0:N:2])
odd = np.array(X[1:N:2])
fft_rec(even)
fft_rec(odd)
for k in xrange(0, N/2):
t = np.exp(np.complex(0, -2 * np.pi * k / N)) * odd[k]
X[k] = even[k] + t
X[N/2 + k] = even[k] - t
x_values = np.arange(0, N, 1)
x = np.sin((2*np.pi*x_values / 32.0)) # 32 - 256Hz
x += np.sin((2*np.pi*x_values / 64.0)) # 64 - 128Hz
X = fft(x)
# Plotting
_, plots = plt.subplots(2)
## Plot in time domain
plots[0].plot(x)
## Plot in frequent domain
powers_all = np.abs(np.divide(X, N/2))
powers = powers_all[0:N/2]
frequencies = np.divide(np.multiply(SAMPLE_RATE, np.arange(0, N/2)), N)
plots[1].plot(frequencies, powers)
## Show plots
plt.show()
@YuePanEdward
Copy link

Really a good job

@skilz80
Copy link

skilz80 commented Sep 28, 2021

I tried implementing this in Python 3+ and I had to change xrange into range which is easy enough. However, when I try to run the script... I'm getting a TypeError message being raised that 'float' object cannot be interpreted as an integer for the code segment: for k in range(0, N/2): on line 42. I like your implementation, have you considered updating it to work with newer versions of Python? I'm currently using Python 3.8.7

@fpeterek
Copy link

I tried implementing this in Python 3+ and I had to change xrange into range which is easy enough. However, when I try to run the script... I'm getting a TypeError message being raised that 'float' object cannot be interpreted as an integer for the code segment: for k in range(0, N/2): on line 42. I like your implementation, have you considered updating it to work with newer versions of Python? I'm currently using Python 3.8.7

The / operator is interpreted as floating point division and returns a float regardless of it's arguments or result. You can use the // operator to perform an integer division. The // operator discards the remainder of the division and converts the result to an int. Alternatively, you can cast the result of the division to int yourself.

range(N//2)/range(int(N/2))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment