now printf is disabled if twiddles arrays are external. Now while used for emdedded there is no prfintf garbage in compiled code.

.project

@@ -0,0 +1,11 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<projectDescription>
+	<name>Theta_fourier</name>
+	<comment></comment>
+	<projects>
+	</projects>
+	<buildSpec>
+	</buildSpec>
+	<natures>
+	</natures>
+</projectDescription>

C/FFT_FP_realisation.c

@@ -0,0 +1,104 @@
+#include "FFT_FP_realisation.h"
+
+#include <math.h>
+#include <stdio.h>
+
+#ifndef FFT_FP_EXTERNAL_TWIDDLES
+static int64_t twiddle_re[TWIDDLE_L] = {0,};
+static int64_t twiddle_im[TWIDDLE_L] = {0,};
+#define PRINTF
+#endif
+
+
+
+void fft_twiddle_gen(int64_t* tw_re, int64_t* tw_im){
+	for (uint32_t k = 0; k < TWIDDLE_L; ++k){
+		double angle = 2.0 * PI * k / DATA_L;
+		tw_re[k] = lround(cos(angle) * FP_acc);
+		tw_im[k] = lround(-sin(angle) * FP_acc);
+	}
+}
+
+void fft_fp_prepare(void){
+	fft_twiddle_gen(twiddle_re, twiddle_im);
+	for (uint32_t k = 0; k < TWIDDLE_L; ++k){
+#ifdef PRINTF
+		printf("k, angle, tw_re, tw_im: %u %g %lld %lld\n", k, 2.0 * PI * k / DATA_L, (long long)twiddle_re[k], (long long)twiddle_im[k]);
+#endif
+	}
+}
+
+void FFT_fp(int64_t* inp, uint32_t inp_L, int64_t* buf){
+
+	    // buf имеет длину inp_L * 2 (Re, Im, Re, Im, ...)
+	    // inp содержит inp_L значений uint32_t
+
+	    uint32_t i, j, bit;
+//	    uint32_t N = inp_L / 2;  // длина комплексного массива (inp содержит только Re-часть)
+	    uint32_t N = inp_L;  // длина комплексного массива (inp содержит только Re-часть)
+
+	    // --- копирование входных данных в буфер ---
+	    for (i = 0; i < inp_L; i++) {
+	        buf[i * 2]     = inp[i];
+	        buf[i * 2 + 1] = 0;
+	    }
+
+	    // --- bit-reversal перестановка ---
+	    j = 0;
+	    for (i = 1; i < N; i++) {
+	        bit = N >> 1;
+	        while (j & bit) {
+	            j ^= bit;
+	            bit >>= 1;
+	        }
+	        j |= bit;
+	        if (i < j) {
+	            int64_t tmp_re = buf[i * 2];
+	            int64_t tmp_im = buf[i * 2 + 1];
+	            buf[i * 2]     = buf[j * 2];
+	            buf[i * 2 + 1] = buf[j * 2 + 1];
+	            buf[j * 2]     = tmp_re;
+	            buf[j * 2 + 1] = tmp_im;
+	        }
+	    }
+
+	    // --- уровни бабочек ---
+	    uint32_t m = 2;
+	    while (m <= N) {
+	        uint32_t half = m >> 1;
+	        uint32_t stride = N / m;
+
+	        for (uint32_t start = 0; start < N; start += m) {
+            for (uint32_t k = 0; k < half; k++) {
+                uint32_t tw_idx = k * stride;
+                int64_t wr = twiddle_re[tw_idx];
+                int64_t wi = twiddle_im[tw_idx];
+                int64_t u_re = buf[(start + k) * 2];
+                int64_t u_im = buf[(start + k) * 2 + 1];
+                int64_t v_re = buf[(start + k + half) * 2];
+                int64_t v_im = buf[(start + k + half) * 2 + 1];
+
+	                // t = w * v  (в фиксированной точке)
+	                int64_t t_re = (wr * v_re - wi * v_im) / FP_acc;
+	                int64_t t_im = (wr * v_im + wi * v_re) / FP_acc;
+
+	                // верх/низ
+	                buf[(start + k) * 2]              = u_re + t_re;
+	                buf[(start + k) * 2 + 1]          = u_im + t_im;
+	                buf[(start + k + half) * 2]       = u_re - t_re;
+	                buf[(start + k + half) * 2 + 1]   = u_im - t_im;
+
+            }
+        }
+        m <<= 1;
+    }
+	    for (uint32_t ii = 0; ii < inp_L; ++ii){
+	    	buf[ii*2] /= inp_L;
+	    	buf[ii*2+1] /= inp_L;
+#ifdef PRINTF
+	    	printf("re,im: %lld, %lld\n", (long long)buf[ii*2], (long long)buf[ii*2 +1]);
+#endif
+
+	    }
+
+}

