]> 4ch.mooo.com Git - 16.git/blob - src/lib/dl/ext/speex/fftwrap.c
meh did some cleanings and i will work on mapread to mm thingy sometime soon! oops...
[16.git] / src / lib / dl / ext / speex / fftwrap.c
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_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
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
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
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
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
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
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
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
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
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
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
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
388 void spx_fft_float(void *table, float *in, float *out)
389 {
390    spx_fft(table, in, out);
391 }
392 void spx_ifft_float(void *table, float *in, float *out)
393 {
394    spx_ifft(table, in, out);
395 }
396
397 #endif