Created
December 19, 2015 12:24
-
-
Save lukicdarkoo/1ab6a9a7b24025cb428a to your computer and use it in GitHub Desktop.
Basic implementation of Cooley-Tukey FFT algorithm in Python
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
''' | |
Basic implementation of Cooley-Tukey FFT algorithm in Python | |
Reference: | |
https://en.wikipedia.org/wiki/Fast_Fourier_transform | |
''' | |
__author__ = 'Darko Lukic' | |
__email__ = '[email protected]' | |
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() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The
/
operator is interpreted as floating point division and returns afloat
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 anint
. Alternatively, you can cast the result of the division toint
yourself.range(N//2)
/range(int(N/2))