C/FFT_FP_realisation.h

@@ -0,0 +1,31 @@
+#ifndef FFT_FP_REALISATION_H
+#define FFT_FP_REALISATION_H
+
+#include <stdint.h>
+
+#ifndef DATA_L
+#define DATA_L (1 << 16)
+#endif
+
+#ifndef FP_acc
+#define FP_acc 1000
+#endif
+
+#ifndef TWIDDLE_L
+#define TWIDDLE_L (DATA_L / 2)
+#endif
+
+#ifndef PI
+#define PI 3.141592653589793238462643383279502884197169399375105820974944
+#endif
+
+#ifdef FFT_FP_EXTERNAL_TWIDDLES
+extern volatile int64_t twiddle_re[TWIDDLE_L];
+extern volatile int64_t twiddle_im[TWIDDLE_L];
+#endif
+
+void fft_twiddle_gen(int64_t* twiddle_re, int64_t* twiddle_im);
+void fft_fp_prepare(void);
+void FFT_fp(int64_t* inp, uint32_t inp_L, int64_t* buf);
+
+#endif

C/FFT_realistaion.h

@@ -0,0 +1,24 @@
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+#define DATA_L (1<<17)
+#define FP_acc (1 << 10)
+#define TWIDDLE_L (DATA_L / 2)
+
+#define PI 3.141592653589793238462643383279502884197169399375105820974944
+
+int64_t twiddle_re[TWIDDLE_L] = {0,};
+int64_t twiddle_im[TWIDDLE_L] = {0,};
+
+void FFT_fp(int64_t* inp, uint32_t inp_L, int64_t* buf){
+
+}
+
+void twiddle_fill(int64_t* tw_re, int64_t* tw_im, uint32_t N){
+
+}
+
+

C/FFT_test.c

@@ -0,0 +1 @@
+ 

C/FP_math.c

@@ -0,0 +1,62 @@
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+#define DATA_L (1<<16)
+#define FP_acc (1<<16)
+
+#include "FFT_FP_realisation.h"
+
+
+void imagery2abs(int64_t* inp, uint32_t inp_L, int64_t *out_abs, int64_t *out_phase){
+	printf("+++++++++++++++++++++\n");
+	printf("calculating abs^2\n");
+	for (uint32_t i = 0; i < inp_L; ++i){
+		int64_t re = inp[i*2];
+		int64_t im = inp[i*2 +1];
+		int64_t val = (re * re + im * im) / FP_acc;
+		int64_t phase = 0;
+		printf("%lld, %ldd\n", (long long)val, (long long) phase);
+		out_abs[i] = val;
+		out_phase[i] = phase;
+
+	}
+}
+void data_generator(int64_t* X, int64_t* Y, uint32_t N){
+	for (int i = 0; i < N; ++i){
+		X[i] = (int64_t)i;
+		Y[i] = lround(sin(((double)i)*2.0* PI/(N/10))*FP_acc);
+		Y[i] += lround(cos(((double)i)*2.0* PI/(N/20))*FP_acc);
+		Y[i] += 1*FP_acc;
+	}
+}
+
+void main(){
+	char* logfilename = "tmp";
+	int64_t data_X[DATA_L] = {0,};
+	int64_t data_Y[DATA_L] = {0,};
+	int64_t data_res_FFT_imag[DATA_L*2] = {0,};
+	int64_t data_res_FFT_abs[DATA_L] = {0,};
+	int64_t data_res_FFT_phase[DATA_L] = {0,};
+
+
+	fft_fp_prepare();
+
+	data_generator(data_X, data_Y, DATA_L);
+
+
+	FFT_fp(data_Y, DATA_L, data_res_FFT_imag);
+
+
+	imagery2abs(data_res_FFT_imag, DATA_L, data_res_FFT_abs, data_res_FFT_phase);
+
+	FILE *logfile_ptr;
+    logfile_ptr = fopen(logfilename, "w");
+    fprintf(logfile_ptr, "X, Y, FFT_val, FFT_phase, FFT_Re, FFT_Im\n");
+	for (uint32_t i = 0; i < DATA_L; ++i){
+		fprintf(logfile_ptr, "%lld, %lg, %lg, %lg, %lg, %lg\n", (long long)data_X[i], ((double)data_Y[i])/(double)FP_acc , (double)data_res_FFT_abs[i]/(double)FP_acc, (double)data_res_FFT_phase[i]/(double)FP_acc,  (double)data_res_FFT_imag[2*i]/(double)FP_acc,  (double)data_res_FFT_imag[2*i +1]/(double)FP_acc     );
+	}
+	fclose(logfile_ptr);
+}

