---
configure.ac | 8 ++-
include/fft.h | 6 +-
src/libfft/fftest.c | 24 ++---
src/libfft/fftfast.c | 263 ++++----------------------------------------------
4 files changed, 41 insertions(+), 260 deletions(-)
diff --git a/configure.ac b/configure.ac
index ad155f6..9708b15 100644
--- a/configure.ac
+++ b/configure.ac
@@ -916,6 +916,9 @@ AM_PROG_MKDIR_P
AM_PROG_LEX
AC_PROG_YACC
+# pkg-config
+PKG_PROG_PKG_CONFIG
+
# libtool's configuration check has a bug that causes a /lib/cpp
# sanity check failure if a C++ compiler is not installed. This makes
# the sanity test pass regardless of whether there is a c++ compiler.
@@ -1568,6 +1571,9 @@ AC_CHECK_LIB(stdc++, main, stdcxx_link_works=yes ;
LIBSTDCXX="-lstdc++",
)]
)
+dnl Search and check for FFTW3
+PKG_CHECK_MODULES([FFTW3], [fftw3])
+
dnl *************************
dnl *** Check for headers ***
dnl *************************
@@ -4095,7 +4101,7 @@ PKG_LIBS="${LIBSOCKET} ${LIBNSL} ${LIBNETWORK}"
SYSV_LIBS=
TERMIO_LIBS="${TERMLIB}"
CURSOR_LIBS="${TERMLIB}"
-FFT_LIBS="${LIBM}"
+FFT_LIBS="${LIBM} ${FFTW3_LIBS}"
ANALYZE_LIBS="${RT} ${RT_LIBS} ${BN} ${BN_LIBS} ${BU} ${BU_LIBS}"
BU_LIBS="${TCL} ${PNG} ${LIBZ} ${LIBM} ${LIBMALLOC} ${LIBTHREAD}"
BN_LIBS="${BU} ${BU_LIBS} ${TCL} ${LIBM}"
diff --git a/include/fft.h b/include/fft.h
index 0cc46e5..c23dd00 100644
--- a/include/fft.h
+++ b/include/fft.h
@@ -24,6 +24,7 @@
#include "common.h"
#include <math.h>
+#include <fftw3.h>
#ifndef M_PI
# define M_PI 3.14159265358979323846264338328
@@ -57,10 +58,11 @@ FFT_EXPORT extern void splitdit(int N, int M);
FFT_EXPORT extern void ditsplit(int n /* length */, int m /* n = 2^m */);
FFT_EXPORT extern void rfft(double *X, int N);
FFT_EXPORT extern void irfft(double *X, int n);
-FFT_EXPORT extern void cfft(COMPLEX *dat, int num);
-FFT_EXPORT extern void icfft(COMPLEX *dat, int num);
FFT_EXPORT extern void cdiv(COMPLEX *result, COMPLEX *val1, COMPLEX *val2);
+FFT_EXPORT extern void cfft(fftw_complex *dat, int num);
+FFT_EXPORT extern void icfft(fftw_complex *dat, int num);
+
/* These should come from a generated header, but until
* CMake is live we'll just add the ones used by our current code */
FFT_EXPORT extern void rfft256(register double X[]);
diff --git a/src/libfft/fftest.c b/src/libfft/fftest.c
index f329118..72571af 100644
--- a/src/libfft/fftest.c
+++ b/src/libfft/fftest.c
@@ -22,7 +22,6 @@
* Complex Number and FFT Library
*
*/
-
#include "common.h"
#include <stdio.h>
@@ -30,36 +29,33 @@
#include "fft.h"
-/***** TEST ROUTINES *****/
-COMPLEX data[64];
-
+#define SIZE 64
void
-display(COMPLEX *dat, int num)
+display(fftw_complex *dat, int num)
{
int i;
-
for (i = 0; i < num; i++) {
printf("%3d : ", i);
- printf("%f, %f\n", dat[i].re, dat[i].im);
+ printf("%f, %f\n", dat[i][0], dat[i][1]);
}
}
int
main(int ac, char *av[])
{
- extern void cfft(COMPLEX *, int);
- extern void icfft(COMPLEX *, int);
-
int i;
+ fftw_complex *data;
+
+ data = fftw_malloc(sizeof(fftw_complex) * SIZE);
if (ac > 1)
printf("Usage: %s\n", av[0]);
for (i = 0; i < 64; i++) {
- data[i].re = sin((double)2.0 * M_PI * i / 64.0);
- data[i].re += 3 * cos((double)2.0 * M_PI * i / 32.0);
- data[i].im = (double)0.0;
+ data[i][0] = sin((double)2.0 * M_PI * i / 64.0);
+ data[i][0] += 3 * cos((double)2.0 * M_PI * i / 32.0);
+ data[i][1] = (double)0.0;
}
printf("Original Data:\n\n");
@@ -73,6 +69,8 @@ main(int ac, char *av[])
printf("\n\nInversed Data:\n\n");
display(data, 64);
+ fftw_free(data);
+
return 0;
}
diff --git a/src/libfft/fftfast.c b/src/libfft/fftfast.c
index b9b1aa7..41543be 100644
--- a/src/libfft/fftfast.c
+++ b/src/libfft/fftfast.c
@@ -35,272 +35,47 @@
#include "common.h"
#include <stdlib.h>
-#include <stdio.h> /* for stderr */
#include "fft.h"
-
-#define MAXSIZE 65536 /* Needed for sin/cos tables */
-int _init_size = 0; /* Internal: shows last initialized size */
-
-void scramble(int numpoints, COMPLEX *dat);
-void butterflies(int numpoints, int inverse, COMPLEX *dat);
-int init_sintab(int size);
-
+static void run_fft(fftw_complex *dat, int size, int dir);
/*
* Forward Complex Fourier Transform
*/
void
-cfft(COMPLEX *dat, int num)
+cfft(fftw_complex *dat, int num)
{
- /* Check for trig table initialization */
- if (num != _init_size) {
- if (init_sintab(num) == 0) {
- /* Can't do requested size */
- return;
- }
- }
-
- scramble(num, dat);
- butterflies(num, -1, dat);
+ run_fft(dat, num, FFTW_FORWARD);
}
/*
* Inverse Complex Fourier Transform
*/
void
-icfft(COMPLEX *dat, int num)
-{
- /* Check for trig table initialization */
- if (num != _init_size) {
- if (init_sintab(num) == 0) {
- /* Can't do requested size */
- return;
- }
- }
-
- scramble(num, dat);
- butterflies(num, 1, dat);
-}
-
-/******************* INTERNAL FFT ROUTINES ********************/
-
-/* The trig tables */
-double *sintab;
-double *costab;
-
-/*
- * I N I T _ S I N T A B
- *
- * Internal routine to initialize the sine/cosine table for
- * transforms of a given size. Checks size for power of two
- * and within table limits.
- *
- * Note that once initialized for one size it ready for one
- * smaller than that also, but it is convenient to do power of
- * two checking here so we change the _init_size every time
- * (We *could* pick up where ever we left off by keeping a
- * _max_init_size but forget that for now).
- *
- * Note that we need sin and cos values for +/- (M_PI * (m / col))
- * where col = 1, 2, 4, ..., N/2
- * m = 0, 1, 2, ..., col-1
- *
- * Thus we can subscript by: table[(m / col) * N/2]
- * or with twice as many values by: table[m + col]
- * We chose the later. (but N.B. this doesn't allow sub
- * _init_size requests to use existing numbers!)
- */
-int
-init_sintab(int size)
-{
- double theta;
- int col, m;
-
- /*
- * Check whether the requested size is within our compiled
- * limit and make sure it's a power of two.
- */
- if (size > MAXSIZE) {
- fprintf(stderr, "fft: Only compiled for max size of %d\n", MAXSIZE);
- fprintf(stderr, "fft: Can't do the requested %d\n", size);
- return 0;
- }
- for (m = size; (m & 1) == 0; m >>= 1)
- ;
- if (m != 1) {
- fprintf(stderr, "fft: Can only do powers of two, not %d\n", size);
- fprintf(stderr, "fft: What do you think this is, a Winograd
transform?\n");
- return 0;
- }
-
- /* Get some buffer space */
- if (sintab != NULL) free(sintab);
- if (costab != NULL) free(costab);
- /* should not use bu_calloc() as libfft is not dependant upon libbu */
- sintab = (double *)calloc(sizeof(*sintab), size);
- costab = (double *)calloc(sizeof(*costab), size);
-
- /*
- * Size is okay. Set up tables.
- */
- for (col = 1; col < size; col <<= 1) {
- for (m = 0; m < col; m++) {
- theta = M_PI * (double)m / (double)col;
- sintab[ m + col ] = sin(theta);
- costab[ m + col ] = cos(theta);
- }
- }
-
- /*
- * Mark size and return success.
- */
- _init_size = size;
-/* fprintf(stderr, "fft: table init, size = %d\n", size);*/
- return 1;
-}
-
-/*
- * This section implements the Cooley-Tukey Complex
- * Fourier transform.
- * Reference: Nov 84 Dr. Dobbs [#97], and
- * "Musical Applications of Microprocessors", Hal Chamberlin.
- */
-
-/*
- * SCRAMBLE - put data in bit reversed order
- *
- * Note: Could speed this up with pointers if necessary,
- * but the butterflies take much longer.
- */
-void
-scramble(int numpoints, COMPLEX *dat)
+icfft(fftw_complex *dat, int num)
{
- register int i, j, m;
- COMPLEX temp;
-
- j = 0;
- for (i = 0; i < numpoints; i++, j += m) {
- if (i < j) {
- /* Switch nodes i and j */
- temp.re = dat[j].re;
- temp.im = dat[j].im;
- dat[j].re = dat[i].re;
- dat[j].im = dat[i].im;
- dat[i].re = temp.re;
- dat[i].im = temp.im;
- }
- m = numpoints/2;
- while (m-1 < j) {
- j -= m;
- m = (m + 1) / 2;
- }
- }
+ run_fft(dat, num, FFTW_BACKWARD);
}
-void
-butterflies(int numpoints, int inverse, COMPLEX *dat)
+static void run_fft(fftw_complex *dat, int num, int dir)
{
- register COMPLEX *node1, *node2;
- register int step, column, m;
- COMPLEX w, temp;
-
- /*
- * For each column of the butterfly
- */
- for (column = 1; column < numpoints; column = step) {
- step = 2 * column; /* step is size of "cross-hatch" */
- /*
- * For each principle value of W (roots on units).
- */
- for (m = 0; m < column; m++) {
- /*
- * Do these by table lookup:
- * theta = M_PI*(inverse*m)/column;
- * w.re = cos(theta);
- * w.im = sin(theta);
- */
- w.re = costab[ column + m ];
- w.im = sintab[ column + m ] * inverse;
- /* Do all pairs of nodes */
- for (node1 = &dat[m]; node1 < &dat[numpoints]; node1 += step) {
- node2 = node1 + column;
- /*
- * Want to compute:
- * dat[node2] = dat[node1] - w * dat[node2];
- * dat[node1] = dat[node1] + w * dat[node2];
- *
- * We do all this with pointers now.
- */
-
- /*cmult(&w, &dat[node2], &temp);*/
- temp.re = w.re*node2->re - w.im*node2->im;
- temp.im = w.im*node2->re + w.re*node2->im;
-
- /*csub(&dat[node1], &temp, &dat[node2]);*/
- node2->re = node1->re - temp.re;
- node2->im = node1->im - temp.im;
-
- /*cadd(&dat[node1], &temp, &dat[node1]);*/
- node1->re += temp.re;
- node1->im += temp.im;
- }
- }
- }
-
- /* Scale Data (on forward transform only) */
- /*
- * Technically speaking this gives us the periodogram. XXX
- * The canonical definition does the scaleing only
- * after the inverse xform. Our method may hurt certain
- * other forms of analysis, e.g. cepstrum.
- * **** We Now Do It The Canonical Way! ****
- */
- if (inverse > 0) {
- for (node1 = &dat[0]; node1 < &dat[numpoints]; node1++) {
- /* cdiv(&dat[i], &const, &dat[i]); */
- node1->re /= (double)numpoints;
- node1->im /= (double)numpoints;
- }
+ int i;
+ fftw_plan plan;
+
+ plan = fftw_plan_dft_1d(num, dat, dat, dir, FFTW_ESTIMATE);
+ fftw_execute(plan);
+
+ /* We need to downscale the output in case of inverse fft */
+ if (dir == FFTW_BACKWARD) {
+ for (i = 0; i < num; i++) {
+ dat[i][0] = dat[i][0] / num;
+ dat[i][1] = dat[i][1] / num;
+ }
}
+ fftw_destroy_plan(plan);
}
-/**** COMPLEX ARITHMETIC ROUTINES ****/
-/**** NO LONGER USED BY TRANSFORMS ***/
-/*
- * CADD - Complex ADDition
- */
-void
-cadd(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
-{
- result->re = val1->re + val2->re;
- result->im = val1->im + val2->im;
-}
-
-/*
- * CSUB - Complex SUBtraction
- */
-void
-csub(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
-{
- result->re = val1->re - val2->re;
- result->im = val1->im - val2->im;
-}
-
-/*
- * CMULT - Complex MULTiply
- */
-void
-cmult(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
-{
- result->re = val1->re*val2->re - val1->im*val2->im;
- result->im = val1->im*val2->re + val1->re*val2->im;
-}
-
-/*
- * CDIV - Complex DIVide
- */
void
cdiv(COMPLEX *result, COMPLEX *val1, COMPLEX *val2)
{
--
1.7.6
------------------------------------------------------------------------------
uberSVN's rich system and user administration capabilities and model
configuration take the hassle out of deploying and managing Subversion and
the tools developers use with it. Learn more about uberSVN and get a free
download at: http://p.sf.net/sfu/wandisco-dev2dev
_______________________________________________
BRL-CAD Developer mailing list
brlcad-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-devel
