From 897f5fbac08776364dc93cd350e0cd6d915960f7 Mon Sep 17 00:00:00 2001
From: Søren Sandmann Pedersen <ssp@redhat.com>
Date: Fri, 24 Sep 2010 23:13:02 -0400
Subject: Add inverse FFT

---
 Makefile   |  8 +++++++-
 fft-test.c | 41 +++++++++++++++++++++++++++++++++++++++++
 fft.c      | 31 +++++++++++++++++++++++++++----
 fft.h      | 11 +++++++++++
 4 files changed, 86 insertions(+), 5 deletions(-)
 create mode 100644 fft-test.c

diff --git a/Makefile b/Makefile
index c0352a8..519e05d 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,11 @@
 CFLAGS=-Wall
 LDFLAGS=-Wall -lm 
 
-filters:	fft.c fft.h filters.c
+all: filters fft-test
+
+filters: fft.c fft.h filters.c
 	$(CC) -o filters fft.c filters.c $(LDFLAGS)
+
+fft-test: fft.c fft.h fft-test.c
+	$(CC) -o fft-test fft.c fft-test.c $(LDFLAGS)
+
diff --git a/fft-test.c b/fft-test.c
new file mode 100644
index 0000000..98667a5
--- /dev/null
+++ b/fft-test.c
@@ -0,0 +1,41 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "fft.h"
+
+#define N_SAMPLES (1 << 14)
+
+int
+main ()
+{
+    complex_t orig[N_SAMPLES];
+    complex_t buffer[N_SAMPLES];
+
+    int i;
+
+    for (i = 0; i < N_SAMPLES; ++i)
+    {
+	orig[i].re = drand48() * 10000;
+	orig[i].im = drand48() * 10000;
+
+	buffer[i] = orig[i];
+    }
+
+    fft (buffer, N_SAMPLES);
+    ifft (buffer, N_SAMPLES);
+    
+    for (i = 0; i < N_SAMPLES; ++i)
+    {
+	if (fabs (buffer[i].re - orig[i].re ) > 0.02)
+	    printf ("at %d: %f vs %f\n", i, buffer[i].re, orig[i].re);
+	
+	if (fabs (buffer[i].im - orig[i].im ) > 0.02)
+	    printf ("at %d: %f vs %f\n", i, buffer[i].im, orig[i].im);
+	
+	assert (fabs (buffer[i].re - orig[i].re) < 0.02);
+	assert (fabs (buffer[i].im - orig[i].im) < 0.02);
+    }
+
+    return 0;
+}
diff --git a/fft.c b/fft.c
index b8af1f4..f903dd4 100644
--- a/fft.c
+++ b/fft.c
@@ -4,12 +4,13 @@
 #include "fft.h"
 
 static void
-do_fft (complex_t *buffer, int n, complex_t principal)
+do_fft (complex_t *buffer, int n, complex_t principal, int is_inverse)
 {
     complex_t *temp, *even, *odd;
     complex_t alpha;
     complex_t p2;
     int i;
+    double inv;
     
     if (n == 1)
 	return;
@@ -26,17 +27,28 @@ do_fft (complex_t *buffer, int n, complex_t principal)
 
     p2 = complex_mul (principal, principal);
     
-    do_fft (even, n/2, p2);
-    do_fft (odd, n/2, p2);
+    do_fft (even, n/2, p2, FALSE);
+    do_fft (odd, n/2, p2, FALSE);
 
     alpha.re = 1.0;
     alpha.im = 0;
+
+    if (is_inverse)
+	inv = 1.0 / n;
     
     for (i = 0; i < n/2; ++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]));
 
+	if (is_inverse)
+	{
+	    buffer[i].re *= inv;
+	    buffer[i].im *= inv;
+	    buffer[i + n/2].re *= inv;
+	    buffer[i + n/2].im *= inv;
+	}
+	
 	alpha = complex_mul (alpha, principal);
     }
 
@@ -51,5 +63,16 @@ fft (complex_t *buffer, int n)
     p.re = cos (2 * M_PI / n);
     p.im = sin (2 * M_PI / n);
     
-    do_fft (buffer, n, p);
+    do_fft (buffer, n, p, FALSE);
+}
+
+void
+ifft (complex_t *buffer, int n)
+{
+    complex_t p;
+
+    p.re = cos (- (2 * M_PI / n));
+    p.im = sin (- (2 * M_PI / n));
+
+    do_fft (buffer, n, p, TRUE);
 }
diff --git a/fft.h b/fft.h
index 2973deb..8a4e20e 100644
--- a/fft.h
+++ b/fft.h
@@ -1,3 +1,11 @@
+#ifndef FALSE
+#define FALSE (0)
+#endif
+
+#ifndef TRUE
+#define TRUE (1)
+#endif
+
 typedef struct
 {
     double re;
@@ -56,3 +64,6 @@ rad_to_degree (double theta)
 
 void
 fft (complex_t *buffer, int n);
+
+void
+ifft (complex_t *buffer, int n);
-- 
cgit v1.2.3