implemented in-place FFT c-style. (by ChatGPT)
This commit is contained in:
171
FFT.py
171
FFT.py
@ -1,5 +1,5 @@
|
|||||||
import numpy as np
|
import numpy as np
|
||||||
|
from math import sin, cos, pi, sqrt
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -12,17 +12,38 @@ def FFT_backup(inp):
|
|||||||
|
|
||||||
|
|
||||||
def FFT_real(inp):
|
def FFT_real(inp):
|
||||||
inp = np.array(inp)
|
|
||||||
out = FFT(inp)
|
mode = 2
|
||||||
out = [val for val in np.abs(out)]
|
if (mode == 0):
|
||||||
|
inp = np.array(inp)
|
||||||
|
out = FFT_np(inp)
|
||||||
|
out = [np.abs(val) for val in out]
|
||||||
|
if (mode == 1):
|
||||||
|
|
||||||
|
out = FFT_arr([[val, 0] for val in inp])
|
||||||
|
#out = [val for val in np.abs(out)]
|
||||||
|
print("FFT_calculated!")
|
||||||
|
print("output:", out)
|
||||||
|
out = [sqrt(val[0]**2 + val[1]**2) for val in out]
|
||||||
|
# out = [2 for val in out]
|
||||||
|
if (mode == 2):
|
||||||
|
|
||||||
|
|
||||||
|
out = FFT_arr_inplace([[val, 0] for val in inp])
|
||||||
|
#out = [val for val in np.abs(out)]
|
||||||
|
print("FFT_calculated!")
|
||||||
|
print("output:", out)
|
||||||
|
out = [sqrt(val[0]**2 + val[1]**2) for val in out]
|
||||||
|
|
||||||
print(out)
|
print(out)
|
||||||
return out
|
return out
|
||||||
|
|
||||||
|
|
||||||
def FFT(inp):
|
def FFT(inp):
|
||||||
return FFT_np(inp)
|
return FFT_arr([[val, 0] for val in inp])
|
||||||
|
|
||||||
def FFT_np(inp):
|
def FFT_np(inp):
|
||||||
|
print("inp:", inp)
|
||||||
inp = np.array(inp)
|
inp = np.array(inp)
|
||||||
N = inp.shape[0]
|
N = inp.shape[0]
|
||||||
if N == 1:
|
if N == 1:
|
||||||
@ -33,3 +54,143 @@ def FFT_np(inp):
|
|||||||
tw = np.exp(-2j * np.pi * k / N) * X_odd
|
tw = np.exp(-2j * np.pi * k / N) * X_odd
|
||||||
return np.concatenate((X_even + tw, X_even - tw))
|
return np.concatenate((X_even + tw, X_even - tw))
|
||||||
|
|
||||||
|
def FFT_arr(x):
|
||||||
|
print("x:", x)
|
||||||
|
|
||||||
|
N = len(x)
|
||||||
|
print("len(x):", N)
|
||||||
|
if N == 1:
|
||||||
|
print("As len(x) == 1, return it`s value")
|
||||||
|
return x
|
||||||
|
X_even = FFT_arr(x[::2])
|
||||||
|
X_odd = FFT_arr(x[1::2])
|
||||||
|
print("X_even:",X_even)
|
||||||
|
|
||||||
|
tw = []
|
||||||
|
for k in range(len(X_odd)):
|
||||||
|
a = cos(2*pi * k/N) #real
|
||||||
|
b = -sin(2*pi * k/N) #imag
|
||||||
|
|
||||||
|
c = X_odd[k][0] # real
|
||||||
|
d = X_odd[k][1] # imag
|
||||||
|
|
||||||
|
print("a,b,c,d:", a,b,c,d)
|
||||||
|
|
||||||
|
tw.append([a*c - b*d, b*c + a*d]) #(a + ib)(c + id) = (ac - bd) + i(bc + ad)
|
||||||
|
res = [0 for i in range(len(X_even)*2)]
|
||||||
|
for i in range(len(X_even)):
|
||||||
|
res[i] = [X_even[i][0] + tw[i][0], X_even[i][1] + tw[i][1]]
|
||||||
|
res[i+len(X_even)] = [X_even[i][0] - tw[i][0], X_even[i][1] - tw[i][1]]
|
||||||
|
return res
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
import math
|
||||||
|
|
||||||
|
def FFT_arr_inplace(buf):
|
||||||
|
"""
|
||||||
|
In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
|
||||||
|
Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
|
||||||
|
"""
|
||||||
|
N = len(buf)
|
||||||
|
# --- bit-reverse перестановка (чтобы бабочки шли последовательно) ---
|
||||||
|
j = 0
|
||||||
|
for i in range(1, N):
|
||||||
|
bit = N >> 1
|
||||||
|
while j & bit:
|
||||||
|
j ^= bit
|
||||||
|
bit >>= 1
|
||||||
|
j |= bit
|
||||||
|
if i < j:
|
||||||
|
buf[i], buf[j] = buf[j], buf[i]
|
||||||
|
|
||||||
|
# --- уровни бабочек ---
|
||||||
|
m = 2
|
||||||
|
while m <= N:
|
||||||
|
half = m // 2
|
||||||
|
# шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw)
|
||||||
|
cΔ = math.cos(2.0 * math.pi / m)
|
||||||
|
sΔ = -math.sin(2.0 * math.pi / m)
|
||||||
|
|
||||||
|
for start in range(0, N, m):
|
||||||
|
# w = e^{-j*0*Δ} = 1 + j0
|
||||||
|
wr, wi = 1.0, 0.0
|
||||||
|
for k in range(half):
|
||||||
|
u_re, u_im = buf[start + k]
|
||||||
|
v_re, v_im = buf[start + k + half]
|
||||||
|
|
||||||
|
# t = w * v
|
||||||
|
t_re = wr * v_re - wi * v_im
|
||||||
|
t_im = wr * v_im + wi * v_re
|
||||||
|
|
||||||
|
# верх/низ
|
||||||
|
buf[start + k][0] = u_re + t_re
|
||||||
|
buf[start + k][1] = u_im + t_im
|
||||||
|
buf[start + k + half][0] = u_re - t_re
|
||||||
|
buf[start + k + half][1] = u_im - t_im
|
||||||
|
|
||||||
|
# w *= w_m (поворот на Δ с помощью рекуррентной формулы)
|
||||||
|
# (wr + j wi) * (cΔ + j sΔ)
|
||||||
|
wr, wi = wr * cΔ - wi * sΔ, wr * sΔ + wi * cΔ
|
||||||
|
|
||||||
|
m <<= 1
|
||||||
|
return buf
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def FFT_arr_for_C(x):
|
||||||
|
print("x:", x)
|
||||||
|
|
||||||
|
N = len(x)
|
||||||
|
print("len(x):", N)
|
||||||
|
if N == 1:
|
||||||
|
print("As len(x) == 1, return it`s value")
|
||||||
|
return x
|
||||||
|
X_even = FFT_arr_for_C(x[::2])
|
||||||
|
X_odd = FFT_arr_for_C(x[1::2])
|
||||||
|
print("X_even:",X_even)
|
||||||
|
|
||||||
|
tw = []
|
||||||
|
for k in range(len(X_odd)):
|
||||||
|
a = cos(2*pi * k/N) #real
|
||||||
|
b = -sin(2*pi * k/N) #imag
|
||||||
|
|
||||||
|
c = X_odd[k][0] # real
|
||||||
|
d = X_odd[k][1] # imag
|
||||||
|
|
||||||
|
print("a,b,c,d:", a,b,c,d)
|
||||||
|
|
||||||
|
tw.append([a*c - b*d, b*c + a*d]) #(a + ib)(c + id) = (ac - bd) + i(bc + ad)
|
||||||
|
res = [0 for i in range(len(X_even)*2)]
|
||||||
|
for i in range(len(X_even)):
|
||||||
|
res[i] = [X_even[i][0] + tw[i][0], X_even[i][1] + tw[i][1]]
|
||||||
|
res[i+len(X_even)] = [X_even[i][0] - tw[i][0], X_even[i][1] - tw[i][1]]
|
||||||
|
return res
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
Binary file not shown.
@ -37,7 +37,7 @@ def DFT_naive(inp, out):
|
|||||||
val_im += inp[n] * cos(phi) /INP_L
|
val_im += inp[n] * cos(phi) /INP_L
|
||||||
|
|
||||||
val_abs = abs_f(val_re, val_im)
|
val_abs = abs_f(val_re, val_im)
|
||||||
print("F, val_abs:",f, val_abs)
|
#print("F, val_abs:",f, val_abs)
|
||||||
out[f] = val_abs
|
out[f] = val_abs
|
||||||
|
|
||||||
|
|
||||||
@ -148,7 +148,7 @@ def DFT_naive_FP(inp_float, out):
|
|||||||
val_im += inp[n] * phi_cos /INP_L
|
val_im += inp[n] * phi_cos /INP_L
|
||||||
|
|
||||||
val_abs = abs_FP(val_re, val_im)
|
val_abs = abs_FP(val_re, val_im)
|
||||||
print("F, val_abs:",f, val_abs)
|
#print("F, val_abs:",f, val_abs)
|
||||||
out[f] = val_abs
|
out[f] = val_abs
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user