Filters

1 files changed, 220 insertions, 10 deletions

diff --git a/fft.c b/fft.c
index 76af8d5..c26df96 100644
--- a/fft.c
+++ b/fft.c
@@ -97,6 +97,26 @@ fft (complex_t *buffer, int n, complex_t principal)
 
 #define N_SAMPLES (1 << 12)
 
+typedef enum
+{
+    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 ()
 {
@@ -105,26 +125,213 @@ main ()
     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 sample_extent = (last - first) / N_SAMPLES;
 	double sn = (i + N_SAMPLES/2) % N_SAMPLES - N_SAMPLES/2;
-	double t = sn * sample_extent;
+	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;
 	
-	if (t == 0.0)
-	    v = 1.0;
-	else
+	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;
@@ -134,8 +341,8 @@ main ()
 
     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_SAMPLES);
+    p.im = sin (2 * M_PI / N_SAMPLES);
     
     fft (test, N_SAMPLES, p);
 #if 0
@@ -152,6 +359,7 @@ main ()
 	test[i] = swapbuffer[i];
 #endif
     
+    sum = 0.0;
     for (i = 0; i < N_SAMPLES; ++i)
     {
 	double im, re;
@@ -164,15 +372,17 @@ main ()
 	re = test[idx].re;
 	double pos = sample_width / 2 + (i * sample_width);
 	double sample_rate = N_SAMPLES / (last - first);
-	printf ("%f %f %f %f %f %f\n",
+	sum += input[idx].re;
+	printf ("%f %f %f %f %f %f %f\n",
 		(i - N_SAMPLES/2)/(last - first),
-		input[i].re,
+		input[idx].re,
 
 		re,
 		im,
 
 		rad_to_degree (complex_arg (test[idx])),
-		complex_mag (test[idx]));
+		complex_mag (test[idx]),
+		sum);
     }
     return 0;
 }