fftwrap.c 10.4 KB
Newer Older
Tristan Matthews's avatar
Tristan Matthews committed
1
/* Copyright (C) 2005-2006 Jean-Marc Valin
2 3
   File: fftwrap.c

Tristan Matthews's avatar
Tristan Matthews committed
4
   Wrapper for various FFTs
5 6 7 8

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
Tristan Matthews's avatar
Tristan Matthews committed
9

10 11
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
Tristan Matthews's avatar
Tristan Matthews committed
12

13 14 15
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
Tristan Matthews's avatar
Tristan Matthews committed
16

17 18 19
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
Tristan Matthews's avatar
Tristan Matthews committed
20

21 22 23 24 25 26 27 28 29 30 31 32 33 34
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

35 36 37 38
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

jm's avatar
jm committed
39
#include "arch.h"
40
#include "os_support.h"
jm's avatar
jm committed
41

42
#define MAX_FFT_SIZE 2048
jm's avatar
jm committed
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

#ifdef FIXED_POINT
static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len)
{
   int i, shift;
   spx_word16_t max_val = 0;
   for (i=0;i<len;i++)
   {
      if (in[i]>max_val)
         max_val = in[i];
      if (-in[i]>max_val)
         max_val = -in[i];
   }
   shift=0;
   while (max_val <= (bound>>1) && max_val != 0)
   {
      max_val <<= 1;
      shift++;
   }
   for (i=0;i<len;i++)
   {
64
      out[i] = SHL16(in[i], shift);
Tristan Matthews's avatar
Tristan Matthews committed
65
   }
jm's avatar
jm committed
66 67 68 69 70 71 72 73
   return shift;
}

static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len)
{
   int i;
   for (i=0;i<len;i++)
   {
74
      out[i] = PSHR16(in[i], shift);
jm's avatar
jm committed
75 76 77
   }
}
#endif
78 79 80 81

#ifdef USE_SMALLFT

#include "smallft.h"
jm's avatar
jm committed
82
#include <math.h>
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

void *spx_fft_init(int size)
{
   struct drft_lookup *table;
   table = speex_alloc(sizeof(struct drft_lookup));
   spx_drft_init((struct drft_lookup *)table, size);
   return (void*)table;
}

void spx_fft_destroy(void *table)
{
   spx_drft_clear(table);
   speex_free(table);
}

void spx_fft(void *table, float *in, float *out)
{
   if (in==out)
   {
      int i;
      float scale = 1./((struct drft_lookup *)table)->n;
jm's avatar
jm committed
104
      speex_warning("FFT should not be done in-place");
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
      for (i=0;i<((struct drft_lookup *)table)->n;i++)
         out[i] = scale*in[i];
   } else {
      int i;
      float scale = 1./((struct drft_lookup *)table)->n;
      for (i=0;i<((struct drft_lookup *)table)->n;i++)
         out[i] = scale*in[i];
   }
   spx_drft_forward((struct drft_lookup *)table, out);
}

void spx_ifft(void *table, float *in, float *out)
{
   if (in==out)
   {
      speex_warning("FFT should not be done in-place");
   } else {
      int i;
      for (i=0;i<((struct drft_lookup *)table)->n;i++)
         out[i] = in[i];
   }
   spx_drft_backward((struct drft_lookup *)table, out);
}

129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
#elif defined(USE_INTEL_MKL)
#include <mkl.h>

struct mkl_config {
  DFTI_DESCRIPTOR_HANDLE desc;
  int N;
};

void *spx_fft_init(int size)
{
  struct mkl_config *table = (struct mkl_config *) speex_alloc(sizeof(struct mkl_config));
  table->N = size;
  DftiCreateDescriptor(&table->desc, DFTI_SINGLE, DFTI_REAL, 1, size);
  DftiSetValue(table->desc, DFTI_PACKED_FORMAT, DFTI_PACK_FORMAT);
  DftiSetValue(table->desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
  DftiSetValue(table->desc, DFTI_FORWARD_SCALE, 1.0f / size);
  DftiCommitDescriptor(table->desc);
  return table;
}

void spx_fft_destroy(void *table)
{
  struct mkl_config *t = (struct mkl_config *) table;
  DftiFreeDescriptor(t->desc);
  speex_free(table);
}

void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
{
  struct mkl_config *t = (struct mkl_config *) table;
  DftiComputeForward(t->desc, in, out);
}

void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
{
  struct mkl_config *t = (struct mkl_config *) table;
  DftiComputeBackward(t->desc, in, out);
}

Jeff Wallace's avatar
Jeff Wallace committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
#elif defined(USE_INTEL_IPP)

#include <ipps.h>

struct ipp_fft_config
{
  IppsDFTSpec_R_32f *dftSpec;
  Ipp8u *buffer;
};

void *spx_fft_init(int size)
{
  int bufferSize = 0;
  int hint;
  struct ipp_fft_config *table;

  table = (struct ipp_fft_config *)speex_alloc(sizeof(struct ipp_fft_config));

  /* there appears to be no performance difference between ippAlgHintFast and
     ippAlgHintAccurate when using the with the floating point version
     of the fft. */
  hint = ippAlgHintAccurate;

  ippsDFTInitAlloc_R_32f(&table->dftSpec, size, IPP_FFT_DIV_FWD_BY_N, hint);

  ippsDFTGetBufSize_R_32f(table->dftSpec, &bufferSize);
  table->buffer = ippsMalloc_8u(bufferSize);

  return table;
}

