libghost/fftwrap.c

   1 /* Copyright (C) 2005 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 #define USE_SMALLFT
  40 /*#define USE_KISS_FFT*/
  41
  42
  43 #include "misc.h"
  44
  45
  46 #ifdef FIXED_POINT
  47 static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len)
  48 {
  49    int i, shift;
  50    spx_word16_t max_val = 0;
  51    for (i=0;i<len;i++)
  52    {
  53       if (in[i]>max_val)
  54          max_val = in[i];
  55       if (-in[i]>max_val)
  56          max_val = -in[i];
  57    }
  58    shift=0;
  59    while (max_val <= (bound>>1) && max_val != 0)
  60    {
  61       max_val <<= 1;
  62       shift++;
  63    }
  64    for (i=0;i<len;i++)
  65    {
  66       out[i] = in[i] << shift;
  67    }
  68    return shift;
  69 }
  70
  71 static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len)
  72 {
  73    int i;
  74    for (i=0;i<len;i++)
  75    {
  76       out[i] = (in[i] + (1<<(shift-1))) >> shift;
  77    }
  78 }
  79 #endif
  80
  81 #ifdef USE_SMALLFT
  82
  83 #include "smallft.h"
  84 #include <math.h>
  85
  86 void *spx_fft_init(int size)
  87 {
  88    struct drft_lookup *table;
  89    table = speex_alloc(sizeof(struct drft_lookup));
  90    spx_drft_init((struct drft_lookup *)table, size);
  91    return (void*)table;
  92 }
  93
  94 void spx_fft_destroy(void *table)
  95 {
  96    spx_drft_clear(table);
  97    speex_free(table);
  98 }
  99
 100 void spx_fft(void *table, float *in, float *out)
 101 {
 102    if (in==out)
 103    {
 104       int i;
 105       speex_warning("FFT should not be done in-place");
 106       float scale = 1./((struct drft_lookup *)table)->n;
 107       for (i=0;i<((struct drft_lookup *)table)->n;i++)
 108          out[i] = scale*in[i];
 109    } else {
 110       int i;
 111       float scale = 1./((struct drft_lookup *)table)->n;
 112       for (i=0;i<((struct drft_lookup *)table)->n;i++)
 113          out[i] = scale*in[i];
 114    }
 115    spx_drft_forward((struct drft_lookup *)table, out);
 116 }
 117
 118 void spx_ifft(void *table, float *in, float *out)
 119 {
 120    if (in==out)
 121    {
 122       int i;
 123       speex_warning("FFT should not be done in-place");
 124    } else {
 125       int i;
 126       for (i=0;i<((struct drft_lookup *)table)->n;i++)
 127          out[i] = in[i];
 128    }
 129    spx_drft_backward((struct drft_lookup *)table, out);
 130 }
 131
 132 #elif defined(USE_KISS_FFT)
 133
 134 #include "kiss_fftr.h"
 135 #include "kiss_fft.h"
 136
 137 struct kiss_config {
 138    kiss_fftr_cfg forward;
 139    kiss_fftr_cfg backward;
 140    kiss_fft_cpx *freq_data;
 141    int N;
 142 };
 143
 144 void *spx_fft_init(int size)
 145 {
 146    struct kiss_config *table;
 147    table = speex_alloc(sizeof(struct kiss_config));
 148    table->freq_data = speex_alloc(sizeof(kiss_fft_cpx)*((size>>1)+1));
 149    table->forward = kiss_fftr_alloc(size,0,NULL,NULL);
 150    table->backward = kiss_fftr_alloc(size,1,NULL,NULL);
 151    table->N = size;
 152    return table;
 153 }
 154
 155 void spx_fft_destroy(void *table)
 156 {
 157    struct kiss_config *t = (struct kiss_config *)table;
 158    kiss_fftr_free(t->forward);
 159    kiss_fftr_free(t->backward);
 160    speex_free(t->freq_data);
 161    speex_free(table);
 162 }
 163
 164 #ifdef FIXED_POINT
 165
 166 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
 167 {
 168    int i;
 169    int shift;
 170    struct kiss_config *t = (struct kiss_config *)table;
 171    shift = maximize_range(in, in, 32000, t->N);
 172    kiss_fftr(t->forward, in, t->freq_data);
 173    out[0] = t->freq_data[0].r;
 174    for (i=1;i<t->N>>1;i++)
 175    {
 176       out[(i<<1)-1] = t->freq_data[i].r;
 177       out[(i<<1)] = t->freq_data[i].i;
 178    }
 179    out[(i<<1)-1] = t->freq_data[i].r;
 180    renorm_range(in, in, shift, t->N);
 181    renorm_range(out, out, shift, t->N);
 182 }
 183
 184 #else
 185
 186 void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
 187 {
 188    int i;
 189    float scale;
 190    struct kiss_config *t = (struct kiss_config *)table;
 191    scale = 1./t->N;
 192    kiss_fftr(t->forward, in, t->freq_data);
 193    out[0] = scale*t->freq_data[0].r;
 194    for (i=1;i<t->N>>1;i++)
 195    {
 196       out[(i<<1)-1] = scale*t->freq_data[i].r;
 197       out[(i<<1)] = scale*t->freq_data[i].i;
 198    }
 199    out[(i<<1)-1] = scale*t->freq_data[i].r;
 200 }
 201 #endif
 202
 203 void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
 204 {
 205    int i;
 206    struct kiss_config *t = (struct kiss_config *)table;
 207    t->freq_data[0].r = in[0];
 208    t->freq_data[0].i = 0;
 209    for (i=1;i<t->N>>1;i++)
 210    {
 211       t->freq_data[i].r = in[(i<<1)-1];
 212       t->freq_data[i].i = in[(i<<1)];
 213    }
 214    t->freq_data[i].r = in[(i<<1)-1];
 215    t->freq_data[i].i = 0;
 216
 217    kiss_fftri(t->backward, t->freq_data, out);
 218 }
 219
 220
 221 #else
 222
 223 #error No other FFT implemented
 224
 225 #endif
 226
 227
 228 int fixed_point = 1;
 229 #ifdef FIXED_POINT
 230 #include "smallft.h"
 231
 232
 233 void spx_fft_float(void *table, float *in, float *out)
 234 {
 235    int i;
 236 #ifdef USE_SMALLFT
 237    int N = ((struct drft_lookup *)table)->n;
 238 #elif defined(USE_KISS_FFT)
 239    int N = ((struct kiss_config *)table)->N;
 240 #else
 241 #endif
 242    spx_word16_t _in[N];
 243    spx_word16_t _out[N];
 244    for (i=0;i<N;i++)
 245       _in[i] = (int)floor(.5+in[i]);
 246    spx_fft(table, _in, _out);
 247    for (i=0;i<N;i++)
 248       out[i] = _out[i];
 249    if (!fixed_point)
 250    {
 251       struct drft_lookup t;
 252       spx_drft_init(&t, ((struct kiss_config *)table)->N);
 253       float scale = 1./((struct kiss_config *)table)->N;
 254       for (i=0;i<((struct kiss_config *)table)->N;i++)
 255          out[i] = scale*in[i];
 256       spx_drft_forward(&t, out);
 257       spx_drft_clear(&t);
 258    }
 259 }
 260
 261 void spx_ifft_float(void *table, float *in, float *out)
 262 {
 263    int i;
 264 #ifdef USE_SMALLFT
 265    int N = ((struct drft_lookup *)table)->n;
 266 #elif defined(USE_KISS_FFT)
 267    int N = ((struct kiss_config *)table)->N;
 268 #else
 269 #endif
 270    spx_word16_t _in[N];
 271    spx_word16_t _out[N];
 272    for (i=0;i<N;i++)
 273       _in[i] = (int)floor(.5+in[i]);
 274    spx_ifft(table, _in, _out);
 275    for (i=0;i<N;i++)
 276       out[i] = _out[i];
 277    if (!fixed_point)
 278    {
 279       int i;
 280       struct drft_lookup t;
 281       spx_drft_init(&t, ((struct kiss_config *)table)->N);
 282       for (i=0;i<((struct kiss_config *)table)->N;i++)
 283          out[i] = in[i];
 284       spx_drft_backward(&t, out);
 285       spx_drft_clear(&t);
 286    }
 287 }
 288
 289 #else
 290
 291 void spx_fft_float(void *table, float *in, float *out)
 292 {
 293    spx_fft(table, in, out);
 294 }
 295 void spx_ifft_float(void *table, float *in, float *out)
 296 {
 297    spx_ifft(table, in, out);
 298 }
 299
 300 #endif