C/FP_math.c_backup

@@ -0,0 +1,152 @@
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+#define DATA_L 1024
+#define FP_acc 1000
+
+#define PI 3.141592653589793238462643383279502884197169399375105820974944
+
+
+
+int cos_delta_FP_arr[DATA_L] = {0,};
+int sin_delta_FP_arr[DATA_L] = {0,};
+
+void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){
+
+	    // buf имеет длину inp_L * 2 (Re, Im, Re, Im, ...)
+	    // inp содержит inp_L значений uint32_t
+
+	    uint32_t i, j, bit;
+//	    uint32_t N = inp_L / 2;  // длина комплексного массива (inp содержит только Re-часть)
+	    uint32_t N = inp_L;  // длина комплексного массива (inp содержит только Re-часть)
+
+	    // --- копирование входных данных в буфер ---
+	    for (i = 0; i < inp_L; i++) {
+	        buf[i * 2]     = (int)inp[i];  // Re
+	        buf[i * 2 + 1] = 0;            // Im = 0
+	    }
+
+	    // --- bit-reversal перестановка ---
+	    j = 0;
+	    for (i = 1; i < N; i++) {
+	        bit = N >> 1;
+	        while (j & bit) {
+	            j ^= bit;
+	            bit >>= 1;
+	        }
+	        j |= bit;
+	        if (i < j) {
+	            int tmp_re = buf[i * 2];
+	            int tmp_im = buf[i * 2 + 1];
+	            buf[i * 2]     = buf[j * 2];
+	            buf[i * 2 + 1] = buf[j * 2 + 1];
+	            buf[j * 2]     = tmp_re;
+	            buf[j * 2 + 1] = tmp_im;
+	        }
+	    }
+
+	    // --- уровни бабочек ---
+	    uint32_t m = 2;
+	    while (m <= N) {
+	        uint32_t half = m >> 1;
+
+	        int c_delta_FP = cos_delta_FP_arr[m - 2];
+	        int s_delta_FP = sin_delta_FP_arr[m - 2];
+
+	        for (uint32_t start = 0; start < N; start += m) {
+	            int wr = FP_acc;
+	            int wi = 0;
+
+	            for (uint32_t k = 0; k < half; k++) {
+	                int u_re = buf[(start + k) * 2];
+	                int u_im = buf[(start + k) * 2 + 1];
+	                int v_re = buf[(start + k + half) * 2];
+	                int v_im = buf[(start + k + half) * 2 + 1];
+
+	                // t = w * v  (в фиксированной точке)
+	                int64_t t_re = ((int64_t)wr * v_re - (int64_t)wi * v_im) / FP_acc;
+	                int64_t t_im = ((int64_t)wr * v_im + (int64_t)wi * v_re) / FP_acc;
+
+	                // верх/низ
+	                buf[(start + k) * 2]              = u_re + (int)t_re;
+	                buf[(start + k) * 2 + 1]          = u_im + (int)t_im;
+	                buf[(start + k + half) * 2]       = u_re - (int)t_re;
+	                buf[(start + k + half) * 2 + 1]   = u_im - (int)t_im;
+
+	                // w *= w_m (обновление twiddle)
+	                int64_t wr_new = ((int64_t)wr * c_delta_FP - (int64_t)wi * s_delta_FP) / FP_acc;
+	                int64_t wi_new = ((int64_t)wr * s_delta_FP + (int64_t)wi * c_delta_FP) / FP_acc;
+	                wr = (int)wr_new;
+	                wi = (int)wi_new;
+	            }
+	        }
+	        m <<= 1;
+	    }
+	    for (uint32_t i = 0; i < inp_L; ++i){
+	    	//buf[i*2] /= inp_L;
+	    	//buf[i*2+1] /= inp_L;
+	    	printf("re,im: %d, %d\n", buf[i*2], buf[i*2 +1]);
+	    }
+
+}
+
+void imagery2abs(uint32_t* inp, uint32_t inp_L, uint32_t *out){
+	printf("+++++++++++++++++++++\n");
+	printf("calculating abs^2\n");
+	for (uint32_t i = 0; i < inp_L; ++i){
+		int val = (inp[i*2]*inp[i*2] + inp[i*2 +1]*inp[i*2 +1])/FP_acc;
+		printf("%d\n", val);
+		out[i] = val;
+	}
+}
+
+
+
+void sin_cos_filler(int* sin_arr, int* cos_arr, uint32_t L){
+	//[int(-sin(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)]
+	for (uint32_t i = 0; i < L; ++i){
+		int m = i + 2;
+		double angle = 2.0* PI/m;
+		sin_arr[i] = lround(-sin(angle) * FP_acc);
+		cos_arr[i] = lround(cos(angle) * FP_acc);
+		printf("i, angle, sin, cos: %d %g %d %d\n", i, angle, sin_arr[i], cos_arr[i]);
+	}
+}
+
+void data_generator(uint32_t* X, uint32_t* Y, uint32_t N){
+	for (int i = 0; i < N; ++i){
+		X[i] = (int) i;
+		Y[i] = lround(sin(((double)i)*2.0* PI/(N/10))*FP_acc);
+//		Y[I] += 10*FP_acc;
+	}
+}
+
+void main(){
+	char* logfilename = "tmp";
+	uint32_t data_X[DATA_L] = {0,};
+	uint32_t data_Y[DATA_L] = {0,};
+	uint32_t data_res_FFT_imag[DATA_L*2] = {0,};
+	uint32_t data_res_FFT_abs[DATA_L] = {0,};
+
+	sin_cos_filler(sin_delta_FP_arr, cos_delta_FP_arr, DATA_L);
+
+	data_generator(data_X, data_Y, DATA_L);
+
+
+	FFT_fp(data_Y, DATA_L, data_res_FFT_imag);
+
+
+	imagery2abs(data_res_FFT_imag, DATA_L, data_res_FFT_abs);
+
+	FILE *logfile_ptr;
+    logfile_ptr = fopen(logfilename, "w");
+    fprintf(logfile_ptr, "X, Y, FFT\n");
+	for (uint32_t i = 0; i < DATA_L; ++i){
+		fprintf(logfile_ptr, "%d, %d, %d \n", data_X[i], data_Y[i],data_res_FFT_abs[i]);
+	}
+	fclose(logfile_ptr);
+}
+

