Restructure: Move fft from common code to 'libfft'
This commit is contained in:
parent
06c07146d9
commit
058f641789
|
@ -0,0 +1,6 @@
|
|||
AM_CPPFLAGS = -Wall -Wextra -g $(all_includes)
|
||||
|
||||
noinst_LIBRARIES = libfft.a
|
||||
|
||||
libfft_a_SOURCES = \
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
|
@ -0,0 +1,3 @@
|
|||
|
||||
void fft_process(int dir, int m, double *x, double *y);
|
||||
|
Loading…
Reference in New Issue