created C folder with FP_math.c, Makefile and plotter. FP_math have a FFT realisation with fixed-point math, and its tester

FFT.py

@@ -179,7 +179,7 @@ def FFT_arr_inplace(buf):
 
 
 
-FP_acc = 1e3
+FP_acc = 2<<16
 
 def FFT_arr_FP(buf):
   """
@@ -202,12 +202,19 @@ def FFT_arr_FP(buf):
       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]
 
   # --- уровни бабочек ---
+  
+  c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
+  s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
   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
+    print("N, m:", N,m)
+    
+    c_delta_FP = c_delta_FP_arr[m -2]
+    s_delta_FP = s_delta_FP_arr[m -2]
+    #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
@@ -238,7 +245,74 @@ def FFT_arr_FP(buf):
   buf = [val/ FP_acc for val in buf] 
   return buf
 
-
+#global arrays. calculated once at compile time. They store sin and cos values as fixed-point
+  c_delta_FP_arr = [int(cos(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)]
+  s_delta_FP_arr = [int(-sin(2.0 * pi / m) * FP_acc) for m in range(2, inp_L//2 + 2)]
+
+def FFT_arr_FP_for_C(inp, inp_L, buf): #version for translation to C directly
+  """
+    In-place radix-2 DIT FFT для списка buf длины N=2^m, где каждый элемент — [re, im].
+    Без массивов twiddle: твиддл на уровне обновляется рекуррентно.
+    """
+    
+  #buf = [int(val*FP_acc) for val in buf]
+  for i in range(inp_L):
+    buf[i*2] = inp[i] #take data from input as real and store it to the buf as Re, Im pairs.
+    buf[i*2 + 1] = 0
+
+  N = inp_L//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]
+
+  # --- уровни бабочек ---
+  
+  c_delta_FP_arr = [cos(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
+  s_delta_FP_arr = [-sin(2.0 * pi / m) * FP_acc for m in range(2,N+2)]
+  m = 2
+  while m <= N:
+    half = m // 2
+    # шаг угла Δ = 2π/m, базовый твиддл w_m = e^{-jΔ} => (cw, sw)
+    #print("N, m:", N,m)
+    
+    c_delta_FP = c_delta_FP_arr[m -2]
+    s_delta_FP = s_delta_FP_arr[m -2]
+    #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
+