void spx_fft_destroy(void *table)
{
  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
  ippsFree(t->buffer);
  ippsDFTFree_R_32f(t->dftSpec);
  speex_free(t);
}

void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
{
  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
  ippsDFTFwd_RToPack_32f(in, out, t->dftSpec, t->buffer);
}

void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
{
  struct ipp_fft_config *t = (struct ipp_fft_config *)table;
  ippsDFTInv_PackToR_32f(in, out, t->dftSpec, t->buffer);
}

219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
#elif defined(USE_GPL_FFTW3)

#include <fftw3.h>

struct fftw_config {
  float *in;
  float *out;
  fftwf_plan fft;
  fftwf_plan ifft;
  int N;
};

void *spx_fft_init(int size)
{
  struct fftw_config *table = (struct fftw_config *) speex_alloc(sizeof(struct fftw_config));
  table->in = fftwf_malloc(sizeof(float) * (size+2));
  table->out = fftwf_malloc(sizeof(float) * (size+2));

  table->fft = fftwf_plan_dft_r2c_1d(size, table->in, (fftwf_complex *) table->out, FFTW_PATIENT);
  table->ifft = fftwf_plan_dft_c2r_1d(size, (fftwf_complex *) table->in, table->out, FFTW_PATIENT);

  table->N = size;
  return table;
}

void spx_fft_destroy(void *table)
{
  struct fftw_config *t = (struct fftw_config *) table;
  fftwf_destroy_plan(t->fft);
  fftwf_destroy_plan(t->ifft);
  fftwf_free(t->in);
  fftwf_free(t->out);
  speex_free(table);
}


void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
{
  int i;
  struct fftw_config *t = (struct fftw_config *) table;
  const int N = t->N;
  float *iptr = t->in;
  float *optr = t->out;
  const float m = 1.0 / N;
  for(i=0;i<N;++i)
    iptr[i]=in[i] * m;

  fftwf_execute(t->fft);

  out[0] = optr[0];
  for(i=1;i<N;++i)
    out[i] = optr[i+1];
}

Tristan Matthews's avatar
Tristan Matthews committed
273
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
274 275 276 277 278 279 280 281 282 283 284
{
  int i;
  struct fftw_config *t = (struct fftw_config *) table;
  const int N = t->N;
  float *iptr = t->in;
  float *optr = t->out;

  iptr[0] = in[0];
  iptr[1] = 0.0f;
  for(i=1;i<N;++i)
    iptr[i+1] = in[i];
jm's avatar
jm committed
285
  iptr[N+1] = 0.0f;
286 287

  fftwf_execute(t->ifft);
Tristan Matthews's avatar
Tristan Matthews committed
288

289 290 291 292
  for(i=0;i<N;++i)
    out[i] = optr[i];
}

293 294 295 296 297 298 299 300 301 302 303 304 305 306
#elif defined(USE_KISS_FFT)

#include "kiss_fftr.h"
#include "kiss_fft.h"

struct kiss_config {
   kiss_fftr_cfg forward;
   kiss_fftr_cfg backward;
   int N;
};

void *spx_fft_init(int size)
{
   struct kiss_config *table;
jm's avatar
jm committed
307
   table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config));
308 309 310 311 312 313 314 315 316 317 318 319 320 321
   table->forward = kiss_fftr_alloc(size,0,NULL,NULL);
   table->backward = kiss_fftr_alloc(size,1,NULL,NULL);
   table->N = size;
   return table;
}

