1 /* Copyright (C) 2005-2006 Jean-Marc Valin 2 File: fftwrap.c 3 4 Wrapper for various FFTs 5 6 Redistribution and use in source and binary forms, with or without 7 modification, are permitted provided that the following conditions 8 are met: 9 10 - Redistributions of source code must retain the above copyright 11 notice, this list of conditions and the following disclaimer. 12 13 - Redistributions in binary form must reproduce the above copyright 14 notice, this list of conditions and the following disclaimer in the 15 documentation and/or other materials provided with the distribution. 16 17 - Neither the name of the Xiph.org Foundation nor the names of its 18 contributors may be used to endorse or promote products derived from 19 this software without specific prior written permission. 20 21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 24 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR 25 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 26 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 27 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 28 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 29 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 30 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 32 33 */ 34 35 #ifdef HAVE_CONFIG_H 36 #include "config.h" 37 #endif 38 39 #include "arch.h" 40 #include "os_support.h" 41 42 #define MAX_FFT_SIZE 2048 43 44 #ifdef FIXED_POINT 45 static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len) 46 { 47 int i, shift; 48 spx_word16_t max_val = 0; 49 for (i=0;i<len;i++) 50 { 51 if (in[i]>max_val) 52 max_val = in[i]; 53 if (-in[i]>max_val) 54 max_val = -in[i]; 55 } 56 shift=0; 57 while (max_val <= (bound>>1) && max_val != 0) 58 { 59 max_val <<= 1; 60 shift++; 61 } 62 for (i=0;i<len;i++) 63 { 64 out[i] = SHL16(in[i], shift); 65 } 66 return shift; 67 } 68 69 static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len) 70 { 71 int i; 72 for (i=0;i<len;i++) 73 { 74 out[i] = PSHR16(in[i], shift); 75 } 76 } 77 #endif 78 79 #ifdef USE_SMALLFT 80 81 #include "smallft.h" 82 #include <math.h> 83 84 void *spx_fft_init(int size) 85 { 86 struct drft_lookup *table; 87 table = speex_alloc(sizeof(struct drft_lookup)); 88 spx_drft_init((struct drft_lookup *)table, size); 89 return (void*)table; 90 } 91 92 void spx_fft_destroy(void *table) 93 { 94 spx_drft_clear(table); 95 speex_free(table); 96 } 97 98 void spx_fft(void *table, float *in, float *out) 99 { 100 if (in==out) 101 { 102 int i; 103 float scale = 1./((struct drft_lookup *)table)->n; 104 speex_warning("FFT should not be done in-place"); 105 for (i=0;i<((struct drft_lookup *)table)->n;i++) 106 out[i] = scale*in[i]; 107 } else { 108 int i; 109 float scale = 1./((struct drft_lookup *)table)->n; 110 for (i=0;i<((struct drft_lookup *)table)->n;i++) 111 out[i] = scale*in[i]; 112 } 113 spx_drft_forward((struct drft_lookup *)table, out); 114 } 115 116 void spx_ifft(void *table, float *in, float *out) 117 { 118 if (in==out) 119 { 120 speex_warning("FFT should not be done in-place"); 121 } else { 122 int i; 123 for (i=0;i<((struct drft_lookup *)table)->n;i++) 124 out[i] = in[i]; 125 } 126 spx_drft_backward((struct drft_lookup *)table, out); 127 } 128 129 #elif defined(USE_INTEL_MKL) 130 #include <mkl.h> 131 132 struct mkl_config { 133 DFTI_DESCRIPTOR_HANDLE desc; 134 int N; 135 }; 136 137 void *spx_fft_init(int size) 138 { 139 struct mkl_config *table = (struct mkl_config *) speex_alloc(sizeof(struct mkl_config)); 140 table->N = size; 141 DftiCreateDescriptor(&table->desc, DFTI_SINGLE, DFTI_REAL, 1, size); 142 DftiSetValue(table->desc, DFTI_PACKED_FORMAT, DFTI_PACK_FORMAT); 143 DftiSetValue(table->desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE); 144 DftiSetValue(table->desc, DFTI_FORWARD_SCALE, 1.0f / size); 145 DftiCommitDescriptor(table->desc); 146 return table; 147 } 148 149 void spx_fft_destroy(void *table) 150 { 151 struct mkl_config *t = (struct mkl_config *) table; 152 DftiFreeDescriptor(t->desc); 153 speex_free(table); 154 } 155 156 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 157 { 158 struct mkl_config *t = (struct mkl_config *) table; 159 DftiComputeForward(t->desc, in, out); 160 } 161 162 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 163 { 164 struct mkl_config *t = (struct mkl_config *) table; 165 DftiComputeBackward(t->desc, in, out); 166 } 167 168 #elif defined(USE_INTEL_IPP) 169 170 #include <ipps.h> 171 172 struct ipp_fft_config 173 { 174 IppsDFTSpec_R_32f *dftSpec; 175 Ipp8u *buffer; 176 }; 177 178 void *spx_fft_init(int size) 179 { 180 int bufferSize = 0; 181 int hint; 182 struct ipp_fft_config *table; 183 184 table = (struct ipp_fft_config *)speex_alloc(sizeof(struct ipp_fft_config)); 185 186 /* there appears to be no performance difference between ippAlgHintFast and 187 ippAlgHintAccurate when using the with the floating point version 188 of the fft. */ 189 hint = ippAlgHintAccurate; 190 191 ippsDFTInitAlloc_R_32f(&table->dftSpec, size, IPP_FFT_DIV_FWD_BY_N, hint); 192 193 ippsDFTGetBufSize_R_32f(table->dftSpec, &bufferSize); 194 table->buffer = ippsMalloc_8u(bufferSize); 195 196 return table; 197 } 198 199 void spx_fft_destroy(void *table) 200 { 201 struct ipp_fft_config *t = (struct ipp_fft_config *)table; 202 ippsFree(t->buffer); 203 ippsDFTFree_R_32f(t->dftSpec); 204 speex_free(t); 205 } 206 207 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 208 { 209 struct ipp_fft_config *t = (struct ipp_fft_config *)table; 210 ippsDFTFwd_RToPack_32f(in, out, t->dftSpec, t->buffer); 211 } 212 213 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 214 { 215 struct ipp_fft_config *t = (struct ipp_fft_config *)table; 216 ippsDFTInv_PackToR_32f(in, out, t->dftSpec, t->buffer); 217 } 218 219 #elif defined(USE_GPL_FFTW3) 220 221 #include <fftw3.h> 222 223 struct fftw_config { 224 float *in; 225 float *out; 226 fftwf_plan fft; 227 fftwf_plan ifft; 228 int N; 229 }; 230 231 void *spx_fft_init(int size) 232 { 233 struct fftw_config *table = (struct fftw_config *) speex_alloc(sizeof(struct fftw_config)); 234 table->in = fftwf_malloc(sizeof(float) * (size+2)); 235 table->out = fftwf_malloc(sizeof(float) * (size+2)); 236 237 table->fft = fftwf_plan_dft_r2c_1d(size, table->in, (fftwf_complex *) table->out, FFTW_PATIENT); 238 table->ifft = fftwf_plan_dft_c2r_1d(size, (fftwf_complex *) table->in, table->out, FFTW_PATIENT); 239 240 table->N = size; 241 return table; 242 } 243 244 void spx_fft_destroy(void *table) 245 { 246 struct fftw_config *t = (struct fftw_config *) table; 247 fftwf_destroy_plan(t->fft); 248 fftwf_destroy_plan(t->ifft); 249 fftwf_free(t->in); 250 fftwf_free(t->out); 251 speex_free(table); 252 } 253 254 255 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 256 { 257 int i; 258 struct fftw_config *t = (struct fftw_config *) table; 259 const int N = t->N; 260 float *iptr = t->in; 261 float *optr = t->out; 262 const float m = 1.0 / N; 263 for(i=0;i<N;++i) 264 iptr[i]=in[i] * m; 265 266 fftwf_execute(t->fft); 267 268 out[0] = optr[0]; 269 for(i=1;i<N;++i) 270 out[i] = optr[i+1]; 271 } 272 273 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 274 { 275 int i; 276 struct fftw_config *t = (struct fftw_config *) table; 277 const int N = t->N; 278 float *iptr = t->in; 279 float *optr = t->out; 280 281 iptr[0] = in[0]; 282 iptr[1] = 0.0f; 283 for(i=1;i<N;++i) 284 iptr[i+1] = in[i]; 285 iptr[N+1] = 0.0f; 286 287 fftwf_execute(t->ifft); 288 289 for(i=0;i<N;++i) 290 out[i] = optr[i]; 291 } 292 293 #elif defined(USE_KISS_FFT) 294 295 #include "kiss_fftr.h" 296 #include "kiss_fft.h" 297 298 struct kiss_config { 299 kiss_fftr_cfg forward; 300 kiss_fftr_cfg backward; 301 int N; 302 }; 303 304 void *spx_fft_init(int size) 305 { 306 struct kiss_config *table; 307 table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config)); 308 table->forward = kiss_fftr_alloc(size,0,NULL,NULL); 309 table->backward = kiss_fftr_alloc(size,1,NULL,NULL); 310 table->N = size; 311 return table; 312 } 313 314 void spx_fft_destroy(void *table) 315 { 316 struct kiss_config *t = (struct kiss_config *)table; 317 kiss_fftr_free(t->forward); 318 kiss_fftr_free(t->backward); 319 speex_free(table); 320 } 321 322 #ifdef FIXED_POINT 323 324 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 325 { 326 int shift; 327 struct kiss_config *t = (struct kiss_config *)table; 328 shift = maximize_range(in, in, 32000, t->N); 329 kiss_fftr2(t->forward, in, out); 330 renorm_range(in, in, shift, t->N); 331 renorm_range(out, out, shift, t->N); 332 } 333 334 #else 335 336 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 337 { 338 int i; 339 float scale; 340 struct kiss_config *t = (struct kiss_config *)table; 341 scale = 1./t->N; 342 kiss_fftr2(t->forward, in, out); 343 for (i=0;i<t->N;i++) 344 out[i] *= scale; 345 } 346 #endif 347 348 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 349 { 350 struct kiss_config *t = (struct kiss_config *)table; 351 kiss_fftri2(t->backward, in, out); 352 } 353 354 355 #else 356 357 #error No other FFT implemented 358 359 #endif 360 361 362 #ifdef FIXED_POINT 363 /*#include "smallft.h"*/ 364 365 366 void spx_fft_float(void *table, float *in, float *out) 367 { 368 int i; 369 #ifdef USE_SMALLFT 370 int N = ((struct drft_lookup *)table)->n; 371 #elif defined(USE_KISS_FFT) 372 int N = ((struct kiss_config *)table)->N; 373 #else 374 #endif 375 #ifdef VAR_ARRAYS 376 spx_word16_t _in[N]; 377 spx_word16_t _out[N]; 378 #else 379 spx_word16_t _in[MAX_FFT_SIZE]; 380 spx_word16_t _out[MAX_FFT_SIZE]; 381 #endif 382 for (i=0;i<N;i++) 383 _in[i] = (int)floor(.5+in[i]); 384 spx_fft(table, _in, _out); 385 for (i=0;i<N;i++) 386 out[i] = _out[i]; 387 #if 0 388 if (!fixed_point) 389 { 390 float scale; 391 struct drft_lookup t; 392 spx_drft_init(&t, ((struct kiss_config *)table)->N); 393 scale = 1./((struct kiss_config *)table)->N; 394 for (i=0;i<((struct kiss_config *)table)->N;i++) 395 out[i] = scale*in[i]; 396 spx_drft_forward(&t, out); 397 spx_drft_clear(&t); 398 } 399 #endif 400 } 401 402 void spx_ifft_float(void *table, float *in, float *out) 403 { 404 int i; 405 #ifdef USE_SMALLFT 406 int N = ((struct drft_lookup *)table)->n; 407 #elif defined(USE_KISS_FFT) 408 int N = ((struct kiss_config *)table)->N; 409 #else 410 #endif 411 #ifdef VAR_ARRAYS 412 spx_word16_t _in[N]; 413 spx_word16_t _out[N]; 414 #else 415 spx_word16_t _in[MAX_FFT_SIZE]; 416 spx_word16_t _out[MAX_FFT_SIZE]; 417 #endif 418 for (i=0;i<N;i++) 419 _in[i] = (int)floor(.5+in[i]); 420 spx_ifft(table, _in, _out); 421 for (i=0;i<N;i++) 422 out[i] = _out[i]; 423 #if 0 424 if (!fixed_point) 425 { 426 int i; 427 struct drft_lookup t; 428 spx_drft_init(&t, ((struct kiss_config *)table)->N); 429 for (i=0;i<((struct kiss_config *)table)->N;i++) 430 out[i] = in[i]; 431 spx_drft_backward(&t, out); 432 spx_drft_clear(&t); 433 } 434 #endif 435 } 436 437 #else 438 439 void spx_fft_float(void *table, float *in, float *out) 440 { 441 spx_fft(table, in, out); 442 } 443 void spx_ifft_float(void *table, float *in, float *out) 444 { 445 spx_ifft(table, in, out); 446 } 447 448 #endif 449