C/Makefile

@@ -0,0 +1,14 @@
+SRCS := FP_math.c FFT_FP_realisation.c
+OBJS := $(SRCS:.c=.o)
+TARGET := FP_math
+
+all: $(TARGET)
+
+$(TARGET): $(OBJS)
+	gcc $(OBJS) -o $@ -lm
+
+%.o: %.c
+	gcc -c $< -o $@
+
+clean:
+	-rm -f $(TARGET) $(OBJS)

C/plotter.py

@@ -0,0 +1,57 @@
+#!/usr/bin/pypy3
+import plotly.graph_objs as go
+from decimal import *
+from sys import argv
+import csv
+
+
+def main():
+  chart = go.Figure()
+  chart.add_trace(go.Scatter(x=[1,2,3], y=[1,2,3]))
+  chart.show()
+
+
+if __name__ == "__main__":
+  if len(argv) == 1:
+    main()
+  else:
+    f = open(argv[1], "rt")
+    #csv_reader = csv.reader(f)
+    data = {}
+#    for line in csv_reader:
+    headers = f.readline()
+    headers = headers.split(",")
+    headers[-1] = headers[-1][:-1]
+    for i in range(len(headers)):
+      headers[i] = headers[i].strip()
+    for header in headers:
+      data[header] = []
+    for line in f:
+      try:
+        line_data = line.split(",")
+        #line_data = line
+        #print(line_data)
+        for i in range(len(line_data)):
+          val = float(line_data[i])
+          if (i < len(headers)):
+            header = headers[i]
+            data[header].append(val)
+          else:
+            print("too much values!")
+            print("headers:", headers)
+            print("values:", line_data)
+      except ValueError:
+        pass
+    f.close()
+#    print(data)
+    chart = go.Figure()
+    keys = [k for k in data.keys()]
+#    print("data.keys:", keys)
+    X_name = keys[0]
+    print("data samples:",len(data[X_name]))
+    for key, val in data.items():
+      if key != X_name:
+#        chart.add_trace(go.Scatter(x=data[X_name], y=val, name=key, mode="markers+lines"))
+        chart.add_trace(go.Scatter(x=data[X_name], y=val, name=key, mode="lines"))
+    chart.show()
+

