From cb30c96bc67d195bde2586161d383cf14d1cae8c Mon Sep 17 00:00:00 2001 From: Andreas Eversberg Date: Thu, 16 Nov 2017 18:25:03 +0100 Subject: Restructure: Move fft from common code to 'libfft' --- src/common/fft.c | 96 -------------------------------------------------- src/common/fft.h | 3 -- src/libfft/Makefile.am | 6 ++++ src/libfft/fft.c | 96 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/libfft/fft.h | 3 ++ 5 files changed, 105 insertions(+), 99 deletions(-) delete mode 100644 src/common/fft.c delete mode 100644 src/common/fft.h create mode 100644 src/libfft/Makefile.am create mode 100644 src/libfft/fft.c create mode 100644 src/libfft/fft.h diff --git a/src/common/fft.c b/src/common/fft.c deleted file mode 100644 index fb819e1..0000000 --- a/src/common/fft.c +++ /dev/null @@ -1,96 +0,0 @@ -/* Fast Fourier Transformation (FFT) - * - * (C) 2016 by Andreas Eversberg - * 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 . - */ - -#include -#include -#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/common/fft.h b/src/common/fft.h deleted file mode 100644 index 8387a8e..0000000 --- a/src/common/fft.h +++ /dev/null @@ -1,3 +0,0 @@ - -void fft_process(int dir, int m, double *x, double *y); - 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 + * 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 . + */ + +#include +#include +#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); + -- cgit v1.2.3