aboutsummaryrefslogtreecommitdiffstats
path: root/src/libfft
diff options
context:
space:
mode:
Diffstat (limited to 'src/libfft')
-rw-r--r--src/libfft/Makefile.am6
-rw-r--r--src/libfft/fft.c96
-rw-r--r--src/libfft/fft.h3
3 files changed, 105 insertions, 0 deletions
diff --git a/src/libfft/Makefile.am b/src/libfft/Makefile.am
new file mode 100644
index 0000000..b6cbfc9
--- /dev/null
+++ b/src/libfft/Makefile.am
@@ -0,0 +1,6 @@
+AM_CPPFLAGS = -Wall -Wextra -g $(all_includes)
+
+noinst_LIBRARIES = libfft.a
+
+libfft_a_SOURCES = \
+ fft.c
diff --git a/src/libfft/fft.c b/src/libfft/fft.c
new file mode 100644
index 0000000..fb819e1
--- /dev/null
+++ b/src/libfft/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;
+ }
+ }
+}
diff --git a/src/libfft/fft.h b/src/libfft/fft.h
new file mode 100644
index 0000000..8387a8e
--- /dev/null
+++ b/src/libfft/fft.h
@@ -0,0 +1,3 @@
+
+void fft_process(int dir, int m, double *x, double *y);
+