FFT.py

@@ -1,32 +0,0 @@
-import numpy as np
-
-
-
-
-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):
-  inp = np.array(inp)
-  out = FFT(inp)
-  out = [val for val in np.abs(out)]
-  print(out)
-  return out
-
-
-def FFT(inp):
-  inp = np.array(inp)
-  N = inp.shape[0]
-  if N == 1:
-    return inp
-  X_even = FFT(inp[::2])
-  X_odd  = FFT(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))
-

python/FFT.py

@@ -0,0 +1,326 @@
+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 = 2<<16
+
+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]
+
+  # --- уровни бабочек ---
+  
+  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
+    
+  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
+    
+
+
+
+
+
+
+
+
+
+
+

python/FP_trigonometry.py

@@ -3,42 +3,18 @@ from math import sin, cos, pi, sqrt
 import plotly.graph_objs as go
 from plotly.subplots import make_subplots
 
-from FFT import FFT_real
+FP_acc = 2<<16
-
-FP_acc = 1e3
-
-INP_L = 1024
-
-
-F_nyquist = INP_L//2
-
-
 pi_FP = 1* FP_acc
 
 
-def abs_f(re, im):
+
-  return sqrt(re*re + im*im)
-  
-  
 def abs_FP(re, im):
   return int(sqrt(re*re + im*im))
   
 def sqrt_FP(val):
   #print(val)
-  return int(sqrt(val))
+  return int(sqrt(val)) 
 
-def DFT_naive(inp, out):
-  for f in range(len(out)):
-    val_re = 0
-    val_im = 0
-    for n in range(len(inp)):
-      phi = 2*pi*f*n/INP_L
-      val_re += inp[n] * sin(phi) /INP_L
-      val_im += inp[n] * cos(phi) /INP_L
-      
-    val_abs = abs_f(val_re, val_im)
-    print("F, val_abs:",f, val_abs)
-    out[f] = val_abs
 
 
 def abs_FP(re, im):
@@ -110,7 +86,7 @@ def sin_FP_constrained(phi_fp):
       print("phi_fp:", phi_fp," >",phi_trh,"...   recusion")
       print("phi_fp/2:", phi_fp/2)
     sin_phi_05 = sin_FP_constrained(phi_fp/2)
-    sin_phi = (2*sin_phi_05 * sqrt_FP(FP_acc**2 - sin_phi_05*sin_phi_05))/FP_acc
+    sin_phi = int((2*sin_phi_05 * sqrt_FP(FP_acc**2 - sin_phi_05*sin_phi_05))/FP_acc)
     if (sin_05_debug):
       print("sin_phi_05:", sin_phi_05)
       print("sin_phi:",sin_phi)
@@ -132,80 +108,8 @@ def sin_FP_constrained(phi_fp):
   
 def cos_FP(phi_fp):
   return sin_FP(phi_fp - pi_FP/2)
