Changeset 01d4d19


Ignore:
Timestamp:
Nov 15, 2018, 3:07:48 AM (12 months ago)
Author:
Paul Brossier <piem@piem.org>
Branches:
feature/autosink, feature/constantq, feature/pitchshift, feature/pydocstrings, feature/timestretch, master
Children:
5999ff2
Parents:
9ef3c6e
Message:

Merge branch 'fix/oddfft' (closes #207)

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • python/tests/test_fft.py

    r9ef3c6e r01d4d19  
    143143        assert_almost_equal ( r[1:], 0)
    144144
     145class aubio_fft_odd_sizes(TestCase):
     146
     147    def test_reconstruct_with_odd_size(self):
     148        win_s = 29
     149        self.recontruct(win_s, 'odd sizes not supported')
     150
     151    def test_reconstruct_with_radix15(self):
     152        win_s = 2 ** 4 * 15
     153        self.recontruct(win_s, 'radix 15 supported')
     154
     155    def test_reconstruct_with_radix5(self):
     156        win_s = 2 ** 4 * 5
     157        self.recontruct(win_s, 'radix 5 supported')
     158
     159    def test_reconstruct_with_radix3(self):
     160        win_s = 2 ** 4 * 3
     161        self.recontruct(win_s, 'radix 3 supported')
     162
     163    def recontruct(self, win_s, skipMessage):
     164        try:
     165            f = fft(win_s)
     166        except RuntimeError:
     167            self.skipTest(skipMessage)
     168        input_signal = fvec(win_s)
     169        input_signal[win_s//2] = 1
     170        c = f(input_signal)
     171        output_signal = f.rdo(c)
     172        assert_almost_equal(input_signal, output_signal)
     173
     174class aubio_fft_wrong_params(TestCase):
     175
    145176    def test_large_input_timegrain(self):
    146177        win_s = 1024
     
    170201        with self.assertRaises(ValueError):
    171202            f.rdo(s)
    172 
    173 class aubio_fft_wrong_params(TestCase):
    174203
    175204    def test_wrong_buf_size(self):
     
    177206        with self.assertRaises(ValueError):
    178207            fft(win_s)
    179 
    180     def test_buf_size_not_power_of_two(self):
    181         # when compiled with fftw3, aubio supports non power of two fft sizes
    182         win_s = 320
    183         try:
    184             with self.assertRaises(RuntimeError):
    185                 fft(win_s)
    186         except AssertionError:
    187             self.skipTest('creating aubio.fft with size %d did not fail' % win_s)
    188208
    189209    def test_buf_size_too_small(self):
  • src/spectral/fft.c

    r9ef3c6e r01d4d19  
    9090#define aubio_vDSP_vsadd               vDSP_vsadd
    9191#define aubio_vDSP_vsmul               vDSP_vsmul
    92 #define aubio_vDSP_create_fftsetup     vDSP_create_fftsetup
    93 #define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetup
    9492#define aubio_DSPComplex               DSPComplex
    9593#define aubio_DSPSplitComplex          DSPSplitComplex
    96 #define aubio_FFTSetup                 FFTSetup
     94#define aubio_vDSP_DFT_Setup           vDSP_DFT_Setup
     95#define aubio_vDSP_DFT_zrop_CreateSetup vDSP_DFT_zrop_CreateSetup
     96#define aubio_vDSP_DFT_Execute         vDSP_DFT_Execute
     97#define aubio_vDSP_DFT_DestroySetup    vDSP_DFT_DestroySetup
    9798#define aubio_vvsqrt                   vvsqrtf
    9899#else
     
    104105#define aubio_vDSP_vsadd               vDSP_vsaddD
    105106#define aubio_vDSP_vsmul               vDSP_vsmulD
    106 #define aubio_vDSP_create_fftsetup     vDSP_create_fftsetupD
    107 #define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetupD
    108107#define aubio_DSPComplex               DSPDoubleComplex
    109108#define aubio_DSPSplitComplex          DSPDoubleSplitComplex
    110 #define aubio_FFTSetup                 FFTSetupD
     109#define aubio_vDSP_DFT_Setup           vDSP_DFT_SetupD
     110#define aubio_vDSP_DFT_zrop_CreateSetup vDSP_DFT_zrop_CreateSetupD
     111#define aubio_vDSP_DFT_Execute         vDSP_DFT_ExecuteD
     112#define aubio_vDSP_DFT_DestroySetup    vDSP_DFT_DestroySetupD
    111113#define aubio_vvsqrt                   vvsqrt
    112114#endif /* HAVE_AUBIO_DOUBLE */
     
    153155
    154156#elif defined HAVE_ACCELERATE  // using ACCELERATE
    155   int log2fftsize;
    156   aubio_FFTSetup fftSetup;
     157  aubio_vDSP_DFT_Setup fftSetupFwd;
     158  aubio_vDSP_DFT_Setup fftSetupBwd;
    157159  aubio_DSPSplitComplex spec;
    158160  smpl_t *in, *out;
     
    211213
    212214#elif defined HAVE_ACCELERATE  // using ACCELERATE
     215  {
     216    uint_t radix = winsize;
     217    uint_t order = 0;
     218    while ((radix / 2) * 2 == radix) {
     219      radix /= 2;
     220      order++;
     221    }
     222    if (order < 4 || (radix != 1 && radix != 3 && radix != 5 && radix != 15)) {
     223      AUBIO_ERR("fft: vDSP/Accelerate supports FFT with sizes = "
     224          "f * 2 ** n, where n > 4 and f in [1, 3, 5, 15], but requested %d. "
     225          "Use the closest power of two, or try recompiling aubio with "
     226          "--enable-fftw3.\n", winsize);
     227      goto beach;
     228    }
     229  }
    213230  s->winsize = winsize;
    214231  s->fft_size = winsize;
    215232  s->compspec = new_fvec(winsize);
    216   s->log2fftsize = aubio_power_of_two_order(s->fft_size);
    217233  s->in = AUBIO_ARRAY(smpl_t, s->fft_size);
    218234  s->out = AUBIO_ARRAY(smpl_t, s->fft_size);
    219235  s->spec.realp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
    220236  s->spec.imagp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
    221   s->fftSetup = aubio_vDSP_create_fftsetup(s->log2fftsize, FFT_RADIX2);
     237  s->fftSetupFwd = aubio_vDSP_DFT_zrop_CreateSetup(NULL,
     238      s->fft_size, vDSP_DFT_FORWARD);
     239  s->fftSetupBwd = aubio_vDSP_DFT_zrop_CreateSetup(s->fftSetupFwd,
     240      s->fft_size, vDSP_DFT_INVERSE);
    222241
    223242#elif defined HAVE_INTEL_IPP  // using Intel IPP
     
    293312  AUBIO_FREE(s->spec.realp);
    294313  AUBIO_FREE(s->spec.imagp);
    295   aubio_vDSP_destroy_fftsetup(s->fftSetup);
     314  aubio_vDSP_DFT_DestroySetup(s->fftSetupBwd);
     315  aubio_vDSP_DFT_DestroySetup(s->fftSetupFwd);
    296316
    297317#elif defined HAVE_INTEL_IPP  // using Intel IPP
     
    351371  aubio_vDSP_ctoz((aubio_DSPComplex*)s->in, 2, &s->spec, 1, s->fft_size/2);
    352372  // compute the FFT
    353   aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
     373  aubio_vDSP_DFT_Execute(s->fftSetupFwd, s->spec.realp, s->spec.imagp,
     374      s->spec.realp, s->spec.imagp);
    354375  // convert from vDSP complex split to [ r0, r1, ..., rN, iN-1, .., i2, i1]
    355376  compspec->data[0] = s->spec.realp[0];
     
    419440  aubio_vDSP_ctoz((aubio_DSPComplex*)s->out, 2, &s->spec, 1, s->fft_size/2);
    420441  // compute the FFT
    421   aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
     442  aubio_vDSP_DFT_Execute(s->fftSetupBwd, s->spec.realp, s->spec.imagp,
     443      s->spec.realp, s->spec.imagp);
    422444  // convert result to real output
    423445  aubio_vDSP_ztoc(&s->spec, 1, (aubio_DSPComplex*)output->data, 2, s->fft_size/2);
     
    494516  }
    495517#endif
    496   if (compspec->data[compspec->length/2] < 0) {
    497     spectrum->phas[spectrum->length - 1] = PI;
     518#ifdef HAVE_FFTW3
     519  // for even length only, make sure last element is 0 or PI
     520  if (2 * (compspec->length / 2) == compspec->length) {
     521#endif
     522    if (compspec->data[compspec->length/2] < 0) {
     523      spectrum->phas[spectrum->length - 1] = PI;
     524    } else {
     525      spectrum->phas[spectrum->length - 1] = 0.;
     526    }
     527#ifdef HAVE_FFTW3
    498528  } else {
    499     spectrum->phas[spectrum->length - 1] = 0.;
    500   }
     529    i = spectrum->length - 1;
     530    spectrum->phas[i] = ATAN2(compspec->data[compspec->length-i],
     531        compspec->data[i]);
     532  }
     533#endif
    501534}
    502535
     
    508541        + SQR(compspec->data[compspec->length - i]) );
    509542  }
    510   spectrum->norm[spectrum->length-1] =
    511     ABS(compspec->data[compspec->length/2]);
     543#ifdef HAVE_FFTW3
     544  // for even length, make sure last element is > 0
     545  if (2 * (compspec->length / 2) == compspec->length) {
     546#endif
     547    spectrum->norm[spectrum->length-1] =
     548      ABS(compspec->data[compspec->length/2]);
     549#ifdef HAVE_FFTW3
     550  } else {
     551    i = spectrum->length - 1;
     552    spectrum->norm[i] = SQRT(SQR(compspec->data[i])
     553        + SQR(compspec->data[compspec->length - i]) );
     554  }
     555#endif
    512556}
    513557
Note: See TracChangeset for help on using the changeset viewer.