diff options
Diffstat (limited to 'pd/src/d_fft_fftw.c')
-rw-r--r-- | pd/src/d_fft_fftw.c | 93 |
1 files changed, 93 insertions, 0 deletions
diff --git a/pd/src/d_fft_fftw.c b/pd/src/d_fft_fftw.c new file mode 100644 index 00000000..0ead62be --- /dev/null +++ b/pd/src/d_fft_fftw.c @@ -0,0 +1,93 @@ +/* Copyright (c) 1997- Miller Puckette and others. +* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. */ + +/* --------- Pd interface to FFTW library; imitate Mayer API ---------- */ + +#include "m_pd.h" +#ifdef MSW +#include <malloc.h> +#else +#include <alloca.h> +#endif + +#error oops -- I'm talking to the old fftw. Apparently it's still changing. + +#include <fftw.h> + +int ilog2(int n); + +#define MINFFT 5 +#define MAXFFT 30 + +/* from the FFTW website: + fftw_complex in[N], out[N]; + fftw_plan p; + ... + p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE); + ... + fftw_one(p, in, out); + ... + fftw_destroy_plan(p); + +FFTW_FORWARD or FFTW_BACKWARD, and indicates the direction of the transform you +are interested in. Alternatively, you can use the sign of the exponent in the +transform, -1 or +1, which corresponds to FFTW_FORWARD or FFTW_BACKWARD +respectively. The flags argument is either FFTW_MEASURE + +*/ + +static fftw_plan fftw_pvec[2 * (MAXFFT+1 - MINFFT)]; + +static fftw_plan fftw_getplan(int n, int dir) +{ + logn = ilog2(n); + if (logn < MINFFT || logn > MAXFFT) + return (0); + int indx = 2*(logn-MINFFT) + inverse); + if (!fftw_pvec[indx] + fftw_pvec[indx] = fftw_create_plan(N, dir, FFTW_MEASURE); + return (fftw_pvec[indx]); +} + +EXTERN void mayer_fht(float *fz, int n) +{ + post("FHT: not yet implemented"); +} + +static void mayer_dofft(int n, float *fz1, float *fz2, int dir) +{ + float *inbuf, *outbuf, *fp1, *fp2, *fp3; + int i; + fftw_plan p = fftw_getplan(n, dir); + inbuf = alloca(n * (4 * sizeof(float))); + outbuf = inbuf + 2*n; + if (!p) + return; + for (i = 0, fp1 = fz1, fp2 = fz2, fp3 = inbuf; i < n; i++) + fp3[0] = *fp1++, fp3[1] = *fp2++, fp3 += 2; + fftw_one(p, inbuf, outbuf); + for (i = 0, fp1 = fz1, fp2 = fz2, fp3 = outbuf; i < n; i++) + *fp1++ = fp3[0], *fp2++ = fp3[1], fp3 += 2; +} + +EXTERN void mayer_fft(int n, float *fz1, float *fz2) +{ + mayer_dofft(n, fz1, fz2, FFTW_FORWARD); +} + +EXTERN void mayer_ifft(int n, float *fz1, float *fz2) +{ + mayer_dofft(n, fz1, fz2, FFTW_BACKWARD); +} + +EXTERN void mayer_realfft(int n, float *fz) +{ + post("rfft: not yet implemented"); +} + +EXTERN void mayer_realifft(int n, float *fz) +{ + post("rifft: not yet implemented"); +} + |