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");
+}
+