Skip to content

Instantly share code, notes, and snippets.

@x42
Created August 29, 2021 19:42
Show Gist options
  • Save x42/17fd61b04d4c4939727dfdccd79f53a5 to your computer and use it in GitHub Desktop.
Save x42/17fd61b04d4c4939727dfdccd79f53a5 to your computer and use it in GitHub Desktop.
// gcc -o /tmp/fft-debug fft-debug.c `pkg-config --libs fftw3f` -lm
// /tmp/fft-debug > /tmp/fft.dat
// gnuplot
// plot '/tmp/fft.dat' u 3 w l, '' u 4 w l, '' u 5 w lp
#ifndef _GNU_SOURCE
#define _GNU_SOURCE /* needed for M_PI */
#endif
#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#define BLKSIZE 2048
#define FFTLEN (2 * BLKSIZE)
#define INSIZE (32 * BLKSIZE)
int
main ()
{
double angle = 90 / 360.0;
const float ca = cos (2 * M_PI * angle);
const float sa = sin (2 * M_PI * angle);
/* input data and expected output */
float input[INSIZE];
float expect[INSIZE];
for (int i = 0; i < INSIZE; i++) {
double twopi = 2.0 * M_PI;
#if 1
double phase = 2.7 * i / (float)FFTLEN;
#else
double phase = 0.37 * i / (float)FFTLEN;
#endif
input[i] = .5 * cosf (phase * twopi);
expect[i] = .5 * cosf (phase * twopi + angle * twopi);
#if 0 /* add DC offset */
input[i] += .1;
expect[i] += .1;
#endif
}
/* normalized Hann window */
float window[FFTLEN];
const double norm = 0.5 / FFTLEN;
for (int i = 0; i < FFTLEN; i++) {
window[i] = norm * (1 - cos (2.0 * M_PI * i / (float)FFTLEN));
}
/* FFT buffers and plan */
float time_data[FFTLEN];
fftwf_complex freq_data[(FFTLEN / 2) + 1];
fftwf_plan plan_r2c = fftwf_plan_dft_r2c_1d (FFTLEN, time_data, freq_data, FFTW_ESTIMATE);
fftwf_plan plan_c2r = fftwf_plan_dft_c2r_1d (FFTLEN, freq_data, time_data, FFTW_ESTIMATE);
/* ***************************************************************************
* iterative process, overlap 1/2 FFTLEN (== BLKSIZE)
*/
int offset = 0;
int remain = INSIZE - FFTLEN;
while (remain > 0) {
/* copy end/overlap of previous iFFT */
float buf_out[BLKSIZE];
memcpy (buf_out, time_data + BLKSIZE, sizeof (float) * BLKSIZE);
memcpy (time_data, &input[offset], FFTLEN * sizeof (float));
#if 0 /* apply window *before* processing */
for (int i = 0; i < FFTLEN; i++) {
time_data[i] *= window[i];
}
#endif
/* forward transform*/
fftwf_execute_dft_r2c (plan_r2c, time_data, freq_data);
/* rotate phase */
for (int i = 1 /*1: skip DC, 0Hz */; i <= FFTLEN / 2; ++i) {
const float re = freq_data[i][0];
const float im = freq_data[i][1];
freq_data[i][0] = (ca * re - sa * im);
freq_data[i][1] = (sa * re + ca * im);
}
/* backward transform*/
fftwf_execute_dft_c2r (plan_c2r, freq_data, time_data);
#if 1 /* apply window *after* processing */
for (int i = 0; i < FFTLEN; i++) {
time_data[i] *= window[i];
}
#endif
for (int i = 0; i < BLKSIZE; ++i) {
buf_out[i] += time_data[i];
}
for (int i = 0; i < BLKSIZE; ++i) {
printf ("%d %f %f %f %f\n",
i + offset,
input[i + offset],
expect[i + offset],
buf_out[i],
(i + offset) < BLKSIZE ? 0 : (buf_out[i] - expect[i + offset]));
}
offset += BLKSIZE;
remain -= BLKSIZE;
}
fftwf_destroy_plan (plan_r2c);
fftwf_destroy_plan (plan_c2r);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment