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