Split fft into its own file

4 files changed, 427 insertions, 345 deletions

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..c0352a8
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,5 @@
+CFLAGS=-Wall
+LDFLAGS=-Wall -lm 
+
+filters:	fft.c fft.h filters.c
+	$(CC) -o filters fft.c filters.c $(LDFLAGS)
diff --git a/fft.c b/fft.c
index c26df96..3f4c6fe 100644
--- a/fft.c
+++ b/fft.c
@@ -1,65 +1,10 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
-
-typedef struct
-{
-    double re;
-    double im;
-} complex_t;
-
-static complex_t
-complex_mul (complex_t a, complex_t b)
-{
-    complex_t r;
-
-    r.re = a.re * b.re - a.im * b.im;
-    r.im = a.re * b.im + a.im * b.re;
-
-    return r;
-}
-
-static complex_t
-complex_add (complex_t a, complex_t b)
-{
-    complex_t r;
-
-    r.re = a.re + b.re;
-    r.im = a.im + b.im;
-
-    return r;
-}
-
-static complex_t
-complex_sub (complex_t a, complex_t b)
-{
-    complex_t r;
-
-    r.re = a.re - b.re;
-    r.im = a.im - b.im;
-
-    return r;
-}
-
-static double
-complex_mag (complex_t a)
-{
-    return sqrt (a.re * a.re + a.im * a.im);
-}
-static double
-complex_arg (complex_t a)
-{
-    return atan2 (a.im, a.re) + M_PI;
-}
-
-static double
-rad_to_degree (double theta)
-{
-    return theta * (360.0 / (2 * M_PI));
-}
+#include "fft.h"
 
 static void
-fft (complex_t *buffer, int n, complex_t principal)
+do_fft (complex_t *buffer, int n, complex_t principal)
 {
     complex_t *temp, *even, *odd;
     complex_t alpha;
@@ -78,8 +23,8 @@ fft (complex_t *buffer, int n, complex_t principal)
 	odd[i] = buffer[2 * i + 1];
     }
 
-    fft (even, n/2, complex_mul (principal, principal));
-    fft (odd, n/2, complex_mul (principal, principal));
+    do_fft (even, n/2, complex_mul (principal, principal));
+    do_fft (odd, n/2, complex_mul (principal, principal));
 
     alpha.re = 1.0;
     alpha.im = 0;
@@ -95,294 +40,13 @@ fft (complex_t *buffer, int n, complex_t principal)
     free (temp);
 }
 
-#define N_SAMPLES (1 << 12)
-
-typedef enum
+void
+fft (complex_t *buffer, int n)
 {
-    MITCHELL,
-    LANCZOS1,
-    LANCZOS2,
-    LANCZOS2_STRETCHED,
-    LANCZOS25,
-    LANCZOS3_SHRINK,
-    LANCZOS4,
-    LANCZOS64,
-    BOX,
-    BOX15,
-    TANH,
-    QUADRATIC,
-    FAST_CUBIC,
-    LINEAR,
-} filter_t;
-
-static filter_t filter = TANH;
-
-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 LANCZOS4:
-	{
-	    if (fabs(t) == 0.0)
-		v = 1;
-	    else if (fabs (t) > 4)
-		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 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), 8));
-	    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;
-	break;
-	}
-
-	test[i].re = v;
-	test[i].im = iv;
-
-	input[i] = test[i];
-    }
-
     complex_t p;
 
-    p.re = cos (2 * M_PI / N_SAMPLES);
-    p.im = sin (2 * M_PI / N_SAMPLES);
+    p.re = cos (2 * M_PI / n);
+    p.im = sin (2 * M_PI / n);
     
-    fft (test, N_SAMPLES, p);
-#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;
-    for (i = 0; i < N_SAMPLES; ++i)
-    {
-	double im, re;
-	int idx = (i + N_SAMPLES/2) % N_SAMPLES;
-	double spacing = 1/((double)N_SAMPLES);
-	double T = (last - first);
-	double period = (last - first)/((double)N_SAMPLES);
-	double sample_width = (M_PI * 2) / N_SAMPLES;
-	im = test[idx].im;
-	re = test[idx].re;
-	double pos = sample_width / 2 + (i * sample_width);
-	double sample_rate = N_SAMPLES / (last - first);
-	sum += input[idx].re;
-	printf ("%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]),
-		sum);
-    }
-    return 0;
+    do_fft (buffer, n, p);
 }
diff --git a/fft.h b/fft.h
new file mode 100644
index 0000000..2973deb
--- /dev/null
+++ b/fft.h
@@ -0,0 +1,58 @@
+typedef struct
+{
+    double re;
+    double im;
+} complex_t;
+
+static inline complex_t
+complex_mul (complex_t a, complex_t b)
+{
+    complex_t r;
+
+    r.re = a.re * b.re - a.im * b.im;
+    r.im = a.re * b.im + a.im * b.re;
+
+    return r;
+}
+
+static inline complex_t
+complex_add (complex_t a, complex_t b)
+{
+    complex_t r;
+
+    r.re = a.re + b.re;
+    r.im = a.im + b.im;
+
+    return r;
+}
+
+static inline complex_t
+complex_sub (complex_t a, complex_t b)
+{
+    complex_t r;
+
+    r.re = a.re - b.re;
+    r.im = a.im - b.im;
+
+    return r;
+}
+
+static inline double
+complex_mag (complex_t a)
+{
+    return sqrt (a.re * a.re + a.im * a.im);
+}
+static inline double
+complex_arg (complex_t a)
+{
+    return atan2 (a.im, a.re) + M_PI;
+}
+
+static inline double
+rad_to_degree (double theta)
+{
+    return theta * (360.0 / (2 * M_PI));
+}
+
+void
+fft (complex_t *buffer, int n);
diff --git a/filters.c b/filters.c
new file mode 100644
index 0000000..17f3537
--- /dev/null
+++ b/filters.c
@@ -0,0 +1,355 @@
+#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;
+}