ChTheo/FFT_and_FP_math
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added precalculated twiddles table, fixed accuracy issues (switched from int32 to int64). Now FFT working correctly

Browse Source
This commit is contained in:
Theodor Chikin
2025-10-08 22:33:32 +03:00
parent 1a9f820724
commit f54daffb72
3 changed files with 1074 additions and 1071 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
C/FP_math
View File

Binary file not shown.

97
C/FP_math.c
View File

@ -11,10 +11,12 @@
int cos_delta_FP_arr[DATA_L] = {0,};
int sin_delta_FP_arr[DATA_L] = {0,};
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(uint32_t* inp, uint32_t inp_L, uint32_t* buf){
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
@ -25,8 +27,8 @@ void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){
// --- копирование входных данных в буфер ---
for (i = 0; i < inp_L; i++) {
buf[i * 2] = (int)inp[i]; // Re
buf[i * 2 + 1] = 0; // Im = 0
buf[i * 2] = inp[i];
buf[i * 2 + 1] = 0;
}
// --- bit-reversal перестановка ---
@ -39,8 +41,8 @@ void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){
}
j |= bit;
if (i < j) {
int tmp_re = buf[i * 2];
int tmp_im = buf[i * 2 + 1];
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;
@ -53,83 +55,85 @@ void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){
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;
int64_t *twiddle_re = wr_arr[m - 2];
int64_t *twiddle_im = wi_arr[m - 2];
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];
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 = ((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;
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 + (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;
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;
// 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;
}
}
}
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]);
printf("re,im: %lld, %lld\n", (long long)buf[i*2], (long long)buf[i*2 +1]);
}
}
void imagery2abs(uint32_t* inp, uint32_t inp_L, uint32_t *out){
void imagery2abs(int64_t* inp, uint32_t inp_L, int64_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);
int64_t re = inp[i*2];
int64_t im = inp[i*2 +1];
int64_t val = (re * re + im * im) / FP_acc;
printf("%lld\n", (long long)val);
out[i] = val;
}
}
void sin_cos_filler(int* sin_arr, int* cos_arr, uint32_t L){
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);
printf("i, angle, sin, cos: %d %g %d %d\n", i, angle, sin_arr[i], cos_arr[i]);
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(uint32_t* X, uint32_t* Y, uint32_t N){
void data_generator(int64_t* X, int64_t* Y, uint32_t N){
for (int i = 0; i < N; ++i){
X[i] = (int) i;
X[i] = (int64_t)i;
Y[i] = lround(sin(((double)i)*2.0* PI/(N/10))*FP_acc);
// Y[I] += 10*FP_acc;
Y[i] += 1*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,};
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,};
sin_cos_filler(sin_delta_FP_arr, cos_delta_FP_arr, DATA_L);
@ -145,8 +149,7 @@ void main(){
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]);
fprintf(logfile_ptr, "%lld, %lg, %lg \n", (long long)data_X[i], ((double)data_Y[i])/(double)FP_acc , (double)data_res_FFT_abs[i]/(double)FP_acc);
}
fclose(logfile_ptr);
}

2048
C/tmp
View File

File diff suppressed because it is too large Load Diff