#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <tgmath.h>
#include <fftw3.h>

#if defined(USE_LONG_DOUBLE)
	typedef long double         real_t;
	typedef long double complex complex_t;
	typedef fftwl_complex       fftw_complex_t;
	typedef fftwl_plan          fftw_plan_t;
	#define MACHINE_EPSILON LDBL_EPSILON
	#define REAL_C(x)       x ## l
	#define PRIr            "Lf"
	#define SCNr            "Lf"
	#define FFTW_MALLOC(size)                        fftwl_malloc(size)
	#define FFTW_PLAN_DFT_1D(n, in, out, sign, flag) fftwl_plan_dft_1d(n, in, out, sign, flag)
	#define FFTW_EXECUTE(plan)                       fftwl_execute(plan)
	#define FFTW_DESTROY_PLAN(plan)                  fftwl_destroy_plan(plan)
	#define FFTW_FREE(p)                             fftwl_free(p)
#elif defined(USE_DOUBLE)
	typedef double         real_t;
	typedef double complex complex_t;
	typedef fftw_complex   fftw_complex_t;
	typedef fftw_plan      fftw_plan_t;
	#define MACHINE_EPSILON DBL_EPSILON
	#define REAL_C(x)       x
	#define PRIr            "f"
	#define SCNr            "lf"
	#define FFTW_MALLOC(size)                        fftw_malloc(size)
	#define FFTW_PLAN_DFT_1D(n, in, out, sign, flag) fftw_plan_dft_1d(n, in, out, sign, flag)
	#define FFTW_EXECUTE(plan)                       fftw_execute(plan)
	#define FFTW_DESTROY_PLAN(plan)                  fftw_destroy_plan(plan)
	#define FFTW_FREE(p)                             fftw_free(p)
#elif defined(USE_FLOAT)
	typedef float         real_t;
	typedef float complex complex_t;
	typedef fftwf_complex fftw_complex_t;
	typedef fftwf_plan    fftw_plan_t;
	#define MACHINE_EPSILON FLT_EPSILON
	#define REAL_C(x)       x ## f
	#define PRIr            "f"
	#define SCNr            "f"
	#define FFTW_MALLOC(size)                        fftwf_malloc(size)
	#define FFTW_PLAN_DFT_1D(n, in, out, sign, flag) fftwf_plan_dft_1d(n, in, out, sign, flag)
	#define FFTW_EXECUTE(plan)                       fftwf_execute(plan)
	#define FFTW_DESTROY_PLAN(plan)                  fftwf_destroy_plan(plan)
	#define FFTW_FREE(p)                             fftwf_free(p)
#else
	#error "Specify USE_* for floating point precision!"
#endif

#define N_FFT 64
#define M_SP  10
#define N_SP  16
#define M_LP   2
#define N_LP  64
#define N_GI  16
#define N_PRE (M_SP * N_SP + M_LP * N_GI + M_LP * N_LP)