-
-
-def DFT_naive_FP(inp_float, out):
-  inp = [val*FP_acc for val in inp_float]
-  for f in range(len(out)):
-    val_re = 0
-    val_im = 0
-    for n in range(len(inp)):
-      phi = 2*pi_FP*f*n/INP_L
-      phi_sin = sin_FP(phi)
-      phi_cos = cos_FP(phi)
-      #print(phi, phi_sin, phi_cos)
-      val_re += inp[n] * phi_sin /INP_L
-      val_im += inp[n] * phi_cos /INP_L
-      
-    val_abs = abs_FP(val_re, val_im)
-    print("F, val_abs:",f, val_abs)
-    out[f] = val_abs
   
-
-
-
   
-def FFT_naive(inp, out):
-  fft_out = FFT_real(inp)
-  for i in range(len(fft_out)):
-    val = fft_out[i]/len(inp)
-    if (i < len(out)):
-      out[i] = val
-    else:
-      pass
-      #out.append(val)
-
-def FFT_tester():
-  inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
-#  inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
-  out_DFT = [0 for i in range(F_nyquist + 1)]
-  out_FFT = [0 for val in range(F_nyquist + 1)]
-  DFT_naive(inp, out_DFT)
-  FFT_naive(inp, out_FFT)
-  Fourier_error = []
-  for a,b in zip(out_FFT, out_DFT):
-    Fourier_error.append(a - b)
-  chart = make_subplots(rows=3, cols=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(out_DFT))], y=out_DFT, name="out_DFT", mode="markers+lines"), row=2, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FFT))], y=out_FFT, name="out_FFT", mode="markers+lines"), row=2, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(Fourier_error))], y=Fourier_error, name="error", mode="markers+lines"), row=3, col=1)
-  chart.show()
-
-
-
-
-
-def DFT_tester():
-  inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
-#  inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
-  out_float = [0 for i in range(F_nyquist + 1)]
-  out_FP = [0 for val in out_float]
-  DFT_naive(inp, out_float)
-  DFT_naive_FP(inp, out_FP)
-  FP_error = []
-  for a,b in zip(out_float, out_FP):
-    FP_error.append(a - b/FP_acc)
-  chart = make_subplots(rows=3, cols=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(out_float))], y=out_float, name="out_float", mode="markers+lines"), row=2, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=[val/FP_acc for val in out_FP], name="out_FP", mode="markers+lines"), row=2, col=1)
-  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=FP_error, name="FP_error", mode="markers+lines"), row=3, col=1)
-  chart.show()
-  
-
-
-
 def sin_tester():
   N = 4000
   angs = [(i - N/2)/1000 for i in range(N)]
@@ -220,19 +124,26 @@ def sin_tester():
     print("angle, sin, cos:",phi, val_fp, val_cos_fp)
     res_FP.append(val_fp)
     res_cos_FP.append(val_cos_fp)
-  chart = go.Figure()
+  chart = make_subplots(rows=2, cols=1)
-  chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_float", mode="markers+lines"))
+  
-  chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_FP], name="sin_FP", mode="markers+lines"))
+  chart.add_trace(go.Scatter(x = angs, y=res_f, name="sin_float", mode="markers+lines"), row=1, col=1)
-  chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_cos_FP], name="cos_FP", mode="markers+lines"))
+  chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_FP], name="sin_FP", mode="markers+lines"), row=1, col=1)
+  chart.add_trace(go.Scatter(x = angs, y=[val/FP_acc for val in res_cos_FP], name="cos_FP", mode="markers+lines"), row=1, col=1)
+  chart.add_trace(go.Scatter(x = angs, y=[a - b/FP_acc for a,b in zip(res_f, res_FP)], name="error", mode="markers+lines"), row=2, col=1)
+  chart.update_xaxes(matches="x1", row=2, col=1)
+
   chart.show()
   
-  
-  
-  
 
 
 if __name__ == "__main__":
-  #main()
+  sin_tester()
-#  DFT_tester()
+  
-  FFT_tester()
+  
-  #sin_tester()
+  
+  
+  
+  
+  
+  
+  

python/naive_DFT.py

