Home | History | Annotate | Download | only in libspeexdsp
      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