1 files changed, 96 insertions, 0 deletions
diff --git a/src/common/fft.c b/src/common/fft.c
new file mode 100644
index 0000000..fb819e1
--- /dev/null
+++ b/src/common/fft.c
@@ -0,0 +1,96 @@
+/* Fast Fourier Transformation (FFT)
+ *
+ * (C) 2016 by Andreas Eversberg <jolly@eversberg.eu>
+ * All Rights Reserved
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+#include <string.h>
+#include "fft.h"
+
+/*
+ * Code based closely to work by Paul Bourke
+ *
+ * This computes an in-place complex-to-complex FFT 
+ * x and y are the real and imaginary arrays of 2^m points.
+ * dir =  1 gives forward transform
+ * dir = -1 gives reverse transform 
+ */
+void fft_process(int dir, int m, double *x, double *y)
+{
+	int	n, i, i1, j, k, i2, l, l1, l2;
+	double	c1, c2, tx, ty, t1, t2, u1, u2, z;
+
+	/* Calculate the number of points */
+	n = 1 << m;
+
+	/* Do the bit reversal */
+	i2 = n >> 1;
+	j = 0;
+	for (i = 0; i < n - 1; i++) {
+		if (i < j) {
+			tx = x[i];
+			ty = y[i];
+			x[i] = x[j];
+			y[i] = y[j];
+			x[j] = tx;
+			y[j] = ty;
+		}
+		k = i2;
+		while (k <= j) {
+			j -= k;
+			k >>= 1;
+		}
+		j += k;
+	}
+
+	/* Compute the FFT */
+	c1 = -1.0; 
+	c2 = 0.0;
+	l2 = 1;
+	for (l = 0; l < m; l++) {
+		l1 = l2;
+		l2 <<= 1;
+		u1 = 1.0; 
+		u2 = 0.0;
+		for (j = 0; j < l1; j++) {
+			for (i = j; i < n; i += l2) {
+				i1 = i + l1;
+				t1 = u1 * x[i1] - u2 * y[i1];
+				t2 = u1 * y[i1] + u2 * x[i1];
+				x[i1] = x[i] - t1; 
+				y[i1] = y[i] - t2;
+				x[i] += t1;
+				y[i] += t2;
+			}
+			z =  u1 * c1 - u2 * c2;
+			u2 = u1 * c2 + u2 * c1;
+			u1 = z;
+		}
+		c2 = sqrt((1.0 - c1) / 2.0);
+		if (dir == 1) 
+			c2 = -c2;
+		c1 = sqrt((1.0 + c1) / 2.0);
+	}
+
+	/* Scaling for forward transform */
+	if (dir == 1) {
+		for (i = 0; i < n; i++) {
+			x[i] /= n;
+			y[i] /= n;
+		}
+	}
+}