fft-debug.c

// 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;
}