diff options
Diffstat (limited to 'pd/src/d_fft_fftw.c')
-rw-r--r-- | pd/src/d_fft_fftw.c | 138 |
1 files changed, 104 insertions, 34 deletions
diff --git a/pd/src/d_fft_fftw.c b/pd/src/d_fft_fftw.c index 0ead62be..04e4729d 100644 --- a/pd/src/d_fft_fftw.c +++ b/pd/src/d_fft_fftw.c @@ -4,31 +4,32 @@ /* --------- 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. +/* changes and additions for FFTW3 by Thomas Grill */ -#include <fftw.h> +#include "m_pd.h" +#include <fftw3.h> int ilog2(int n); -#define MINFFT 5 +#define MINFFT 0 #define MAXFFT 30 /* from the FFTW website: - fftw_complex in[N], out[N]; + #include <fftw3.h> + ... + { + fftw_complex *in, *out; fftw_plan p; ... - p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE); + in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); + out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); + p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); ... - fftw_one(p, in, out); + fftw_execute(p); ... fftw_destroy_plan(p); + fftw_free(in); fftw_free(out); + } 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 @@ -37,57 +38,126 @@ respectively. The flags argument is either FFTW_MEASURE */ -static fftw_plan fftw_pvec[2 * (MAXFFT+1 - MINFFT)]; +/* complex stuff */ + +typedef struct { + fftwf_plan plan; + fftwf_complex *in,*out; +} cfftw_info; + +static cfftw_info cfftw_fwd[MAXFFT+1 - MINFFT],cfftw_bwd[MAXFFT+1 - MINFFT]; + +static cfftw_info *cfftw_getplan(int n,int fwd) +{ + cfftw_info *info; + int logn = ilog2(n); + if (logn < MINFFT || logn > MAXFFT) + return (0); + info = (fwd?cfftw_fwd:cfftw_bwd)+(logn-MINFFT); + if (!info->plan) + { + info->in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * n); + info->out = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * n); + info->plan = fftwf_plan_dft_1d(n, info->in, info->out, fwd?FFTW_FORWARD:FFTW_BACKWARD, FFTW_MEASURE); + } + return info; +} + -static fftw_plan fftw_getplan(int n, int dir) +/* real stuff */ + +typedef struct { + fftwf_plan plan; + float *in,*out; +} rfftw_info; + +static rfftw_info rfftw_fwd[MAXFFT+1 - MINFFT],rfftw_bwd[MAXFFT+1 - MINFFT]; + +static rfftw_info *rfftw_getplan(int n,int fwd) { - logn = ilog2(n); + rfftw_info *info; + int 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]); + info = (fwd?rfftw_fwd:rfftw_bwd)+(logn-MINFFT); + if (!info->plan) + { + info->in = (float*) fftwf_malloc(sizeof(float) * n); + info->out = (float*) fftwf_malloc(sizeof(float) * n); + info->plan = fftwf_plan_r2r_1d(n, info->in, info->out, fwd?FFTW_R2HC:FFTW_HC2R, FFTW_MEASURE); + } + return info; } + + 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) +static void mayer_do_cfft(int n, float *fz1, float *fz2, int fwd) { - 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; + float *fz; + cfftw_info *p = cfftw_getplan(n, fwd); 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; + + for (i = 0, fz = (float *)p->in; i < n; i++) + fz[i*2] = fz1[i], fz[i*2+1] = fz2[i]; + + fftwf_execute(p->plan); + + for (i = 0, fz = (float *)p->out; i < n; i++) + fz1[i] = fz[i*2], fz2[i] = fz[i*2+1]; } EXTERN void mayer_fft(int n, float *fz1, float *fz2) { - mayer_dofft(n, fz1, fz2, FFTW_FORWARD); + mayer_do_cfft(n, fz1, fz2, 1); } EXTERN void mayer_ifft(int n, float *fz1, float *fz2) { - mayer_dofft(n, fz1, fz2, FFTW_BACKWARD); + mayer_do_cfft(n, fz1, fz2, 0); } +/* + in the following the sign flips are done to + be compatible with the mayer_fft implementation, + but it's probably the mayer_fft that should be corrected... +*/ + EXTERN void mayer_realfft(int n, float *fz) { - post("rfft: not yet implemented"); + int i; + rfftw_info *p = rfftw_getplan(n, 1); + if (!p) + return; + + for (i = 0; i < n; i++) + p->in[i] = fz[i]; + fftwf_execute(p->plan); + for (i = 0; i < n/2+1; i++) + fz[i] = p->out[i]; + for (; i < n; i++) + fz[i] = -p->out[i]; } EXTERN void mayer_realifft(int n, float *fz) { - post("rifft: not yet implemented"); + int i; + rfftw_info *p = rfftw_getplan(n, 0); + if (!p) + return; + + for (i = 0; i < n/2+1; i++) + p->in[i] = fz[i]; + for (; i < n; i++) + p->in[i] = -fz[i]; + fftwf_execute(p->plan); + for (i = 0; i < n; i++) + fz[i] = p->out[i]; } |