Changeset 54dd945 for src/spectral

Timestamp:

Apr 11, 2013, 2:31:40 AM (12 years ago)

Author:

Paul Brossier <piem@piem.org>

Branches:

feature/autosink, feature/cnn, feature/cnn_org, feature/constantq, feature/crepe, feature/crepe_org, feature/pitchshift, feature/pydocstrings, feature/timestretch, fix/ffmpeg5, master, pitchshift, sampler, timestretch, yinfft+

Children:

ac67de7

Parents:

e8bc8e9

Message:

src/spectral/fft.c: add vDSP Accelerate

File:

• src/spectral/fft.c (modified) (11 diffs)

src/spectral/fft.c

@@ -25 +25 @@
 #include "spectral/fft.h"
 
-#ifdef HAVE_FFTW3
+#ifdef HAVE_FFTW3             // using FFTW3
 /* note that <complex.h> is not included here but only in aubio_priv.h, so that
  * c++ projects can still use their own complex definition. */
@@ -78 +78 @@
 pthread_mutex_t aubio_fftw_mutex = PTHREAD_MUTEX_INITIALIZER;
 
+#else
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+// https://developer.apple.com/library/mac/#documentation/Accelerate/Reference/vDSPRef/Reference/reference.html
+#include <Accelerate/Accelerate.h>
+
+#else                         // using OOURA
+// let's use ooura instead
+extern void rdft(int, int, double *, int *, double *);
+
+#endif /* HAVE_ACCELERATE */
 #endif /* HAVE_FFTW3 */
 
@@ -83 +93 @@
   uint_t winsize;
   uint_t fft_size;
-#ifdef HAVE_FFTW3
+#ifdef HAVE_FFTW3             // using FFTW3
   real_t *in, *out;
   fftw_plan pfw, pbw;
   fft_data_t * specdata;     /* complex spectral data */
 #else
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+  int log2fftsize;
+  FFTSetup fftSetup;
+  DSPSplitComplex spec;
+  float *in, *out;
+#else                         // using OOURA
   double *in, *out;
   double *w;
   int *ip;
+#endif /* HAVE_ACCELERATE */
 #endif /* HAVE_FFTW3 */
   fvec_t * compspec;
 };
-
-#ifndef HAVE_FFTW3
-// let's use ooura instead
-extern void rdft(int, int, double *, int *, double *);
-#endif
 
 aubio_fft_t * new_aubio_fft (uint_t winsize) {
@@ -131 +143 @@
   }
 #else
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+  s->winsize = winsize;
+  s->fft_size = winsize;
+  s->compspec = new_fvec(winsize);
+  s->log2fftsize = (uint_t)log2f(s->fft_size);
+  s->in = AUBIO_ARRAY(float, s->fft_size);
+  s->out = AUBIO_ARRAY(float, s->fft_size);
+  s->spec.realp = AUBIO_ARRAY(float, s->fft_size/2);
+  s->spec.imagp = AUBIO_ARRAY(float, s->fft_size/2);
+  s->fftSetup = vDSP_create_fftsetup(s->log2fftsize, FFT_RADIX2);
+#else                         // using OOURA
   s->winsize = winsize;
   s->fft_size = winsize / 2 + 1;
@@ -139 +162 @@
   s->w     = AUBIO_ARRAY(double, s->fft_size);
   s->ip[0] = 0;
-#endif
+#endif /* HAVE_ACCELERATE */
+#endif /* HAVE_FFTW3 */
   return s;
 }
@@ -146 +170 @@
   /* destroy data */
   del_fvec(s->compspec);
-#ifdef HAVE_FFTW3
+#ifdef HAVE_FFTW3             // using FFTW3
   fftw_destroy_plan(s->pfw);
   fftw_destroy_plan(s->pbw);
   fftw_free(s->specdata);
 #else /* HAVE_FFTW3 */
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+  AUBIO_FREE(s->spec.realp);
+  AUBIO_FREE(s->spec.imagp);
+#else                         // using OOURA
   AUBIO_FREE(s->w);
   AUBIO_FREE(s->ip);
+#endif /* HAVE_ACCELERATE */
 #endif /* HAVE_FFTW3 */
   AUBIO_FREE(s->out);
@@ -174 +203 @@
     s->in[i] = input->data[i];
   }
-#ifdef HAVE_FFTW3
+#ifdef HAVE_FFTW3             // using FFTW3
   fftw_execute(s->pfw);
 #ifdef HAVE_COMPLEX_H
@@ -189 +218 @@
 #endif /* HAVE_COMPLEX_H */
 #else /* HAVE_FFTW3 */
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+  // convert real data to even/odd format used in vDSP
+  vDSP_ctoz((COMPLEX*)s->in, 2, &s->spec, 1, s->fft_size/2);
+  // compute the FFT
+  vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
+  // convert from vDSP complex split to [ r0, r1, ..., rN, iN-1, .., i2, i1]
+  compspec->data[0] = s->spec.realp[0];
+  compspec->data[s->fft_size / 2] = s->spec.imagp[0];
+  for (i = 1; i < s->fft_size / 2; i++) {
+    compspec->data[i] = s->spec.realp[i];
+    compspec->data[s->fft_size - i] = s->spec.imagp[i];
+  }
+  // apply scaling
+  smpl_t scale = 1./2.;
+  vDSP_vsmul(compspec->data, 1, &scale, compspec->data, 1, s->fft_size);
+#else                         // using OOURA
   rdft(s->winsize, 1, s->in, s->ip, s->w);
   compspec->data[0] = s->in[0];
@@ -196 +241 @@
     compspec->data[s->winsize - i] = - s->in[2 * i + 1];
   }
+#endif /* HAVE_ACCELERATE */
 #endif /* HAVE_FFTW3 */
 }
@@ -220 +266 @@
   }
 #else /* HAVE_FFTW3 */
+#ifdef HAVE_ACCELERATE        // using ACCELERATE
+  // convert from real imag  [ r0, r1, ..., rN, iN-1, .., i2, i1]
+  // to vDSP packed format   [ r0, rN, r1, i1, ..., rN-1, iN-1 ]
+  s->out[0] = compspec->data[0];
+  s->out[1] = compspec->data[s->winsize / 2];
+  for (i = 1; i < s->fft_size / 2; i++) {
+    s->out[2 * i] = compspec->data[i];
+    s->out[2 * i + 1] = compspec->data[s->winsize - i];
+  }
+  // convert to split complex format used in vDSP
+  vDSP_ctoz((COMPLEX*)s->out, 2, &s->spec, 1, s->fft_size/2);
+  // compute the FFT
+  vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
+  // convert result to real output
+  vDSP_ztoc(&s->spec, 1, (COMPLEX*)output->data, 2, s->fft_size/2);
+  // apply scaling
+  smpl_t scale = 1.0 / s->winsize;
+  vDSP_vsmul(output->data, 1, &scale, output->data, 1, s->fft_size);
+#else                         // using OOURA
   smpl_t scale = 2.0 / s->winsize;
   s->out[0] = compspec->data[0];
@@ -231 +296 @@
     output->data[i] = s->out[i] * scale;
   }
+#endif /* HAVE_ACCELERATE */
 #endif /* HAVE_FFTW3 */
 }