initial

1 files changed, 108 insertions, 0 deletions

diff --git a/fft.c b/fft.c
new file mode 100644
index 0000000..08e999d
--- /dev/null
+++ b/fft.c
@@ -0,0 +1,108 @@
+#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 void
+fft (complex_t *buffer, int n)
+{
+    complex_t *temp, *even, *odd;
+    complex_t principal, omega;
+    int i;
+    
+    if (n == 1)
+	return;
+
+    temp = malloc (n * sizeof (complex_t));
+    even = temp;
+    odd = temp + n/2;
+    
+    for (i = 0; i < n/2; ++i)
+    {
+	even[i] = buffer[2 * i];
+	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);
+
+    omega.re = 1.0;
+    omega.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]));
+
+	omega = complex_mul (omega, principal);
+    }
+
+    free (temp);
+}
+
+#define N_SAMPLES 128
+
+int
+main ()
+{
+    complex_t test[N_SAMPLES];
+    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);
+
+	test[i].re = v;
+	test[i].im = 0.0;
+    }
+
+    fft (test, N_SAMPLES);
+
+    for (i = 0; i < N_SAMPLES; ++i)
+    {
+	printf ("%f %f\n", test[i].re, test[i].im);
+    }
+    return 0;
+}