diff options
author | Andreas Eversberg <jolly@eversberg.eu> | 2017-11-16 18:25:03 +0100 |
---|---|---|
committer | Andreas Eversberg <jolly@eversberg.eu> | 2017-11-25 19:23:57 +0100 |
commit | cb30c96bc67d195bde2586161d383cf14d1cae8c (patch) | |
tree | 8e494fc2dc38be58a3e4c2198cd14b96cdb8f08d /src/libfft | |
parent | 40c24a70b2a295073428ec78c281e86a7a84f7ca (diff) |
Restructure: Move fft from common code to 'libfft'
Diffstat (limited to 'src/libfft')
-rw-r--r-- | src/libfft/Makefile.am | 6 | ||||
-rw-r--r-- | src/libfft/fft.c | 96 | ||||
-rw-r--r-- | src/libfft/fft.h | 3 |
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); + |