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_GPL_FFTW3) 169 170 #include <fftw3.h> 171 172 struct fftw_config { 173 float *in; 174 float *out; 175 fftwf_plan fft; 176 fftwf_plan ifft; 177 int N; 178 }; 179 180 void *spx_fft_init(int size) 181 { 182 struct fftw_config *table = (struct fftw_config *) speex_alloc(sizeof(struct fftw_config)); 183 table->in = fftwf_malloc(sizeof(float) * (size+2)); 184 table->out = fftwf_malloc(sizeof(float) * (size+2)); 185 186 table->fft = fftwf_plan_dft_r2c_1d(size, table->in, (fftwf_complex *) table->out, FFTW_PATIENT); 187 table->ifft = fftwf_plan_dft_c2r_1d(size, (fftwf_complex *) table->in, table->out, FFTW_PATIENT); 188 189 table->N = size; 190 return table; 191 } 192 193 void spx_fft_destroy(void *table) 194 { 195 struct fftw_config *t = (struct fftw_config *) table; 196 fftwf_destroy_plan(t->fft); 197 fftwf_destroy_plan(t->ifft); 198 fftwf_free(t->in); 199 fftwf_free(t->out); 200 speex_free(table); 201 } 202 203 204 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 205 { 206 int i; 207 struct fftw_config *t = (struct fftw_config *) table; 208 const int N = t->N; 209 float *iptr = t->in; 210 float *optr = t->out; 211 const float m = 1.0 / N; 212 for(i=0;i<N;++i) 213 iptr[i]=in[i] * m; 214 215 fftwf_execute(t->fft); 216 217 out[0] = optr[0]; 218 for(i=1;i<N;++i) 219 out[i] = optr[i+1]; 220 } 221 222 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 223 { 224 int i; 225 struct fftw_config *t = (struct fftw_config *) table; 226 const int N = t->N; 227 float *iptr = t->in; 228 float *optr = t->out; 229 230 iptr[0] = in[0]; 231 iptr[1] = 0.0f; 232 for(i=1;i<N;++i) 233 iptr[i+1] = in[i]; 234 iptr[N+1] = 0.0f; 235 236 fftwf_execute(t->ifft); 237 238 for(i=0;i<N;++i) 239 out[i] = optr[i]; 240 } 241 242 #elif defined(USE_KISS_FFT) 243 244 #include "kiss_fftr.h" 245 #include "kiss_fft.h" 246 247 struct kiss_config { 248 kiss_fftr_cfg forward; 249 kiss_fftr_cfg backward; 250 int N; 251 }; 252 253 void *spx_fft_init(int size) 254 { 255 struct kiss_config *table; 256 table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config)); 257 table->forward = kiss_fftr_alloc(size,0,NULL,NULL); 258 table->backward = kiss_fftr_alloc(size,1,NULL,NULL); 259 table->N = size; 260 return table; 261 } 262 263 void spx_fft_destroy(void *table) 264 { 265 struct kiss_config *t = (struct kiss_config *)table; 266 kiss_fftr_free(t->forward); 267 kiss_fftr_free(t->backward); 268 speex_free(table); 269 } 270 271 #ifdef FIXED_POINT 272 273 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 274 { 275 int shift; 276 struct kiss_config *t = (struct kiss_config *)table; 277 shift = maximize_range(in, in, 32000, t->N); 278 kiss_fftr2(t->forward, in, out); 279 renorm_range(in, in, shift, t->N); 280 renorm_range(out, out, shift, t->N); 281 } 282 283 #else 284 285 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) 286 { 287 int i; 288 float scale; 289 struct kiss_config *t = (struct kiss_config *)table; 290 scale = 1./t->N; 291 kiss_fftr2(t->forward, in, out); 292 for (i=0;i<t->N;i++) 293 out[i] *= scale; 294 } 295 #endif 296 297 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) 298 { 299 struct kiss_config *t = (struct kiss_config *)table; 300 kiss_fftri2(t->backward, in, out); 301 } 302 303 304 #else 305 306 #error No other FFT implemented 307 308 #endif 309 310 311 #ifdef FIXED_POINT 312 /*#include "smallft.h"*/ 313 314 315 void spx_fft_float(void *table, float *in, float *out) 316 { 317 int i; 318 #ifdef USE_SMALLFT 319 int N = ((struct drft_lookup *)table)->n; 320 #elif defined(USE_KISS_FFT) 321 int N = ((struct kiss_config *)table)->N; 322 #else 323 #endif 324 #ifdef VAR_ARRAYS 325 spx_word16_t _in[N]; 326 spx_word16_t _out[N]; 327 #else 328 spx_word16_t _in[MAX_FFT_SIZE]; 329 spx_word16_t _out[MAX_FFT_SIZE]; 330 #endif 331 for (i=0;i<N;i++) 332 _in[i] = (int)floor(.5+in[i]); 333 spx_fft(table, _in, _out); 334 for (i=0;i<N;i++) 335 out[i] = _out[i]; 336 #if 0 337 if (!fixed_point) 338 { 339 float scale; 340 struct drft_lookup t; 341 spx_drft_init(&t, ((struct kiss_config *)table)->N); 342 scale = 1./((struct kiss_config *)table)->N; 343 for (i=0;i<((struct kiss_config *)table)->N;i++) 344 out[i] = scale*in[i]; 345 spx_drft_forward(&t, out); 346 spx_drft_clear(&t); 347 } 348 #endif 349 } 350 351 void spx_ifft_float(void *table, float *in, float *out) 352 { 353 int i; 354 #ifdef USE_SMALLFT 355 int N = ((struct drft_lookup *)table)->n; 356 #elif defined(USE_KISS_FFT) 357 int N = ((struct kiss_config *)table)->N; 358 #else 359 #endif 360 #ifdef VAR_ARRAYS 361 spx_word16_t _in[N]; 362 spx_word16_t _out[N]; 363 #else 364 spx_word16_t _in[MAX_FFT_SIZE]; 365 spx_word16_t _out[MAX_FFT_SIZE]; 366 #endif 367 for (i=0;i<N;i++) 368 _in[i] = (int)floor(.5+in[i]); 369 spx_ifft(table, _in, _out); 370 for (i=0;i<N;i++) 371 out[i] = _out[i]; 372 #if 0 373 if (!fixed_point) 374 { 375 int i; 376 struct drft_lookup t; 377 spx_drft_init(&t, ((struct kiss_config *)table)->N); 378 for (i=0;i<((struct kiss_config *)table)->N;i++) 379 out[i] = in[i]; 380 spx_drft_backward(&t, out); 381 spx_drft_clear(&t); 382 } 383 #endif 384 } 385 386 #else 387 388 void spx_fft_float(void *table, float *in, float *out) 389 { 390 spx_fft(table, in, out); 391 } 392 void spx_ifft_float(void *table, float *in, float *out) 393 { 394 spx_ifft(table, in, out); 395 } 396 397 #endif 398
