fft.py
51 lines
| 1020 B
| text/x-python
|
PythonLexer
/ SRC / fft.py
Alexis Jeandet
|
r0 | import math | |
import numpy as np | |||
from bin16 import * | |||
def power_of_2(number): | |||
if number == 0: | |||
return 0 | |||
while(number % 2 == 0): | |||
number = number / 2 | |||
if number > 1: | |||
return 0 | |||
return 1 | |||
#here we simply do the FFT algorithm, using the Cooley-Tukey idea of divide and conquer and quantizing the values (working with integers) | |||
#note that the output result will be already divided by N, which is not what happens for the np.fft routine | |||
def fft_CT(x): | |||
N = len(x) | |||
if N == 1: | |||
return x | |||
x_e = x[::2] | |||
x_o = x[1::2] | |||
X_e = fft_CT(x_e) | |||
X_o = fft_CT(x_o) | |||
X = np.zeros(N,'complex') | |||
M = N/2 | |||
for k in range(M): | |||
X[k] = X_e[k] + X_o[k] * pow(math.e,-2j*math.pi*k/N) | |||
a = X[k].real/2 | |||
b = X[k].imag/2 | |||
a = int(a) | |||
a = trunc16(a) | |||
b = int(b) | |||
b = trunc16(b) | |||
c = complex(a,b) | |||
X[k] = c | |||
for k in range(M,N): | |||
X[k] = X_e[k-M] - X_o[k-M] * pow(math.e,-2j*math.pi*(k-M)/N) | |||
a = X[k].real/2 | |||
b = X[k].imag/2 | |||
a = int(a) | |||
a = trunc16(a) | |||
b = int(b) | |||
b = trunc16(b) | |||
c = complex(a,b) | |||
X[k] = c | |||
return X | |||