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