import numpy as np from math import sin, cos, pi, sqrt def FFT_backup(inp): inp = np.array(inp) out = np.fft.fft(inp) out = [val for val in np.abs(out)] print(out) return out def FFT_real(inp): 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_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: return inp X_even = FFT_np(inp[::2]) X_odd = FFT_np(inp[1::2]) k = np.arange(N // 2) 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)/2 # --- 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*2], buf[i*2+1], buf[j*2], buf[j*2+1] = buf[j*2], buf[j*2+1], buf[i*2], buf[i*2+1] # --- уровни бабочек --- 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)*2], buf[(start + k)*2 +1] v_re, v_im = buf[(start + k + half)*2], buf[(start + k + half)*2 +1] # t = w * v t_re = wr * v_re - wi * v_im t_im = wr * v_im + wi * v_re # верх/низ buf[(start + k)*2 +0] = u_re + t_re buf[(start + k)*2 +1] = u_im + t_im buf[(start + k + half)*2 + 0] = u_re - t_re buf[(start + k + half)*2 + 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