FFT_and_FP_math/FFT.py

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 = 3
  if (mode == 0):
    inp = np.array(inp)
    out = FFT_np(inp)
    out = [np.abs(val) for val in out]
  if (mode == 1): # simple with recursion
    
    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): #no internal buffs
  
    tmp = []
    for re in inp:
      tmp.append(re)
      tmp.append(0)
      
    out = FFT_arr_inplace(tmp)
    
    tmp = out.copy()
    out = []
    for i in range(len(tmp)//2):
      out.append([tmp[2*i], tmp[2*i+1]])
    
    
    #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]
    
  if (mode == 3): #fixed-point
  
    tmp = []
    for re in inp:
      tmp.append(re)
      tmp.append(0)
      
    out = FFT_arr_FP(tmp)
    
    tmp = out.copy()
    out = []
    for i in range(len(tmp)//2):
      out.append([tmp[2*i], tmp[2*i+1]])
    
    
    #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




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_delta = cos(2.0 * math.pi / m)
    s_delta = -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
        print("wr, v_re, wi, v_im:",wr, v_re, wi, v_im)
        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_delta - wi * s_delta, wr * s_delta + wi * c_delta

    m <<= 1
  return buf










FP_acc = 1e3

def FFT_arr_FP(buf):
  """
    In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
    Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
    """
    
  buf = [int(val*FP_acc) for val in buf]

  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_delta_FP = cos(2.0 * pi / m) * FP_acc
    s_delta_FP = -sin(2.0 * pi / m) * FP_acc

    for start in range(0, N, m):
      # w = e^{-j*0*Δ} = 1 + j0
      wr, wi = 1.0 * FP_acc, 0.0 * FP_acc
      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
        print("wr, v_re, wi, v_im:", wr, v_re, wi, v_im)
        t_re = (wr * v_re - wi * v_im) / FP_acc
        t_im = (wr * v_im + wi * v_re) /  FP_acc

        # верх/низ
        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_delta_FP - wi * s_delta_FP)/ FP_acc, (wr * s_delta_FP + wi * c_delta_FP)/ FP_acc
        

    m <<= 1
    
  buf = [val/ FP_acc for val in buf] 
  return buf