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
maximize_range(spx_word16_t * in,spx_word16_t * out,spx_word16_t bound,int len)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 
renorm_range(spx_word16_t * in,spx_word16_t * out,int shift,int len)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 
spx_fft_init(int size)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 
spx_fft_destroy(void * table)92 void spx_fft_destroy(void *table)
93 {
94    spx_drft_clear(table);
95    speex_free(table);
96 }
97 
spx_fft(void * table,float * in,float * out)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 
spx_ifft(void * table,float * in,float * out)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 
spx_fft_init(int size)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 
spx_fft_destroy(void * table)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 
spx_fft(void * table,spx_word16_t * in,spx_word16_t * out)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 
spx_ifft(void * table,spx_word16_t * in,spx_word16_t * out)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_GPL_FFTW3)
169 
170 #include <fftw3.h>
171 
172 struct fftw_config {
173   float *in;
174   float *out;
175   fftwf_plan fft;
176   fftwf_plan ifft;
177   int N;
178 };
179 
spx_fft_init(int size)180 void *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 
spx_fft_destroy(void * table)193 void 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 
spx_fft(void * table,spx_word16_t * in,spx_word16_t * out)204 void 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 
spx_ifft(void * table,spx_word16_t * in,spx_word16_t * out)222 void 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 
247 struct kiss_config {
248    kiss_fftr_cfg forward;
249    kiss_fftr_cfg backward;
250    int N;
251 };
252 
spx_fft_init(int size)253 void *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 
spx_fft_destroy(void * table)263 void 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 
spx_fft(void * table,spx_word16_t * in,spx_word16_t * out)273 void 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 
spx_fft(void * table,spx_word16_t * in,spx_word16_t * out)285 void 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 
spx_ifft(void * table,spx_word16_t * in,spx_word16_t * out)297 void 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 
spx_fft_float(void * table,float * in,float * out)315 void 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 
spx_ifft_float(void * table,float * in,float * out)351 void 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 
spx_fft_float(void * table,float * in,float * out)388 void spx_fft_float(void *table, float *in, float *out)
389 {
390    spx_fft(table, in, out);
391 }
spx_ifft_float(void * table,float * in,float * out)392 void spx_ifft_float(void *table, float *in, float *out)
393 {
394    spx_ifft(table, in, out);
395 }
396 
397 #endif
398