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