void spx_fft_destroy(void *table)
{
   struct kiss_config *t = (struct kiss_config *)table;
   kiss_fftr_free(t->forward);
   kiss_fftr_free(t->backward);
   speex_free(table);
}

jm's avatar
jm committed
322 323
#ifdef FIXED_POINT

324 325 326 327
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
{
   int shift;
   struct kiss_config *t = (struct kiss_config *)table;
jm's avatar
jm committed
328
   shift = maximize_range(in, in, 32000, t->N);
329
   kiss_fftr2(t->forward, in, out);
jm's avatar
jm committed
330 331 332 333
   renorm_range(in, in, shift, t->N);
   renorm_range(out, out, shift, t->N);
}

334
#else
jm's avatar
jm committed
335 336 337 338 339 340

void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
{
   int i;
   float scale;
   struct kiss_config *t = (struct kiss_config *)table;
341
   scale = 1./t->N;
342 343 344
   kiss_fftr2(t->forward, in, out);
   for (i=0;i<t->N;i++)
      out[i] *= scale;
345
}
jm's avatar
jm committed
346
#endif
347 348 349 350

void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
{
   struct kiss_config *t = (struct kiss_config *)table;
351
   kiss_fftri2(t->backward, in, out);
352 353 354
}


355 356 357 358 359 360 361
#else

#error No other FFT implemented

#endif


362
#ifdef FIXED_POINT
jm's avatar
jm committed
363
/*#include "smallft.h"*/
364 365 366 367 368 369 370 371 372 373 374


void spx_fft_float(void *table, float *in, float *out)
{
   int i;
#ifdef USE_SMALLFT
   int N = ((struct drft_lookup *)table)->n;
#elif defined(USE_KISS_FFT)
   int N = ((struct kiss_config *)table)->N;
#else
#endif
375
#ifdef VAR_ARRAYS
376 377
   spx_word16_t _in[N];
   spx_word16_t _out[N];
378 379 380 381
#else
   spx_word16_t _in[MAX_FFT_SIZE];
   spx_word16_t _out[MAX_FFT_SIZE];
#endif
382 383 384 385 386
   for (i=0;i<N;i++)
      _in[i] = (int)floor(.5+in[i]);
   spx_fft(table, _in, _out);
   for (i=0;i<N;i++)
      out[i] = _out[i];
jm's avatar
jm committed
387
#if 0
jm's avatar
jm committed
388 389
   if (!fixed_point)
   {
390
      float scale;
jm's avatar
jm committed
391 392
      struct drft_lookup t;
      spx_drft_init(&t, ((struct kiss_config *)table)->N);
393
      scale = 1./((struct kiss_config *)table)->N;
jm's avatar
jm committed
394 395 396 397 398
      for (i=0;i<((struct kiss_config *)table)->N;i++)
         out[i] = scale*in[i];
      spx_drft_forward(&t, out);
      spx_drft_clear(&t);
   }
jm's avatar
jm committed
399
#endif
400 401 402 403 404 405 406 407 408 409 410
}

void spx_ifft_float(void *table, float *in, float *out)
{
   int i;
#ifdef USE_SMALLFT
   int N = ((struct drft_lookup *)table)->n;
#elif defined(USE_KISS_FFT)
   int N = ((struct kiss_config *)table)->N;
#else
#endif
411
#ifdef VAR_ARRAYS
412 413
   spx_word16_t _in[N];
   spx_word16_t _out[N];
414 415 416 417
#else
   spx_word16_t _in[MAX_FFT_SIZE];
   spx_word16_t _out[MAX_FFT_SIZE];
#endif
418 419 420 421 422
   for (i=0;i<N;i++)
      _in[i] = (int)floor(.5+in[i]);
   spx_ifft(table, _in, _out);
   for (i=0;i<N;i++)
      out[i] = _out[i];
jm's avatar
jm committed
423
#if 0
jm's avatar
jm committed
424 425 426 427 428 429 430 431 432 433
   if (!fixed_point)
   {
      int i;
      struct drft_lookup t;
      spx_drft_init(&t, ((struct kiss_config *)table)->N);
      for (i=0;i<((struct kiss_config *)table)->N;i++)
         out[i] = in[i];
      spx_drft_backward(&t, out);
      spx_drft_clear(&t);
   }
jm's avatar
jm committed
434
#endif
435 436 437
}

#else
438 439 440 441 442 443 444 445 446

void spx_fft_float(void *table, float *in, float *out)
{
   spx_fft(table, in, out);
}
void spx_ifft_float(void *table, float *in, float *out)
{
   spx_ifft(table, in, out);
}
447 448

#endif