misc

1 files changed, 108 insertions, 27 deletions

diff --git a/fft.c b/fft.c
index e160644..52dbb0f 100644
--- a/fft.c
+++ b/fft.c
@@ -41,11 +41,28 @@ complex_sub (complex_t a, complex_t b)
     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));
+}
+
 static void
-fft (complex_t *buffer, int n)
+fft (complex_t *buffer, int n, complex_t principal)
 {
     complex_t *temp, *even, *odd;
-    complex_t principal, omega;
+    complex_t omega;
     int i;
     
     if (n == 1)
@@ -61,11 +78,8 @@ fft (complex_t *buffer, int n)
 	odd[i] = buffer[2 * i + 1];
     }
 
-    fft (even, n/2);
-    fft (odd, n/2);
-
-    principal.re = cos ( 2 * M_PI / n);
-    principal.im = sin ( 2 * M_PI / n);
+    fft (even, n/2, complex_mul (principal, principal));
+    fft (odd, n/2, complex_mul (principal, principal));
 
     omega.re = 1.0;
     omega.im = 0;
@@ -75,14 +89,6 @@ fft (complex_t *buffer, int n)
 	buffer[i]	= complex_add (even[i], complex_mul (omega, odd[i]));
 	buffer[i + n/2] = complex_sub (even[i], complex_mul (omega, odd[i]));
 
-#if 0
-	if (i%2 == -0)
-	{
-	    buffer[i].re = - buffer[i].re;
-	    buffer[i + n/2].re = buffer[i + n/2].re;
-	}
-#endif
-    
 	omega = complex_mul (omega, principal);
     }
 
@@ -96,15 +102,15 @@ main ()
 {
     complex_t input[N_SAMPLES];
     complex_t test[N_SAMPLES];
-    double first = -512;
-    double last = 512;
+    double first = -64;
+    double last = 64;
     int i;
 
     for (i = 0; i < N_SAMPLES; ++i)
     {
 	double sample_extent = (last - first) / N_SAMPLES;
 	double t = first + i * sample_extent + sample_extent / 2;
-
+	double iv = 0.0;
 	double v;
 
 	t = t * 1.0;
@@ -114,12 +120,33 @@ main ()
 	else
 	{
 #if 0
-	    /* lanczos 3 */
+	    /* 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)) * (sin (M_PI * t / 2.5) / (M_PI * t / 2.5));
+		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)
@@ -127,18 +154,60 @@ main ()
 	    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
-	    /* slightly stretched lanczos 2 */
-	    if (fabs (t) > 2.5)
+	    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/5.0 * M_PI * t) / (4/5.0 * M_PI * t)) * (sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2));
+		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)
@@ -173,12 +242,17 @@ main ()
 	}
 
 	test[i].re = v;
-	test[i].im = 0.0;
+	test[i].im = iv;
 
 	input[i] = test[i];
     }
 
-    fft (test, N_SAMPLES);
+    complex_t p;
+
+    p.re = cos ( 2 * M_PI / N_SAMPLES);
+    p.im = sin ( 2 * M_PI / N_SAMPLES);
+    
+    fft (test, N_SAMPLES, p);
 
     for (i = 0; i < N_SAMPLES; ++i)
     {
@@ -192,8 +266,15 @@ 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\n", (i - N_SAMPLES/2)/(last - first), re , im,
-		((atan2(im, re) + M_PI)/(2 * M_PI)) * 360, sqrt (re * re + im * im));
+	printf ("%f %f %f %f %f %f\n",
+		(i - N_SAMPLES/2)/(last - first),
+		input[i].re,
+
+		re,
+		im,
+
+		rad_to_degree (complex_arg (test[idx])),
+		complex_mag (test[idx]));
     }
     return 0;
 }