#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fft.h"

#define N_SAMPLES (1 << 14)

typedef enum
{
    MITCHELL,
    LANCZOS1,
    LANCZOS2,
    LANCZOS2_STRETCHED,
    LANCZOS25,
    LANCZOS3,
    LANCZOS3_SHRINK,
    LANCZOS4,
    LANCZOS64,
    BOX,
    BOX15,
    TANH,
    TANH_SHRINK,
    QUADRATIC,
    FAST_CUBIC,
    LINEAR,
    GAMMA,
    NONGAMMA,
} filter_t;

static filter_t filter = GAMMA;

static double
MIN (double a, double b)
{
    return a < b ? a : b;
}

int
main ()
{
    complex_t input[N_SAMPLES];
    complex_t test[N_SAMPLES];
    double first = -64;
    double last = 64;
    int i;
    double sample_extent = (last - first) / N_SAMPLES;
    double n_samples = 1 / sample_extent;
    double sum = 0.0;

    for (i = 0; i < N_SAMPLES; ++i)
    {
	double sn = (i + N_SAMPLES/2) % N_SAMPLES - N_SAMPLES/2;
	double t = (sn * sample_extent);
	double iv = 0.0;
	double v;

	t = t * 1.0;

	switch (filter)
	{
	case MITCHELL:
	{
	    const double B = 0.5;
	    const double C = 0.6;
	    
	    if (fabs (t) < 1)
	    {
		v = (12 - 9 * B - 6 * C) * fabs (t) * fabs (t) * fabs (t)	+
		    (-18 + 12 * B + 6 * C) * fabs (t) * fabs (t)		+
		    (6 - 2 * B);
	    }
	    else if (fabs (t) < 2)
	    {
		v = (-B - 6 * C) * fabs (t) * fabs(t) * fabs(t)			+
		    (6 * B + 30 * C) * fabs(t) * fabs(t)			+
		    (-12 * B - 48 * C) * fabs(t)				+
		    (8 * B + 24 * C);
	    }
	    else
	    {
		v = 0.0;
	    }

	    v = v / 6.0;
	}
	break;

	case LANCZOS2:
	{
	    if (fabs(t) == 0.0)
	    {
		v = 1;
	    }
	    else if (fabs (t) > 2)
		v = 0.0;
	    else
		v = sin (M_PI * t) / (M_PI * t) * sin (M_PI * t / 2.0) / (M_PI * t / 2.0);
	}
	break;

	case LANCZOS2_STRETCHED:
	{
	    if (fabs(t) == 0.0)
		v = 1;
	    else if (fabs (t) > 2.5)
		v = 0.0;
	    else
	    {
		double tt = t * 2 / 2.5;
		v = sin (M_PI * tt) / (M_PI * tt) * sin (M_PI * tt / 2.0) / (M_PI * tt / 2.0);
	    }
	}
	break;
	
	case LANCZOS25:
	{
	    if (fabs(t) == 0.0)
		v = 1;
	    else if (fabs (t) > 2.5)
		v = 0.0;
	    else
	    {
		v = sin (M_PI * t) / (M_PI * t) * sin (M_PI * t / 2.5) / (M_PI * t / 2.5);
	    }
	}
	break;

	case LANCZOS3_SHRINK:
	{
	    if (fabs(t) == 0.0)
		v = 1;
	    else if (fabs (t) > 2.5)
		v = 0.0;
	    else
	    {
		double tt = t * 3 / 2.5;
		v = sin (M_PI * tt) / (M_PI * tt) * sin (M_PI * tt / 3.0) / (M_PI * tt / 3.0);
	    }
	}
	break;

	case LANCZOS3:
	{
	    if (fabs(t) == 0.0)
		v = 1;
	    else if (fabs (t) > 3)
		v = 0.0;
	    else
	    {
		double tt = t / 3;
		v = sin (M_PI * t) / (M_PI * t) * sin (M_PI * tt) / (M_PI * tt);
	    }
	}
	break;

	case LANCZOS4:
	{
	    if (fabs(t) == 0.0)
		v = 1;
	    else if (fabs (t) > 4)
		v = 0.0;
	    else
	    {
		double tt = t;
		v = sin (M_PI * tt) / (M_PI * tt) * sin (M_PI * tt / 4.0) / (M_PI * tt / 4.0);
	    }
	}
	break;

	case LANCZOS64:
	{
	    int i;
	    double tt;
	    double avg = 0;

	    for (i = 0; i < n_samples; ++i)
	    {
		tt = t - 0.5 +  i/(double)n_samples;
		
		if (fabs(tt) == 0.0)
		    avg += 1;
		else if (fabs (tt) > 2)
		    avg += 0.0;
		else
		{
		    avg += sin (M_PI * tt) / (M_PI * tt) * sin (M_PI * tt / 2.0) / (M_PI * tt / 2.0);
		}
	    }

	    v = avg / n_samples;

#if 0
	    if (v > 32.0)
		v = 32.0;
	    if (v < 0)
		v = 0;
#endif
	}
	break;
	
	case BOX:
	{
	    if (fabs (t) > 0.5)
		v = 0.0;
	    else
		v = 1.0;
	}
	break;

	case BOX15:
	{
	    if (fabs (t) > 0.75)
		v = 0.0;
	    else
		v = 0.66666667;
	}
	break;

	case LANCZOS1:
	{
	    if (fabs (t) == 0)
		v = 1;
	    else if (fabs(t) > 1.0)
		v = 0;
	    else
	    {
		double tt = t ;
		v = sin (M_PI * tt) / (M_PI * tt) * sin (M_PI * tt / 1.0) / (M_PI * tt / 1.0);
	    }
	    break;
	}

	case TANH:
	    v = (1 - pow (tanh (t), 16));
	    break;

	case TANH_SHRINK:
	    v = (1 - pow (tanh (1.1 * t), 128));
	    break;

	case QUADRATIC:
	    if (fabs (t) <= 0.5)
	    {
		v = -fabs(t) * fabs(t) + 3/4.0;
	    }
	    else if (fabs(t) <= 3/2.0)
	    {
		v = 0.5 * fabs(t) * fabs(t) - 3/2.0 * fabs(t) + 9/8.0;
	    }
	    else
		v = 0;
	    break;

	case FAST_CUBIC:
	    if (t < 0 && t >= -1)
		v = -2 * t * t * t - 3 * t * t + 1;
	    else if (t >= 0 && t <= 1)
		v = 2 * t * t * t - 3 * t * t + 1;
	    else
		v = 0;
	    break;

	case LINEAR:
	    if (fabs (t) < 1)
		v = 1 - fabs(t);
	    else
		v = 0;
	    break;

	case GAMMA:
	    if (fabs (t) >= 1)
		v = 1;
	    else 
		v = 0.5;

	    if (v > 1)
		printf ("wtf: %f\n", v);

	    if (fabs (t) > 0.5)
		v = 1 - v;
	    else
		v = 0 - v;
	    break;
	    
	case NONGAMMA:
	    if (fabs (t) >= 1)
		v = 1;
	    else
		v = 0.5;

	    v = pow (v, 2.2);

	    if (fabs (t) > 0.5)
		v = 1 - v;
	    else
		v = 0 - v;
	    break;
	    
	break;
	}

	test[i].re = v;
	test[i].im = iv;

	input[i] = test[i];
    }

    fft (test, N_SAMPLES);
#if 0
    complex_t *swapbuffer = malloc (N_SAMPLES * sizeof (complex_t));
    for (i = 0; i < N_SAMPLES; ++i)
    {
	swapbuffer[i] = test[i];

	if (i % 2 == 1)
	    swapbuffer[i].re = -swapbuffer[i].re;
    }

    for (i = 0; i < N_SAMPLES; ++i)
	test[i] = swapbuffer[i];
#endif
    
    sum = 0.0;
    double rsum = 0.0;
    double qs = 0.0;
    for (i = 0; i < N_SAMPLES; ++i)
    {
	double im, re;
	int idx = (i + N_SAMPLES/2) % N_SAMPLES;
	im = test[idx].im;
	re = test[idx].re;
	sum += input[idx].re;
	rsum += re;
	double dc = complex_mag (test[0]);
	dc = ((N_SAMPLES) / (last - first));
	double qnorm = ((re / dc) * (re / dc) + (im / dc) * (im/dc));
	double f = 16 * ((i - N_SAMPLES/2)/(last - first));
	double w = (f == 0)? 1 : MIN (1, pow (8/f, 8));
	qs += w * qnorm;
	printf ("%f %f %f %f %f %f %f %f %f %f\n",
		(i - N_SAMPLES/2)/(last - first),
		input[idx].re,

		re,
		im,

		rad_to_degree (complex_arg (test[idx])),
		complex_mag (test[idx]) / dc,
		sum,
		rsum,
	        w * qnorm,
		qs);
    }
    return 0;
}