fftwrap.c revision 98913fed6520d8849fb2e246be943e04474aefa4 
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
45static 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
69static 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
84void *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
92void spx_fft_destroy(void *table)
93{
94   spx_drft_clear(table);
95   speex_free(table);
96}
97
98void 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
116void 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
132struct mkl_config {
133  DFTI_DESCRIPTOR_HANDLE desc;
134  int N;
135};
136
137void *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
149void 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
156void 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
162void 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
172struct fftw_config {
173  float *in;
174  float *out;
175  fftwf_plan fft;
176  fftwf_plan ifft;
177  int N;
178};
179
180void *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
193void 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
204void 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
222void 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
247struct kiss_config {
248   kiss_fftr_cfg forward;
249   kiss_fftr_cfg backward;
250   int N;
251};
252
253void *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
263void 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
273void 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
285void 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
297void 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
315void 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
351void 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
388void spx_fft_float(void *table, float *in, float *out)
389{
390   spx_fft(table, in, out);
391}
392void spx_ifft_float(void *table, float *in, float *out)
393{
394   spx_ifft(table, in, out);
395}
396
397#endif
398