Shift to make it symmetric around sample 0

1 files changed, 21 insertions, 123 deletions

diff --git a/fft.c b/fft.c
index 52dbb0f..76af8d5 100644
--- a/fft.c
+++ b/fft.c
@@ -62,7 +62,7 @@ static void
 fft (complex_t *buffer, int n, complex_t principal)
 {
     complex_t *temp, *even, *odd;
-    complex_t omega;
+    complex_t alpha;
     int i;
     
     if (n == 1)
@@ -81,15 +81,15 @@ fft (complex_t *buffer, int n, complex_t principal)
     fft (even, n/2, complex_mul (principal, principal));
     fft (odd, n/2, complex_mul (principal, principal));
 
-    omega.re = 1.0;
-    omega.im = 0;
+    alpha.re = 1.0;
+    alpha.im = 0;
     
     for (i = 0; i < n/2; ++i)
     {
-	buffer[i]	= complex_add (even[i], complex_mul (omega, odd[i]));
-	buffer[i + n/2] = complex_sub (even[i], complex_mul (omega, odd[i]));
+	buffer[i]	= complex_add (even[i], complex_mul (alpha, odd[i]));
+	buffer[i + n/2] = complex_sub (even[i], complex_mul (alpha, odd[i]));
 
-	omega = complex_mul (omega, principal);
+	alpha = complex_mul (alpha, principal);
     }
 
     free (temp);
@@ -109,7 +109,8 @@ main ()
     for (i = 0; i < N_SAMPLES; ++i)
     {
 	double sample_extent = (last - first) / N_SAMPLES;
-	double t = first + i * sample_extent + sample_extent / 2;
+	double sn = (i + N_SAMPLES/2) % N_SAMPLES - N_SAMPLES/2;
+	double t = sn * sample_extent;
 	double iv = 0.0;
 	double v;
 
@@ -119,126 +120,10 @@ main ()
 	    v = 1.0;
 	else
 	{
-#if 0
-	    /* slightly stretched lanczos 2 */
-	    if (fabs (t) > 2.5)
-		v = 0.0;
-	    else
-		v = (sin (4/5.0 * M_PI * t) / (4/5.0 * M_PI * t)) * pow ((sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2)), 1);
-#endif
-
-#if 0
 	    if (fabs (t) > 0.5)
 		v = 0.0;
 	    else
 		v = 1.0;
-#endif
-	    
-#if 0
-	    /* lanczos 3 */
-	    if (fabs (t) > 4)
-		v = 0.0;
-	    else
-		v = (sin (3.0/4 * M_PI * t) / (3.0/4 * M_PI * t)) * (sin (M_PI * 3.0/4 * t / 3) / (M_PI * 3.0/4 * t / 3));
-#endif
-
-	    /* gaussian window */
-	    if (fabs (t) > 4)
-		v = 0.0;
-	    else
-		v = (sin (M_PI * t) / (M_PI * t)) * (exp (-0.5 * t * t));
-#if 0
-	    /* raised cosine */
-	    if (fabs (t) > 2.5)
-		v = 0.0;
-	    else
-		v = (sin (4/5.0 * M_PI * t) / (4/5.0 * M_PI * t)) * (pow (0.5 + 0.5 * cos (M_PI * t / 2.5), 0.55));
-#endif
-#if 0
-	    v = pow (0.5 + 0.5 * cos (M_PI * t / 2.5), 0.55);
-	    if (fabs (t) > 1.5)
-		v = 0.0;
-	    else
-		v = 1.0;
-#endif
-	    
-#if 0
-	    v = sin (2 * M_PI * t * 1/100000.0 * pow (drand48(), 0.5));
-#endif
-#if 0
-	    int j;
-
-	    v = 0.0;
-	    t = fmod ((drand48() * 1024 - 512.0), 17) * (512.0/17);
-	    v = t; //512 * pow (fabs (t)/512.0, 0.5 * drand48()) * sin (2 * M_PI * (fabs (t)/512.0) * t);
-#endif
-
-#if 0
-	    if (fabs (t) > 10)
-		v = 0.0;
-	    else
-		v = (sin (4/4.5 * M_PI * t) / (4/4.5 * M_PI * t)) * (sin (4/4.5 * M_PI * t / 2.0) / (4/4.5 * M_PI * t / 2.0));
-#endif
-#if 0
-	    if (fabs (t) < 0.01)
-		0.01;
-#endif
-
-	    double lol = drand48();
-	    
-#if 0
-	    if (fabs (t) > 4)
-		v = 0;
-	    else
-		v = sin (2 * M_PI * t);
-#endif
-	    
-#if 0
-#define EXPO	0.6
-	    
-	    v = (sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2));
-	    if (v < 0)
-	    {
-		v = cos (M_PI * EXPO) * pow (-v, EXPO);
-
-		iv = sin (M_PI * EXPO);
-		
-	    }
-	    else
-	    {
-		v = pow (v, EXPO);
-	    }
-#endif
-#if 0
-	    if (fabs (t) > 2.5)
-		v = 0.0;
-	    else
-		v = (sin (M_PI * t) / (M_PI * t)) * (0.5 * (1 - cos (((t*M_PI)/2.5 - M_PI))));
-#endif
-#if 0
-	    if (fabs (t) > 1.5)
-		v = 0.0;
-	    else
-		v = (sin (M_PI * t) / (M_PI * t)) * (sin (M_PI * (t/1.5)) / M_PI / (t/1.5));
-#endif
-#if 0
-	    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));
-#endif
-#if 0
-	    if (fabs (t) > 3)
-		v = 0.0;
-	    else
-		v = (sin (M_PI * t) / (M_PI * t)) * (sin (M_PI * (t/3)) / M_PI / (t/3));
-#endif
-#if 0
-	    if (fabs (t) > 2)
-		v = 0.0;
-	    else
-		v = (sin (M_PI * t) / (M_PI * t));
-#endif
 	}
 
 	test[i].re = v;
@@ -253,8 +138,21 @@ main ()
     p.im = sin ( 2 * M_PI / N_SAMPLES);
     
     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
+    
+    for (i = 0; i < N_SAMPLES; ++i)
     {
 	double im, re;
 	int idx = (i + N_SAMPLES/2) % N_SAMPLES;