src/lib/dl/ext/speex/fftwrap.c

   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