diff --git a/FFT.py b/FFT.py index 25e6b14..4620b9a 100644 --- a/FFT.py +++ b/FFT.py @@ -1,5 +1,5 @@ import numpy as np - +from math import sin, cos, pi, sqrt @@ -12,17 +12,38 @@ def FFT_backup(inp): def FFT_real(inp): - inp = np.array(inp) - out = FFT(inp) - out = [val for val in np.abs(out)] + + mode = 2 + 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) return out def FFT(inp): - return FFT_np(inp) + return FFT_arr([[val, 0] for val in inp]) def FFT_np(inp): + print("inp:", inp) inp = np.array(inp) N = inp.shape[0] if N == 1: @@ -33,3 +54,143 @@ def FFT_np(inp): tw = np.exp(-2j * np.pi * k / N) * X_odd 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 + + + + + + + + + \ No newline at end of file diff --git a/__pycache__/FFT.cpython-312.pyc b/__pycache__/FFT.cpython-312.pyc index 7ed90aa..c7e41bb 100644 Binary files a/__pycache__/FFT.cpython-312.pyc and b/__pycache__/FFT.cpython-312.pyc differ diff --git a/naive_DFT.py b/naive_DFT.py index f00483e..6b3b917 100755 --- a/naive_DFT.py +++ b/naive_DFT.py @@ -37,7 +37,7 @@ def DFT_naive(inp, out): val_im += inp[n] * cos(phi) /INP_L 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 @@ -148,7 +148,7 @@ def DFT_naive_FP(inp_float, out): val_im += inp[n] * phi_cos /INP_L 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