97 lines
3.2 KiB
C
97 lines
3.2 KiB
C
#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,};
|
|
#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){
|
|
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]);
|
|
}
|
|
}
|
|
|
|
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]);
|
|
}
|
|
|
|
}
|