@@ -0,0 +1,124 @@
+#!/usr/bin/python3
+from math import sin, cos, pi, sqrt
+import plotly.graph_objs as go
+from plotly.subplots import make_subplots
+
+from FFT import FFT_real
+from FP_trigonometry import *
+
+
+INP_L = 1024
+
+
+F_nyquist = INP_L//2
+
+
+
+
+
+def abs_f(re, im):
+  return sqrt(re*re + im*im)
+  
+  
+
+
+def DFT_naive(inp, out):
+  for f in range(len(out)):
+    val_re = 0
+    val_im = 0
+    for n in range(len(inp)):
+      phi = 2*pi*f*n/INP_L
+      val_re += inp[n] * sin(phi) /INP_L
+      val_im += inp[n] * cos(phi) /INP_L
+      
+    val_abs = abs_f(val_re, val_im)
+    #print("F, val_abs:",f, val_abs)
+    out[f] = val_abs
+
+
+
+
+def DFT_naive_FP(inp_float, out):
+  inp = [val*FP_acc for val in inp_float]
+  for f in range(len(out)):
+    val_re = 0
+    val_im = 0
+    for n in range(len(inp)):
+      phi = 2*pi_FP*f*n/INP_L
+      phi_sin = sin_FP(phi)
+      phi_cos = cos_FP(phi)
+      #print(phi, phi_sin, phi_cos)
+      val_re += inp[n] * phi_sin /INP_L
+      val_im += inp[n] * phi_cos /INP_L
+      
+    val_abs = abs_FP(val_re, val_im)
+    #print("F, val_abs:",f, val_abs)
+    out[f] = val_abs
+  
+
+
+
+  
+def FFT_naive(inp, out):
+  fft_out = FFT_real(inp)
+  for i in range(len(fft_out)):
+    val = fft_out[i]/len(inp)
+    if (i < len(out)):
+      out[i] = val
+    else:
+      pass
+      #out.append(val)
+
+def FFT_tester():
+  inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
+#  inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
+  out_DFT = [0 for i in range(F_nyquist + 1)]
+  out_FFT = [0 for val in range(F_nyquist + 1)]
+  DFT_naive(inp, out_DFT)
+  FFT_naive(inp, out_FFT)
+  Fourier_error = []
+  for a,b in zip(out_FFT, out_DFT):
+    Fourier_error.append(a - b)
+  chart = make_subplots(rows=3, cols=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(out_DFT))], y=out_DFT, name="out_DFT", mode="markers+lines"), row=2, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FFT))], y=out_FFT, name="out_FFT", mode="markers+lines"), row=2, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(Fourier_error))], y=Fourier_error, name="error", mode="markers+lines"), row=3, col=1)
+  chart.update_xaxes(matches="x2", row=3, col=1)
+  chart.show()
+
+
+
+
+
+def DFT_tester():
+  inp = [-1 + 0.01*i + sin(2*pi*i/10) + cos(2*pi*i/20) + sin(2*pi*i/250) + sin(2*pi*i/2.001) for i in range(INP_L)]
+#  inp = [sin(2*pi*i/2.001)for i in range(INP_L)]
+  out_float = [0 for i in range(F_nyquist + 1)]
+  out_FP = [0 for val in out_float]
+  DFT_naive(inp, out_float)
+  DFT_naive_FP(inp, out_FP)
+  FP_error = []
+  for a,b in zip(out_float, out_FP):
+    FP_error.append(a - b/FP_acc)
+  chart = make_subplots(rows=3, cols=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(inp))], y=inp, name="inp", mode="markers+lines"), row=1, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(out_float))], y=out_float, name="out_float", mode="markers+lines"), row=2, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=[val/FP_acc for val in out_FP], name="out_FP", mode="markers+lines"), row=2, col=1)
+  chart.add_trace(go.Scatter(x=[i for i in range(len(out_FP))], y=FP_error, name="FP_error", mode="markers+lines"), row=3, col=1)
+  chart.update_xaxes(matches="x2", row=3, col=1)
+  chart.show()
+  
+
+
+
+
+  
+  
+
+
+if __name__ == "__main__":
+  #main()
+#  DFT_tester()
+  FFT_tester()
+  #sin_tester()