int main(int argc, char *argv[]) {
	fftw_complex_t *pre;
	fftw_complex_t *in, *out;
	fftw_plan_t plan;

	pre  = (fftw_complex_t *)FFTW_MALLOC(sizeof(fftw_complex_t) * N_PRE);
	in   = (fftw_complex_t *)FFTW_MALLOC(sizeof(fftw_complex_t) * N_FFT);
	out  = (fftw_complex_t *)FFTW_MALLOC(sizeof(fftw_complex_t) * N_FFT);
	plan = FFTW_PLAN_DFT_1D(N_FFT, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);

	// SP
	{
		for (int f = -N_FFT/2; f < N_FFT/2; f++) {
			// NOTE: FFTW maps positive frequencies [0, N_FFT/2) to array index [0, N_FFT/2),
			//       and negative frequencies [-N_FFT/2, 0) to array index [N_FFT/2, N_FFT).
			int i =
				(f < 0) ? f + N_FFT:
				f;
			in[i] =
				(f == -24 ||
				 f == -16 ||
				 f ==  -4 ||
				 f ==  12 ||
				 f ==  16 ||
				 f ==  20 ||
				 f ==  24) ?  sqrt(REAL_C(13.0) / REAL_C(6.0)) * (REAL_C(1.0) + REAL_C(1.0) * I) :
				(f == -20 ||
				 f == -12 ||
				 f ==  -8 ||
				 f ==   4 ||
				 f ==   8) ? -sqrt(REAL_C(13.0) / REAL_C(6.0)) * (REAL_C(1.0) + REAL_C(1.0) * I) :
				REAL_C(0.0) + REAL_C(0.0) * I;
		}

		FFTW_EXECUTE(plan);

		// NOTE: FFTW calculate backward FFT w/o 1/N_FFT scaling.
		for (int i = 0; i < N_FFT; i++) out[i] /= (real_t)N_FFT;

		for (int m = 0; m < M_SP; m++) {
			for (int n = 0; n < N_SP; n++) {
				// NOTE: FFTW maps positive samples [0, N_FFT/2) to array index [0, N_FFT/2),
				//       and negative samples [-N_FFT/2, 0) to array index [N_FFT/2, N_FFT).
				int i = n;
				pre[m * N_SP + n] = out[i];
			}
		}
	}

	// LP
	{
		for (int f = -N_FFT/2; f < N_FFT/2; f++) {
			// NOTE: FFTW maps positive frequencies [0, N_FFT/2) to array index [0, N_FFT/2),
			//       and negative frequencies [-N_FFT/2, 0) to array index [N_FFT/2, N_FFT).
			int i =
				(f < 0) ? f + N_FFT :
				f;
			in[i] =
				(f == -26 ||
				 f == -25 ||
				 f == -22 ||
				 f == -21 ||
				 f == -19 ||
				 f == -17 ||
				 f == -16 ||
				 f == -15 ||
				 f == -14 ||
				 f == -13 ||
				 f == -12 ||
				 f ==  -9 ||
				 f ==  -8 ||
				 f ==  -6 ||
				 f ==  -4 ||
				 f ==  -3 ||
				 f ==  -2 ||
				 f ==  -1 ||
				 f ==   1 ||
				 f ==   4 ||
				 f ==   5 ||
				 f ==   7 ||
				 f ==   9 ||
				 f ==  15 ||
				 f ==  16 ||
				 f ==  19 ||
				 f ==  21 ||
				 f ==  23 ||
				 f ==  24 ||
				 f ==  25 ||
				 f ==  26) ?  REAL_C(1.0) + REAL_C(0.0) * I :
				(f == -24 ||
				 f == -23 ||
				 f == -20 ||
				 f == -18 ||
				 f == -11 ||
				 f == -10 ||
				 f ==  -7 ||
				 f ==  -5 ||
				 f ==   2 ||
				 f ==   3 ||
				 f ==   6 ||
				 f ==   8 ||
				 f ==  10 ||
				 f ==  11 ||
				 f ==  12 ||
				 f ==  13 ||
				 f ==  14 ||
				 f ==  17 ||
				 f ==  18 ||
				 f ==  20 ||
				 f ==  22) ? REAL_C(-1.0) + REAL_C(0.0) * I :
				REAL_C(0.0) + REAL_C(0.0) * I;
		}

		FFTW_EXECUTE(plan);

		// NOTE: FFTW calculate backward FFT w/o 1/N_FFT scaling.
		for (int i = 0; i < N_FFT; i++) out[i] /= (real_t)N_FFT;

		for (int m = 0; m < M_LP; m++) {
			for (int n = 0; n < N_LP; n++) {
				// NOTE: FFTW maps positive samples [0, N_FFT/2) to array index [0, N_FFT/2),
				//       and negative samples [-N_FFT/2, 0) to array index [N_FFT/2, N_FFT).
				int i = n;
				pre[M_SP * N_SP + M_LP * N_GI + m * N_LP + n] = out[i];
			}
		}
	}

	// 2GI
	{
		for (int m = 0; m < M_LP; m++) {
			for (int n = 0; n < N_GI; n++) {
				pre[M_SP * N_SP + (m * N_GI + n)] = pre[N_PRE - M_LP * N_GI + m * N_GI + n];
			}
		}
	}

#if defined(USE_SWS)
	// symbol wave shaping
	// FIXME: This should be performed in reverse order, ATM...
	pre[M_SP * N_SP] += pre[0];
	pre[M_SP * N_SP] /= REAL_C(2.0);
	pre[0] /= REAL_C(2.0);
#endif

	// dump
	printf("# PRE \n");
	for (int n = 0; n < N_PRE; n++) {
		printf("%d %" PRIr " %" PRIr "\n",
			n,
			creal(pre[n]),
			cimag(pre[n])
		);
	}

	FFTW_DESTROY_PLAN(plan);
	FFTW_FREE(out);
	FFTW_FREE(in);
	FFTW_FREE(pre);

	return 0;
}