Changeset 695e171 for src/spectral/fft.c


Ignore:
Timestamp:
Sep 6, 2015, 10:40:14 AM (9 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, sampler
Children:
3d30b90
Parents:
65c352e (diff), 827267b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'develop' into awhitening

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/spectral/fft.c

    r65c352e r695e171  
    6767#warning "Using aubio in double precision with fftw3 in single precision"
    6868#endif /* HAVE_AUBIO_DOUBLE */
    69 #define real_t float 
     69#define real_t float
    7070#else /* HAVE_FFTW3F */
    7171#if !HAVE_AUBIO_DOUBLE
    7272#warning "Using aubio in single precision with fftw3 in double precision"
    7373#endif /* HAVE_AUBIO_DOUBLE */
    74 #define real_t double 
     74#define real_t double
    7575#endif /* HAVE_FFTW3F */
    7676
    … …  
    8282// https://developer.apple.com/library/mac/#documentation/Accelerate/Reference/vDSPRef/Reference/reference.html
    8383#include <Accelerate/Accelerate.h>
     84
     85#if !HAVE_AUBIO_DOUBLE
     86#define aubio_vDSP_ctoz                vDSP_ctoz
     87#define aubio_vDSP_fft_zrip            vDSP_fft_zrip
     88#define aubio_vDSP_ztoc                vDSP_ztoc
     89#define aubio_vDSP_zvmags              vDSP_zvmags
     90#define aubio_vDSP_zvphas              vDSP_zvphas
     91#define aubio_vDSP_vsadd               vDSP_vsadd
     92#define aubio_vDSP_vsmul               vDSP_vsmul
     93#define aubio_vDSP_create_fftsetup     vDSP_create_fftsetup
     94#define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetup
     95#define aubio_DSPComplex               DSPComplex
     96#define aubio_DSPSplitComplex          DSPSplitComplex
     97#define aubio_FFTSetup                 FFTSetup
     98#define aubio_vvsqrt                   vvsqrtf
     99#else
     100#define aubio_vDSP_ctoz                vDSP_ctozD
     101#define aubio_vDSP_fft_zrip            vDSP_fft_zripD
     102#define aubio_vDSP_ztoc                vDSP_ztocD
     103#define aubio_vDSP_zvmags              vDSP_zvmagsD
     104#define aubio_vDSP_zvphas              vDSP_zvphasD
     105#define aubio_vDSP_vsadd               vDSP_vsaddD
     106#define aubio_vDSP_vsmul               vDSP_vsmulD
     107#define aubio_vDSP_create_fftsetup     vDSP_create_fftsetupD
     108#define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetupD
     109#define aubio_DSPComplex               DSPDoubleComplex
     110#define aubio_DSPSplitComplex          DSPDoubleSplitComplex
     111#define aubio_FFTSetup                 FFTSetupD
     112#define aubio_vvsqrt                   vvsqrt
     113#endif /* HAVE_AUBIO_DOUBLE */
    84114
    85115#else                         // using OOURA
    … …  
    96126  real_t *in, *out;
    97127  fftw_plan pfw, pbw;
    98   fft_data_t * specdata;     /* complex spectral data */
     128  fft_data_t * specdata;      /* complex spectral data */
    99129#else
    100130#ifdef HAVE_ACCELERATE        // using ACCELERATE
    101131  int log2fftsize;
    102 #if !HAVE_AUBIO_DOUBLE
    103   FFTSetup fftSetup;
    104   DSPSplitComplex spec;
    105   float *in, *out;
    106 #else
    107   FFTSetupD fftSetup;
    108   DSPDoubleSplitComplex spec;
    109   double *in, *out;
    110 #endif
     132  aubio_FFTSetup fftSetup;
     133  aubio_DSPSplitComplex spec;
     134  smpl_t *in, *out;
    111135#else                         // using OOURA
    112136  smpl_t *in, *out;
    … …  
    120144aubio_fft_t * new_aubio_fft (uint_t winsize) {
    121145  aubio_fft_t * s = AUBIO_NEW(aubio_fft_t);
     146  if ((sint_t)winsize < 1) {
     147    AUBIO_ERR("fft: got winsize %d, but can not be < 1\n", winsize);
     148    goto beach;
     149  }
    122150#ifdef HAVE_FFTW3
    123151  uint_t i;
    … …  
    154182  s->compspec = new_fvec(winsize);
    155183  s->log2fftsize = (uint_t)log2f(s->fft_size);
    156 #if !HAVE_AUBIO_DOUBLE
    157   s->in = AUBIO_ARRAY(float, s->fft_size);
    158   s->out = AUBIO_ARRAY(float, s->fft_size);
    159   s->spec.realp = AUBIO_ARRAY(float, s->fft_size/2);
    160   s->spec.imagp = AUBIO_ARRAY(float, s->fft_size/2);
    161   s->fftSetup = vDSP_create_fftsetup(s->log2fftsize, FFT_RADIX2);
    162 #else
    163   s->in = AUBIO_ARRAY(double, s->fft_size);
    164   s->out = AUBIO_ARRAY(double, s->fft_size);
    165   s->spec.realp = AUBIO_ARRAY(double, s->fft_size/2);
    166   s->spec.imagp = AUBIO_ARRAY(double, s->fft_size/2);
    167   s->fftSetup = vDSP_create_fftsetupD(s->log2fftsize, FFT_RADIX2);
    168 #endif
    169 #else                         // using OOURA
     184  s->in = AUBIO_ARRAY(smpl_t, s->fft_size);
     185  s->out = AUBIO_ARRAY(smpl_t, s->fft_size);
     186  s->spec.realp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
     187  s->spec.imagp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
     188  s->fftSetup = aubio_vDSP_create_fftsetup(s->log2fftsize, FFT_RADIX2);
     189#else                         // using OOURA
     190  if (aubio_is_power_of_two(winsize) != 1) {
     191    AUBIO_ERR("fft: can only create with sizes power of two,"
     192              " requested %d\n", winsize);
     193    goto beach;
     194  }
    170195  s->winsize = winsize;
    171196  s->fft_size = winsize / 2 + 1;
    … …  
    179204#endif /* HAVE_FFTW3 */
    180205  return s;
     206beach:
     207  AUBIO_FREE(s);
     208  return NULL;
    181209}
    182210
    … …  
    192220  AUBIO_FREE(s->spec.realp);
    193221  AUBIO_FREE(s->spec.imagp);
    194 #if !HAVE_AUBIO_DOUBLE
    195   vDSP_destroy_fftsetup(s->fftSetup);
    196 #else
    197   vDSP_destroy_fftsetupD(s->fftSetup);
    198 #endif
     222  aubio_vDSP_destroy_fftsetup(s->fftSetup);
    199223#else                         // using OOURA
    200224  AUBIO_FREE(s->w);
    … …  
    219243void aubio_fft_do_complex(aubio_fft_t * s, fvec_t * input, fvec_t * compspec) {
    220244  uint_t i;
     245#ifndef HAVE_MEMCPY_HACKS
    221246  for (i=0; i < s->winsize; i++) {
    222247    s->in[i] = input->data[i];
    223248  }
     249#else
     250  memcpy(s->in, input->data, s->winsize * sizeof(smpl_t));
     251#endif /* HAVE_MEMCPY_HACKS */
    224252#ifdef HAVE_FFTW3             // using FFTW3
    225253  fftw_execute(s->pfw);
    … …  
    238266#else /* HAVE_FFTW3 */
    239267#ifdef HAVE_ACCELERATE        // using ACCELERATE
    240 #if !HAVE_AUBIO_DOUBLE
    241268  // convert real data to even/odd format used in vDSP
    242   vDSP_ctoz((DSPComplex*)s->in, 2, &s->spec, 1, s->fft_size/2);
     269  aubio_vDSP_ctoz((aubio_DSPComplex*)s->in, 2, &s->spec, 1, s->fft_size/2);
    243270  // compute the FFT
    244   vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
    245 #else
    246   // convert real data to even/odd format used in vDSP
    247   vDSP_ctozD((DSPDoubleComplex*)s->in, 2, &s->spec, 1, s->fft_size/2);
    248   // compute the FFT
    249   vDSP_fft_zripD(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
    250 #endif
     271  aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
    251272  // convert from vDSP complex split to [ r0, r1, ..., rN, iN-1, .., i2, i1]
    252273  compspec->data[0] = s->spec.realp[0];
    … …  
    258279  // apply scaling
    259280  smpl_t scale = 1./2.;
    260 #if !HAVE_AUBIO_DOUBLE
    261   vDSP_vsmul(compspec->data, 1, &scale, compspec->data, 1, s->fft_size);
    262 #else
    263   vDSP_vsmulD(compspec->data, 1, &scale, compspec->data, 1, s->fft_size);
    264 #endif
     281  aubio_vDSP_vsmul(compspec->data, 1, &scale, compspec->data, 1, s->fft_size);
    265282#else                         // using OOURA
    266283  rdft(s->winsize, 1, s->in, s->ip, s->w);
    … …  
    305322    s->out[2 * i + 1] = compspec->data[s->winsize - i];
    306323  }
    307 #if !HAVE_AUBIO_DOUBLE
    308324  // convert to split complex format used in vDSP
    309   vDSP_ctoz((DSPComplex*)s->out, 2, &s->spec, 1, s->fft_size/2);
     325  aubio_vDSP_ctoz((aubio_DSPComplex*)s->out, 2, &s->spec, 1, s->fft_size/2);
    310326  // compute the FFT
    311   vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
     327  aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
    312328  // convert result to real output
    313   vDSP_ztoc(&s->spec, 1, (DSPComplex*)output->data, 2, s->fft_size/2);
     329  aubio_vDSP_ztoc(&s->spec, 1, (aubio_DSPComplex*)output->data, 2, s->fft_size/2);
    314330  // apply scaling
    315331  smpl_t scale = 1.0 / s->winsize;
    316   vDSP_vsmul(output->data, 1, &scale, output->data, 1, s->fft_size);
    317 #else
    318   // convert to split complex format used in vDSP
    319   vDSP_ctozD((DSPDoubleComplex*)s->out, 2, &s->spec, 1, s->fft_size/2);
    320   // compute the FFT
    321   vDSP_fft_zripD(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
    322   // convert result to real output
    323   vDSP_ztocD(&s->spec, 1, (DSPDoubleComplex*)output->data, 2, s->fft_size/2);
    324   // apply scaling
    325   smpl_t scale = 1.0 / s->winsize;
    326   vDSP_vsmulD(output->data, 1, &scale, output->data, 1, s->fft_size);
    327 #endif
     332  aubio_vDSP_vsmul(output->data, 1, &scale, output->data, 1, s->fft_size);
    328333#else                         // using OOURA
    329334  smpl_t scale = 2.0 / s->winsize;
    … …  
    392397  uint_t i;
    393398  for (i = 0; i < compspec->length / 2 + 1; i++) {
    394     compspec->data[i] = 
     399    compspec->data[i] =
    395400      spectrum->norm[i]*COS(spectrum->phas[i]);
    396401  }
Note: See TracChangeset for help on using the changeset viewer.