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_INTEL_IPP)
169
170#include <ipps.h>
171
172struct ipp_fft_config
173{
174  IppsDFTSpec_R_32f *dftSpec;
175  Ipp8u *buffer;
176};
177
178void *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
199void 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
207void 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
213void 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
223struct fftw_config {
224  float *in;
225  float *out;
226  fftwf_plan fft;
227  fftwf_plan ifft;
228  int N;
229};
230
231void *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
244void 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
255void 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
273void 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
298struct kiss_config {
299   kiss_fftr_cfg forward;
300   kiss_fftr_cfg backward;
301   int N;
302};
303
304void *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
314void 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
324void 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
336void 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
348void 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
366void 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
402void 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
439void spx_fft_float(void *table, float *in, float *out)
440{
441   spx_fft(table, in, out);
442}
443void spx_ifft_float(void *table, float *in, float *out)
444{
445   spx_ifft(table, in, out);
446}
447
448#endif
449