replaced huge 2d twiddle arrays by smaller 1d arrays. Now it uses very low memory

C/FFT_FP_realisation.c

@@ -0,0 +1,88 @@
+#include "FFT_FP_realisation.h"
+
+#include <math.h>
+#include <stdio.h>
+
+static int64_t twiddle_re[TWIDDLE_L] = {0,};
+static int64_t twiddle_im[TWIDDLE_L] = {0,};
+
+void fft_fp_prepare(void){
+	for (uint32_t k = 0; k < TWIDDLE_L; ++k){
+		double angle = 2.0 * PI * k / DATA_L;
+		twiddle_re[k] = lround(cos(angle) * FP_acc);
+		twiddle_im[k] = lround(-sin(angle) * FP_acc);
+		printf("k, angle, tw_re, tw_im: %u %g %lld %lld\n", k, angle, (long long)twiddle_re[k], (long long)twiddle_im[k]);
+	}
+}
+
+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;
+	    	printf("re,im: %lld, %lld\n", (long long)buf[ii*2], (long long)buf[ii*2 +1]);
+	    }
+
+}

C/FFT_FP_realisation.h

@@ -0,0 +1,25 @@
+#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
+
+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

@@ -4,91 +4,12 @@
 #include <stdint.h>
 #include <math.h>
 
-#define DATA_L 1024
+#define DATA_L (1<<16)
 #define FP_acc 1000
 
-#define PI 3.141592653589793238462643383279502884197169399375105820974944
+#include "FFT_FP_realisation.h"
 
 
-
-int64_t cos_delta_FP_arr[DATA_L] = {0,};
-int64_t sin_delta_FP_arr[DATA_L] = {0,};
-int64_t wr_arr[DATA_L][DATA_L / 2] = {0,};
-int64_t wi_arr[DATA_L][DATA_L / 2] = {0,};
-
-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;
-
-	        for (uint32_t start = 0; start < N; start += m) {
-            int64_t *twiddle_re = wr_arr[m - 2];
-            int64_t *twiddle_im = wi_arr[m - 2];
-
-            for (uint32_t k = 0; k < half; k++) {
-                int64_t wr = twiddle_re[k];
-                int64_t wi = twiddle_im[k];
-                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 i = 0; i < inp_L; ++i){
-	    	//buf[i*2] /= inp_L;
-	    	//buf[i*2+1] /= inp_L;
-	    	printf("re,im: %lld, %lld\n", (long long)buf[i*2], (long long)buf[i*2 +1]);
-	    }
-
-}
-
 void imagery2abs(int64_t* inp, uint32_t inp_L, int64_t *out){
 	printf("+++++++++++++++++++++\n");
 	printf("calculating abs^2\n");
@@ -100,26 +21,6 @@ void imagery2abs(int64_t* inp, uint32_t inp_L, int64_t *out){
 		out[i] = val;
 	}
 }
-
-
-
-void sin_cos_filler(int64_t* sin_arr, int64_t* 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);
-		uint32_t half = m >> 1;
-		for (uint32_t k = 0; k < half; ++k){
-			double tw_angle = angle * k;
-			wr_arr[i][k] = lround(cos(tw_angle) * FP_acc);
-			wi_arr[i][k] = lround(-sin(tw_angle) * FP_acc);
-		}
-		printf("i, angle, sin, cos: %d %g %lld %lld\n", i, angle, (long long)sin_arr[i], (long long)cos_arr[i]);
-	}
-}
-
 void data_generator(int64_t* X, int64_t* Y, uint32_t N){
 	for (int i = 0; i < N; ++i){
 		X[i] = (int64_t)i;
@@ -135,7 +36,7 @@ void main(){
 	int64_t data_res_FFT_imag[DATA_L*2] = {0,};
 	int64_t data_res_FFT_abs[DATA_L] = {0,};
 
-	sin_cos_filler(sin_delta_FP_arr, cos_delta_FP_arr, DATA_L);
+	fft_fp_prepare();
 
 	data_generator(data_X, data_Y, DATA_L);
 

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

@@ -1,5 +1,14 @@
-all:
-	gcc FP_math.c -o FP_math -lm
+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 FP_math
+	-rm -f $(TARGET) $(OBJS)

C/plotter.py

@@ -51,6 +51,7 @@ if __name__ == "__main__":
     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="markers+lines"))
+        chart.add_trace(go.Scatter(x=data[X_name], y=val, name=key, mode="lines"))
     chart.show()