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];
 }