Changeset 633400d for src/spectral/fft.c
- Timestamp:
- Dec 5, 2018, 10:34:39 PM (5 years ago)
- Branches:
- feature/cnn, feature/crepe, feature/pitchshift, feature/timestretch, fix/ffmpeg5, master
- Children:
- 283a619a
- Parents:
- 5b46bc3 (diff), f19db54 (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. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/spectral/fft.c
r5b46bc3 r633400d 78 78 pthread_mutex_t aubio_fftw_mutex = PTHREAD_MUTEX_INITIALIZER; 79 79 80 #else 81 #ifdef HAVE_ACCELERATE // using ACCELERATE 80 #elif defined HAVE_ACCELERATE // using ACCELERATE 82 81 // https://developer.apple.com/library/mac/#documentation/Accelerate/Reference/vDSPRef/Reference/reference.html 83 82 #include <Accelerate/Accelerate.h> … … 91 90 #define aubio_vDSP_vsadd vDSP_vsadd 92 91 #define aubio_vDSP_vsmul vDSP_vsmul 93 #define aubio_vDSP_create_fftsetup vDSP_create_fftsetup94 #define aubio_vDSP_destroy_fftsetup vDSP_destroy_fftsetup95 92 #define aubio_DSPComplex DSPComplex 96 93 #define aubio_DSPSplitComplex DSPSplitComplex 97 #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 98 98 #define aubio_vvsqrt vvsqrtf 99 99 #else … … 105 105 #define aubio_vDSP_vsadd vDSP_vsaddD 106 106 #define aubio_vDSP_vsmul vDSP_vsmulD 107 #define aubio_vDSP_create_fftsetup vDSP_create_fftsetupD108 #define aubio_vDSP_destroy_fftsetup vDSP_destroy_fftsetupD109 107 #define aubio_DSPComplex DSPDoubleComplex 110 108 #define aubio_DSPSplitComplex DSPDoubleSplitComplex 111 #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 112 113 #define aubio_vvsqrt vvsqrt 113 114 #endif /* HAVE_AUBIO_DOUBLE */ 114 115 115 #else // using OOURA 116 #elif defined HAVE_INTEL_IPP // using INTEL IPP 117 118 #if !HAVE_AUBIO_DOUBLE 119 #define aubio_IppFloat Ipp32f 120 #define aubio_IppComplex Ipp32fc 121 #define aubio_FFTSpec FFTSpec_R_32f 122 #define aubio_ippsMalloc_complex ippsMalloc_32fc 123 #define aubio_ippsFFTInit_R ippsFFTInit_R_32f 124 #define aubio_ippsFFTGetSize_R ippsFFTGetSize_R_32f 125 #define aubio_ippsFFTInv_CCSToR ippsFFTInv_CCSToR_32f 126 #define aubio_ippsFFTFwd_RToCCS ippsFFTFwd_RToCCS_32f 127 #define aubio_ippsAtan2 ippsAtan2_32f_A21 128 #else /* HAVE_AUBIO_DOUBLE */ 129 #define aubio_IppFloat Ipp64f 130 #define aubio_IppComplex Ipp64fc 131 #define aubio_FFTSpec FFTSpec_R_64f 132 #define aubio_ippsMalloc_complex ippsMalloc_64fc 133 #define aubio_ippsFFTInit_R ippsFFTInit_R_64f 134 #define aubio_ippsFFTGetSize_R ippsFFTGetSize_R_64f 135 #define aubio_ippsFFTInv_CCSToR ippsFFTInv_CCSToR_64f 136 #define aubio_ippsFFTFwd_RToCCS ippsFFTFwd_RToCCS_64f 137 #define aubio_ippsAtan2 ippsAtan2_64f_A50 138 #endif 139 140 141 #else // using OOURA 116 142 // let's use ooura instead 117 143 extern void aubio_ooura_rdft(int, int, smpl_t *, int *, smpl_t *); 118 144 119 #endif /* HAVE_ACCELERATE */ 120 #endif /* HAVE_FFTW3 */ 145 #endif 121 146 122 147 struct _aubio_fft_t { 123 148 uint_t winsize; 124 149 uint_t fft_size; 150 125 151 #ifdef HAVE_FFTW3 // using FFTW3 126 152 real_t *in, *out; 127 153 fftw_plan pfw, pbw; 128 fft_data_t * specdata; 129 #else 130 # ifdef HAVE_ACCELERATE// using ACCELERATE131 int log2fftsize;132 aubio_ FFTSetup fftSetup;154 fft_data_t * specdata; /* complex spectral data */ 155 156 #elif defined HAVE_ACCELERATE // using ACCELERATE 157 aubio_vDSP_DFT_Setup fftSetupFwd; 158 aubio_vDSP_DFT_Setup fftSetupBwd; 133 159 aubio_DSPSplitComplex spec; 134 160 smpl_t *in, *out; 161 162 #elif defined HAVE_INTEL_IPP // using Intel IPP 163 smpl_t *in, *out; 164 Ipp8u* memSpec; 165 Ipp8u* memInit; 166 Ipp8u* memBuffer; 167 struct aubio_FFTSpec* fftSpec; 168 aubio_IppComplex* complexOut; 135 169 #else // using OOURA 136 170 smpl_t *in, *out; 137 171 smpl_t *w; 138 172 int *ip; 139 #endif /* HAVE_ACCELERATE*/140 #endif /* HAVE_FFTW3 */ 173 #endif /* using OOURA */ 174 141 175 fvec_t * compspec; 142 176 }; … … 148 182 goto beach; 149 183 } 184 150 185 #ifdef HAVE_FFTW3 151 186 uint_t i; … … 176 211 s->specdata[i] = 0.; 177 212 } 178 #else 179 #ifdef HAVE_ACCELERATE // using ACCELERATE 213 214 #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 } 180 230 s->winsize = winsize; 181 231 s->fft_size = winsize; 182 232 s->compspec = new_fvec(winsize); 183 s->log2fftsize = (uint_t)log2f(s->fft_size);184 233 s->in = AUBIO_ARRAY(smpl_t, s->fft_size); 185 234 s->out = AUBIO_ARRAY(smpl_t, s->fft_size); 186 235 s->spec.realp = AUBIO_ARRAY(smpl_t, s->fft_size/2); 187 236 s->spec.imagp = AUBIO_ARRAY(smpl_t, s->fft_size/2); 188 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); 241 242 #elif defined HAVE_INTEL_IPP // using Intel IPP 243 const IppHintAlgorithm qualityHint = ippAlgHintAccurate; // OR ippAlgHintFast; 244 const int flags = IPP_FFT_NODIV_BY_ANY; // we're scaling manually afterwards 245 int order = aubio_power_of_two_order(winsize); 246 int sizeSpec, sizeInit, sizeBuffer; 247 IppStatus status; 248 249 if (winsize <= 4 || aubio_is_power_of_two(winsize) != 1) 250 { 251 AUBIO_ERR("intel IPP fft: can only create with sizes > 4 and power of two, requested %d," 252 " try recompiling aubio with --enable-fftw3\n", winsize); 253 goto beach; 254 } 255 256 status = aubio_ippsFFTGetSize_R(order, flags, qualityHint, 257 &sizeSpec, &sizeInit, &sizeBuffer); 258 if (status != ippStsNoErr) { 259 AUBIO_ERR("fft: failed to initialize fft. IPP error: %d\n", status); 260 goto beach; 261 } 262 s->fft_size = s->winsize = winsize; 263 s->compspec = new_fvec(winsize); 264 s->in = AUBIO_ARRAY(smpl_t, s->winsize); 265 s->out = AUBIO_ARRAY(smpl_t, s->winsize); 266 s->memSpec = ippsMalloc_8u(sizeSpec); 267 s->memBuffer = ippsMalloc_8u(sizeBuffer); 268 if (sizeInit > 0 ) { 269 s->memInit = ippsMalloc_8u(sizeInit); 270 } 271 s->complexOut = aubio_ippsMalloc_complex(s->fft_size / 2 + 1); 272 status = aubio_ippsFFTInit_R( 273 &s->fftSpec, order, flags, qualityHint, s->memSpec, s->memInit); 274 if (status != ippStsNoErr) { 275 AUBIO_ERR("fft: failed to initialize. IPP error: %d\n", status); 276 goto beach; 277 } 278 189 279 #else // using OOURA 190 280 if (aubio_is_power_of_two(winsize) != 1) { … … 201 291 s->w = AUBIO_ARRAY(smpl_t, s->fft_size); 202 292 s->ip[0] = 0; 203 #endif /* HAVE_ACCELERATE*/204 #endif /* HAVE_FFTW3 */ 293 #endif /* using OOURA */ 294 205 295 return s; 296 206 297 beach: 207 298 AUBIO_FREE(s); … … 211 302 void del_aubio_fft(aubio_fft_t * s) { 212 303 /* destroy data */ 213 del_fvec(s->compspec);214 304 #ifdef HAVE_FFTW3 // using FFTW3 215 305 pthread_mutex_lock(&aubio_fftw_mutex); … … 218 308 fftw_free(s->specdata); 219 309 pthread_mutex_unlock(&aubio_fftw_mutex); 220 #else /* HAVE_FFTW3 */ 221 # ifdef HAVE_ACCELERATE// using ACCELERATE310 311 #elif defined HAVE_ACCELERATE // using ACCELERATE 222 312 AUBIO_FREE(s->spec.realp); 223 313 AUBIO_FREE(s->spec.imagp); 224 aubio_vDSP_destroy_fftsetup(s->fftSetup); 314 aubio_vDSP_DFT_DestroySetup(s->fftSetupBwd); 315 aubio_vDSP_DFT_DestroySetup(s->fftSetupFwd); 316 317 #elif defined HAVE_INTEL_IPP // using Intel IPP 318 ippFree(s->memSpec); 319 ippFree(s->memInit); 320 ippFree(s->memBuffer); 321 ippFree(s->complexOut); 322 225 323 #else // using OOURA 226 324 AUBIO_FREE(s->w); 227 325 AUBIO_FREE(s->ip); 228 #endif /* HAVE_ACCELERATE */ 229 #endif /* HAVE_FFTW3 */ 326 #endif 327 328 del_fvec(s->compspec); 329 AUBIO_FREE(s->in); 230 330 AUBIO_FREE(s->out); 231 AUBIO_FREE(s->in);232 331 AUBIO_FREE(s); 233 332 } … … 252 351 memcpy(s->in, input->data, s->winsize * sizeof(smpl_t)); 253 352 #endif /* HAVE_MEMCPY_HACKS */ 353 254 354 #ifdef HAVE_FFTW3 // using FFTW3 255 355 fftw_execute(s->pfw); … … 266 366 } 267 367 #endif /* HAVE_COMPLEX_H */ 268 #else /* HAVE_FFTW3 */ 269 # ifdef HAVE_ACCELERATE// using ACCELERATE368 369 #elif defined HAVE_ACCELERATE // using ACCELERATE 270 370 // convert real data to even/odd format used in vDSP 271 371 aubio_vDSP_ctoz((aubio_DSPComplex*)s->in, 2, &s->spec, 1, s->fft_size/2); 272 372 // compute the FFT 273 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); 274 375 // convert from vDSP complex split to [ r0, r1, ..., rN, iN-1, .., i2, i1] 275 376 compspec->data[0] = s->spec.realp[0]; … … 282 383 smpl_t scale = 1./2.; 283 384 aubio_vDSP_vsmul(compspec->data, 1, &scale, compspec->data, 1, s->fft_size); 385 386 #elif defined HAVE_INTEL_IPP // using Intel IPP 387 388 // apply fft 389 aubio_ippsFFTFwd_RToCCS(s->in, (aubio_IppFloat*)s->complexOut, s->fftSpec, s->memBuffer); 390 // convert complex buffer to [ r0, r1, ..., rN, iN-1, .., i2, i1] 391 compspec->data[0] = s->complexOut[0].re; 392 compspec->data[s->fft_size / 2] = s->complexOut[s->fft_size / 2].re; 393 for (i = 1; i < s->fft_size / 2; i++) { 394 compspec->data[i] = s->complexOut[i].re; 395 compspec->data[s->fft_size - i] = s->complexOut[i].im; 396 } 397 284 398 #else // using OOURA 285 399 aubio_ooura_rdft(s->winsize, 1, s->in, s->ip, s->w); … … 290 404 compspec->data[s->winsize - i] = - s->in[2 * i + 1]; 291 405 } 292 #endif /* HAVE_ACCELERATE */ 293 #endif /* HAVE_FFTW3 */ 406 #endif /* using OOURA */ 294 407 } 295 408 … … 314 427 output->data[i] = s->out[i]*renorm; 315 428 } 316 #else /* HAVE_FFTW3 */ 317 # ifdef HAVE_ACCELERATE// using ACCELERATE429 430 #elif defined HAVE_ACCELERATE // using ACCELERATE 318 431 // convert from real imag [ r0, r1, ..., rN, iN-1, .., i2, i1] 319 432 // to vDSP packed format [ r0, rN, r1, i1, ..., rN-1, iN-1 ] … … 327 440 aubio_vDSP_ctoz((aubio_DSPComplex*)s->out, 2, &s->spec, 1, s->fft_size/2); 328 441 // compute the FFT 329 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); 330 444 // convert result to real output 331 445 aubio_vDSP_ztoc(&s->spec, 1, (aubio_DSPComplex*)output->data, 2, s->fft_size/2); … … 333 447 smpl_t scale = 1.0 / s->winsize; 334 448 aubio_vDSP_vsmul(output->data, 1, &scale, output->data, 1, s->fft_size); 449 450 #elif defined HAVE_INTEL_IPP // using Intel IPP 451 452 // convert from real imag [ r0, 0, ..., rN, iN-1, .., i2, i1, iN-1] to complex format 453 s->complexOut[0].re = compspec->data[0]; 454 s->complexOut[0].im = 0; 455 s->complexOut[s->fft_size / 2].re = compspec->data[s->fft_size / 2]; 456 s->complexOut[s->fft_size / 2].im = 0.0; 457 for (i = 1; i < s->fft_size / 2; i++) { 458 s->complexOut[i].re = compspec->data[i]; 459 s->complexOut[i].im = compspec->data[s->fft_size - i]; 460 } 461 // apply fft 462 aubio_ippsFFTInv_CCSToR((const aubio_IppFloat *)s->complexOut, output->data, s->fftSpec, s->memBuffer); 463 // apply scaling 464 aubio_ippsMulC(output->data, 1.0 / s->winsize, output->data, s->fft_size); 465 335 466 #else // using OOURA 336 467 smpl_t scale = 2.0 / s->winsize; … … 345 476 output->data[i] = s->out[i] * scale; 346 477 } 347 #endif /* HAVE_ACCELERATE */ 348 #endif /* HAVE_FFTW3 */ 478 #endif 349 479 } 350 480 … … 366 496 spectrum->phas[0] = 0.; 367 497 } 498 #if defined(HAVE_INTEL_IPP) 499 // convert from real imag [ r0, r1, ..., rN, iN-1, ..., i2, i1, i0] 500 // to [ r0, r1, ..., rN, i0, i1, i2, ..., iN-1] 501 for (i = 1; i < spectrum->length / 2; i++) { 502 ELEM_SWAP(compspec->data[compspec->length - i], 503 compspec->data[spectrum->length + i - 1]); 504 } 505 aubio_ippsAtan2(compspec->data + spectrum->length, 506 compspec->data + 1, spectrum->phas + 1, spectrum->length - 1); 507 // revert the imaginary part back again 508 for (i = 1; i < spectrum->length / 2; i++) { 509 ELEM_SWAP(compspec->data[spectrum->length + i - 1], 510 compspec->data[compspec->length - i]); 511 } 512 #else 368 513 for (i=1; i < spectrum->length - 1; i++) { 369 514 spectrum->phas[i] = ATAN2(compspec->data[compspec->length-i], 370 515 compspec->data[i]); 371 516 } 372 if (compspec->data[compspec->length/2] < 0) { 373 spectrum->phas[spectrum->length - 1] = PI; 517 #endif 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 374 528 } else { 375 spectrum->phas[spectrum->length - 1] = 0.; 376 } 529 i = spectrum->length - 1; 530 spectrum->phas[i] = ATAN2(compspec->data[compspec->length-i], 531 compspec->data[i]); 532 } 533 #endif 377 534 } 378 535 … … 384 541 + SQR(compspec->data[compspec->length - i]) ); 385 542 } 386 spectrum->norm[spectrum->length-1] = 387 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 388 556 } 389 557
Note: See TracChangeset
for help on using the changeset viewer.