From b9d6f6cdcc66710f3b321f6a08159f186d71684d Mon Sep 17 00:00:00 2001
From: Ivica Ico Bukvic <ico@vt.edu>
Date: Tue, 19 Feb 2013 18:20:27 -0500
Subject: [PATCH] improved fftw3 support
---
pd/src/configure.in | 38 ++++++------
pd/src/d_fft_fftw.c | 138 +++++++++++++++++++++++++++++++++-----------
2 files changed, 125 insertions(+), 51 deletions(-)
diff --git a/pd/src/configure.in b/pd/src/configure.in
index ca68c45ff..b378e0761 100644
--- a/pd/src/configure.in
+++ b/pd/src/configure.in
@@ -47,8 +47,10 @@ AC_ARG_ENABLE(setuid, [ --enable-setuid install as setuid (linux)],
setuid=$enableval)
AC_ARG_ENABLE(fftw, [ --enable-fftw use FFTW package],
fftw=$enableval)
-AC_ARG_ENABLE(fat, [ --disable-fat build fat binary on Mac OS X],
+AC_ARG_ENABLE(fat, [ --disable-fat build fat binary on Mac OS X],
fat=$enableval, fat="yes")
+AC_ARG_ENABLE(tk, [ --disable-tk build without tcl/tk-GUI],
+ tk=$enableval)
dnl Checks for programs.
AC_PROG_CC
@@ -96,8 +98,8 @@ AC_CHECK_LIB(pthread, pthread_create,PDLIB="$PDLIB -lpthread",
dnl Check for fftw package
if test x$fftw = "xyes";
then
-AC_CHECK_LIB(fftw, fftw_one,PDLIB="$PDLIB -lfftw",
- echo "fftw package not found - using built-in FFT"; fftw=no)
+AC_CHECK_LIB(fftw3f, fftwf_version,PDLIB="$PDLIB -lfftw3f",
+ echo "fftw3 package not found - using built-in FFT"; fftw=no)
fi
dnl look for tcl 8.x... do I really have to go through all this!?
@@ -176,13 +178,10 @@ else
GUISRC=
fi
-# leave the executable as 'pd' on Mac OS X and Windows, but call it
-# 'pdl2ork' on GNU/Linux
-PDEXEC=pd
+PDEXEC=pd-l2ork
-if test `uname -s` = Linux;
+if test `uname -s` = Linux -o `uname -s` = "GNU/kFreeBSD" -o `uname -s` = "GNU";
then
- PDEXEC=pd-l2ork
dnl Ckecking for ALSA
echo .................... alsa= $alsa
dnl This should be fixed so Pd can use ALSA shared libraries where appropriate.
@@ -206,7 +205,16 @@ dnl This should be fixed so Pd can use ALSA shared libraries where appropriate.
CPPFLAGS="-DDL_OPEN -DPA_USE_OSS -DUNIX -DUNISTD\
-DUSEAPI_OSS \
-fno-strict-aliasing"
- SYSSRC="s_midi_oss.c s_audio_oss.c"
+dnl No OSS on hurd
+ if test `uname -s` = "GNU";
+ then
+ SYSSRC="s_midi_dummy.c"
+ else
+ SYSSRC="s_midi_oss.c s_audio_oss.c"
+ CPPFLAGS=$CPPFLAGS" -DPA_USE_OSS -DUSEAPI_OSS"
+ fi
+
+
if test x$alsa = "xyes";
then
SYSSRC=$SYSSRC" s_audio_alsa.c s_audio_alsamm.c s_midi_alsa.c"
@@ -223,7 +231,6 @@ dnl This should be fixed so Pd can use ALSA shared libraries where appropriate.
-I../portmidi/pm_common \
-I../portmidi/pm_linux"
SYSSRC="s_audio_pa.c \
- s_audio_pablio.c \
s_audio_paring.c \
../portaudio/src/common/pa_allocation.c \
../portaudio/src/common/pa_converters.c \
@@ -281,7 +288,7 @@ then
-framework AudioUnit -framework AudioToolbox \
-framework Carbon -framework CoreMIDI"
EXT=pd_darwin
- CPPFLAGS="-DDL_OPEN -DMACOSX -DUNISTD -I/usr/X11R6/include \
+ CPPFLAGS="-DHAVE_LIBDL -DMACOSX -DHAVE_UNISTD_H -I/usr/X11R6/include \
-I../portaudio/include -I../portaudio/src/common \
-I../portaudio/src/os/mac_osx/ \
-I../portmidi/pm_common -I../portmidi/pm_mac \
@@ -303,7 +310,6 @@ then
EXTERNTARGET=d_ppc
fi
SYSSRC="s_midi_pm.c s_audio_pa.c \
- s_audio_pablio.c \
s_audio_paring.c \
../portaudio/src/common/pa_allocation.c \
../portaudio/src/common/pa_converters.c \
@@ -323,8 +329,6 @@ then
../portaudio/src/hostapi/coreaudio/pa_mac_core_utilities.c \
../portmidi/pm_mac/pmmac.c \
../portmidi/pm_mac/pmmacosxcm.c \
- ../portmidi/pm_mac/finddefault.c \
- ../portmidi/pm_mac/readbinaryplist.c \
../portmidi/pm_common/pmutil.c \
../portmidi/pm_common/portmidi.c \
../portmidi/porttime/ptmacosx_cf.c "
@@ -357,7 +361,7 @@ then
OSNUMBER=2
if test x$jack = "xyes";
then
- LDFLAGS=$LDFLAGS" -weak_framework Jackmp"
+ LDFLAGS=$LDFLAGS" -weak_framework Jackmp"
fi
if test x$jack = "xrun";
then
@@ -377,7 +381,7 @@ then
-mwindows -mms-bitfields "$MORECFLAGS
PDLIB=$PDLIB" -lwsock32 -lwinmm -lole32 -lstdc++"
- SYSSRC="s_audio_pa.c s_audio_pablio.c s_audio_paring.c \
+ SYSSRC="s_audio_pa.c s_audio_paring.c \
s_audio_mmio.c s_midi_mmio.c \
../portaudio/src/common/pa_allocation.c \
../portaudio/src/common/pa_converters.c \
@@ -418,7 +422,7 @@ fi
if test x$fftw = "xyes";
then
SYSSRC=$SYSSRC" d_fft_fftw.c d_fftroutine.c"
- LDFLAGS=$LDFLAGS" -lfftw"
+ LDFLAGS=$LDFLAGS" -lfftw3"
else
SYSSRC=$SYSSRC" d_fft_mayer.c d_fftroutine.c"
fi
diff --git a/pd/src/d_fft_fftw.c b/pd/src/d_fft_fftw.c
index 0ead62be4..04e4729de 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];
}
--
GitLab