various stuff

1 files changed, 97 insertions, 6 deletions

diff --git a/fft.c b/fft.c
index 08e999d..e160644 100644
--- a/fft.c
+++ b/fft.c
@@ -64,8 +64,8 @@ fft (complex_t *buffer, int n)
     fft (even, n/2);
     fft (odd, n/2);
 
-    principal.re = cos (2 * M_PI / n);
-    principal.im = sin (2 * M_PI / n);
+    principal.re = cos ( 2 * M_PI / n);
+    principal.im = sin ( 2 * M_PI / n);
 
     omega.re = 1.0;
     omega.im = 0;
@@ -75,34 +75,125 @@ 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);
     }
 
     free (temp);
 }
 
-#define N_SAMPLES 128
+#define N_SAMPLES (1 << 12)
 
 int
 main ()
 {
+    complex_t input[N_SAMPLES];
     complex_t test[N_SAMPLES];
+    double first = -512;
+    double last = 512;
     int i;
 
     for (i = 0; i < N_SAMPLES; ++i)
     {
-	double t = -10 + i * (20 / (double)N_SAMPLES);
-	double v = (t == 0.0) ? 1.0 : sin (5 * t * M_PI) / (5 * t * M_PI);
+	double sample_extent = (last - first) / N_SAMPLES;
+	double t = first + i * sample_extent + sample_extent / 2;
+
+	double v;
+
+	t = t * 1.0;
+	
+	if (t == 0.0)
+	    v = 1.0;
+	else
+	{
+#if 0
+	    /* lanczos 3 */
+	    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));
+#endif
+#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
+	    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;
+	    
+#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)) * (sin (4 / 5.0 * M_PI * t / 2) / (4/5.0 * M_PI * t / 2));
+#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;
 	test[i].im = 0.0;
+
+	input[i] = test[i];
     }
 
     fft (test, N_SAMPLES);
 
     for (i = 0; i < N_SAMPLES; ++i)
     {
-	printf ("%f %f\n", test[i].re, test[i].im);
+	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);
+	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));
     }
     return 0;
 }