spectrogram: add FFTW3 support Sox specrogram has two different algorithms to do spectrograms: lsx_safe_rdft() for dft_size of powers of two (-ysize=2^n+1) and rdft_p(), private to spectrogram.c, which does any size but is from 150 to 250 times slower. This adds FFTW3 support, which is about the same speed as the old lsx routine but works at any size: - stuff in "configure" to autodetect it and "./configure --without-fftw" to forcibly disable it - code in src/spectrogram.c to use FFTW if it is enabled or the old routines otherwise - changes to debian/control to build with FFTW (though the debian package build dies for other reasons not concerning this patch). The output from the old algorithms and FFTW (the PNGs) are identical. diff --git a/configure.ac b/configure.ac index 23138a9..4bc3064 100644 --- a/configure.ac +++ b/configure.ac @@ -332,6 +332,28 @@ AC_SUBST(PNG_LIBS) +dnl Check for FFTW3 libraries +AC_ARG_WITH(fftw, + AS_HELP_STRING([--without-fftw], + [Don't try to use FFTW])) +using_fftw=no +if test "$with_fftw" != "no"; then + AC_CHECK_HEADERS(fftw3.h, using_fftw=yes) + if test "$using_fftw" = "yes"; then + AC_CHECK_LIB(fftw3, fftw_execute, FFTW_LIBS="$FFTW_LIBS -lfftw3", using_fftw=no) + fi + if test "$with_fftw" = "yes" -a "$using_fftw" = "no"; then + AC_MSG_FAILURE([cannot find FFTW3]) + fi +fi +if test "$using_fftw" = yes; then + AC_DEFINE(HAVE_FFTW, 1, [Define to 1 if you have FFTW3.]) +fi +AM_CONDITIONAL(HAVE_FFTW, test x$using_fftw = xyes) +AC_SUBST(FFTW_LIBS) + + + dnl Test for LADSPA AC_ARG_WITH(ladspa, AS_HELP_STRING([--without-ladspa], [Don't try to use LADSPA])) @@ -756,6 +778,7 @@ echo "OTHER OPTIONS" echo "ladspa effects.............$using_ladspa" echo "magic support..............$using_magic" echo "png support................$using_png" +echo "FFTW support...............$using_fftw" if test "x$OPENMP_CFLAGS" = "x"; then echo "OpenMP support.............no" else diff --git a/debian/control b/debian/control index dc3c34b..94681e8 100644 --- a/debian/control +++ b/debian/control @@ -8,6 +8,7 @@ Build-Depends: dh-autoreconf, ladspa-sdk, libao-dev, libasound2-dev [linux-any], + libfftw3-dev, libgsm1-dev, libid3tag0-dev, libltdl3-dev, diff --git a/src/Makefile.am b/src/Makefile.am index 7cceaaf..5838f1f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -90,6 +90,9 @@ libsox_la_LIBADD = @PNG_LIBS@ if HAVE_MAGIC libsox_la_LIBADD += @MAGIC_LIBS@ endif +if HAVE_FFTW +libsox_la_LIBADD += @FFTW_LIBS@ +endif libsox_la_LIBADD += @GOMP_LIBS@ diff --git a/src/spectrogram.c b/src/spectrogram.c index 34457f2..3d8c208 100644 --- a/src/spectrogram.c +++ b/src/spectrogram.c @@ -30,6 +30,10 @@ #endif #include +#ifdef HAVE_FFTW3_H +#include +#endif + /* For SET_BINARY_MODE: */ #include #ifdef HAVE_IO_H @@ -80,6 +84,9 @@ typedef struct { double block_norm, max; double * magnitudes; /* [(dft_size / 2) + 1] */ float * dBfs; +#if HAVE_FFTW + fftw_plan fftw_plan; /* Used if FFT_type == FFT_fftw */ +#endif } priv_t; #define secs(cols) \ @@ -197,6 +204,8 @@ static double make_window(priv_t * p, int end) return sum; } +#if !HAVE_FFTW + static double * rdft_init(size_t n) { double * q = lsx_malloc(2 * (n / 2 + 1) * n * sizeof(*q)), * p = q; @@ -218,6 +227,8 @@ static void rdft_p(double const * q, double const * in, double * out, int n) } } +#endif /* HAVE_FFTW */ + static int start(sox_effect_t * effp) { priv_t * p = (priv_t *)effp->priv; @@ -266,8 +277,10 @@ static int start(sox_effect_t * effp) if (p->y_size) { p->dft_size = 2 * (p->y_size - 1); +#if !HAVE_FFTW if (!is_p2(p->dft_size) && !effp->flow) p->shared = rdft_init((size_t)(p->dft_size)); +#endif } else { int y = max(32, (p->Y_size? p->Y_size : 550) / effp->in_signal.channels - 2); for (p->dft_size = 128; p->dft_size <= y; p->dft_size <<= 1); @@ -279,12 +292,20 @@ static int start(sox_effect_t * effp) p->window = lsx_calloc(p->dft_size + 1, sizeof(*(p->window))); p->magnitudes = lsx_calloc((p->dft_size / 2) + 1, sizeof(*(p->magnitudes))); + /* Initialize the FFT routine */ +#if HAVE_FFTW + /* We have one FFT plan per flow because the input/output arrays differ. */ + p->fftw_plan = fftw_plan_r2r_1d(p->dft_size, p->dft_buf, p->dft_buf, + FFTW_R2HC, FFTW_MEASURE); +#else if (is_p2(p->dft_size) && !effp->flow) lsx_safe_rdft(p->dft_size, 1, p->dft_buf); +#endif + lsx_debug("duration=%g x_size=%i pixels_per_sec=%g dft_size=%i", duration, p->x_size, pixels_per_sec, p->dft_size); p->end = p->dft_size; - p->rows = (p->dft_size >> 1) + 1; + p->rows = (p->dft_size / 2) + 1; actual = make_window(p, p->last_end = 0); lsx_debug("window_density=%g", actual / p->dft_size); p->step_size = (p->slack_overlap? sqrt(actual * p->dft_size) : actual) + .5; @@ -359,6 +380,19 @@ static int flow(sox_effect_t * effp, if ((p->end = max(p->end, p->end_min)) != p->last_end) make_window(p, p->last_end = p->end); for (i = 0; i < p->dft_size; ++i) p->dft_buf[i] = p->buf[i] * p->window[i]; +#if HAVE_FFTW + fftw_execute(p->fftw_plan); + /* Convert from FFTW's "half complex" format to an array of magnitudes. + * In HC format, the values are stored: + * r0, r1, r2 ... r(n/2), i(n+1)/2-1 .. i2, i1 + */ + p->magnitudes[0] += sqr(p->dft_buf[0]); + for (i = 1; i < p->dft_size / 2; ++i) { + p->magnitudes[i] += sqr(p->dft_buf[i]) + sqr(p->dft_buf[p->dft_size - i +]); + } + p->magnitudes[p->dft_size / 2] += sqr(p->dft_buf[p->dft_size / 2]); +#else if (is_p2(p->dft_size)) { lsx_safe_rdft(p->dft_size, 1, p->dft_buf); p->magnitudes[0] += sqr(p->dft_buf[0]); @@ -367,6 +401,8 @@ static int flow(sox_effect_t * effp, p->magnitudes[p->dft_size >> 1] += sqr(p->dft_buf[1]); } else rdft_p(*p->shared_ptr, p->dft_buf, p->magnitudes, p->dft_size); +#endif + if (++p->block_num == p->block_steps && do_column(effp) == SOX_EOF) return SOX_EOF; }