|
1 /*Copyright (c) 2003-2004, Mark Borgerding |
|
2 Lots of modifications by Jean-Marc Valin |
|
3 Copyright (c) 2005-2007, Xiph.Org Foundation |
|
4 Copyright (c) 2008, Xiph.Org Foundation, CSIRO |
|
5 |
|
6 All rights reserved. |
|
7 |
|
8 Redistribution and use in source and binary forms, with or without |
|
9 modification, are permitted provided that the following conditions are met: |
|
10 |
|
11 * Redistributions of source code must retain the above copyright notice, |
|
12 this list of conditions and the following disclaimer. |
|
13 * Redistributions in binary form must reproduce the above copyright notice, |
|
14 this list of conditions and the following disclaimer in the |
|
15 documentation and/or other materials provided with the distribution. |
|
16 |
|
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
|
18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
|
19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
|
20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
|
21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
|
22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
|
23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
|
24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
|
25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
|
26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
|
27 POSSIBILITY OF SUCH DAMAGE.*/ |
|
28 |
|
29 /* This code is originally from Mark Borgerding's KISS-FFT but has been |
|
30 heavily modified to better suit Opus */ |
|
31 |
|
32 #ifndef SKIP_CONFIG_H |
|
33 # ifdef HAVE_CONFIG_H |
|
34 # include "config.h" |
|
35 # endif |
|
36 #endif |
|
37 |
|
38 #include "_kiss_fft_guts.h" |
|
39 #include "arch.h" |
|
40 #include "os_support.h" |
|
41 #include "mathops.h" |
|
42 #include "stack_alloc.h" |
|
43 |
|
44 /* The guts header contains all the multiplication and addition macros that are defined for |
|
45 complex numbers. It also delares the kf_ internal functions. |
|
46 */ |
|
47 |
|
48 static void kf_bfly2( |
|
49 kiss_fft_cpx * Fout, |
|
50 const size_t fstride, |
|
51 const kiss_fft_state *st, |
|
52 int m, |
|
53 int N, |
|
54 int mm |
|
55 ) |
|
56 { |
|
57 kiss_fft_cpx * Fout2; |
|
58 const kiss_twiddle_cpx * tw1; |
|
59 int i,j; |
|
60 kiss_fft_cpx * Fout_beg = Fout; |
|
61 for (i=0;i<N;i++) |
|
62 { |
|
63 Fout = Fout_beg + i*mm; |
|
64 Fout2 = Fout + m; |
|
65 tw1 = st->twiddles; |
|
66 for(j=0;j<m;j++) |
|
67 { |
|
68 kiss_fft_cpx t; |
|
69 Fout->r = SHR32(Fout->r, 1);Fout->i = SHR32(Fout->i, 1); |
|
70 Fout2->r = SHR32(Fout2->r, 1);Fout2->i = SHR32(Fout2->i, 1); |
|
71 C_MUL (t, *Fout2 , *tw1); |
|
72 tw1 += fstride; |
|
73 C_SUB( *Fout2 , *Fout , t ); |
|
74 C_ADDTO( *Fout , t ); |
|
75 ++Fout2; |
|
76 ++Fout; |
|
77 } |
|
78 } |
|
79 } |
|
80 |
|
81 static void ki_bfly2( |
|
82 kiss_fft_cpx * Fout, |
|
83 const size_t fstride, |
|
84 const kiss_fft_state *st, |
|
85 int m, |
|
86 int N, |
|
87 int mm |
|
88 ) |
|
89 { |
|
90 kiss_fft_cpx * Fout2; |
|
91 const kiss_twiddle_cpx * tw1; |
|
92 kiss_fft_cpx t; |
|
93 int i,j; |
|
94 kiss_fft_cpx * Fout_beg = Fout; |
|
95 for (i=0;i<N;i++) |
|
96 { |
|
97 Fout = Fout_beg + i*mm; |
|
98 Fout2 = Fout + m; |
|
99 tw1 = st->twiddles; |
|
100 for(j=0;j<m;j++) |
|
101 { |
|
102 C_MULC (t, *Fout2 , *tw1); |
|
103 tw1 += fstride; |
|
104 C_SUB( *Fout2 , *Fout , t ); |
|
105 C_ADDTO( *Fout , t ); |
|
106 ++Fout2; |
|
107 ++Fout; |
|
108 } |
|
109 } |
|
110 } |
|
111 |
|
112 static void kf_bfly4( |
|
113 kiss_fft_cpx * Fout, |
|
114 const size_t fstride, |
|
115 const kiss_fft_state *st, |
|
116 int m, |
|
117 int N, |
|
118 int mm |
|
119 ) |
|
120 { |
|
121 const kiss_twiddle_cpx *tw1,*tw2,*tw3; |
|
122 kiss_fft_cpx scratch[6]; |
|
123 const size_t m2=2*m; |
|
124 const size_t m3=3*m; |
|
125 int i, j; |
|
126 |
|
127 kiss_fft_cpx * Fout_beg = Fout; |
|
128 for (i=0;i<N;i++) |
|
129 { |
|
130 Fout = Fout_beg + i*mm; |
|
131 tw3 = tw2 = tw1 = st->twiddles; |
|
132 for (j=0;j<m;j++) |
|
133 { |
|
134 C_MUL4(scratch[0],Fout[m] , *tw1 ); |
|
135 C_MUL4(scratch[1],Fout[m2] , *tw2 ); |
|
136 C_MUL4(scratch[2],Fout[m3] , *tw3 ); |
|
137 |
|
138 Fout->r = PSHR32(Fout->r, 2); |
|
139 Fout->i = PSHR32(Fout->i, 2); |
|
140 C_SUB( scratch[5] , *Fout, scratch[1] ); |
|
141 C_ADDTO(*Fout, scratch[1]); |
|
142 C_ADD( scratch[3] , scratch[0] , scratch[2] ); |
|
143 C_SUB( scratch[4] , scratch[0] , scratch[2] ); |
|
144 C_SUB( Fout[m2], *Fout, scratch[3] ); |
|
145 tw1 += fstride; |
|
146 tw2 += fstride*2; |
|
147 tw3 += fstride*3; |
|
148 C_ADDTO( *Fout , scratch[3] ); |
|
149 |
|
150 Fout[m].r = scratch[5].r + scratch[4].i; |
|
151 Fout[m].i = scratch[5].i - scratch[4].r; |
|
152 Fout[m3].r = scratch[5].r - scratch[4].i; |
|
153 Fout[m3].i = scratch[5].i + scratch[4].r; |
|
154 ++Fout; |
|
155 } |
|
156 } |
|
157 } |
|
158 |
|
159 static void ki_bfly4( |
|
160 kiss_fft_cpx * Fout, |
|
161 const size_t fstride, |
|
162 const kiss_fft_state *st, |
|
163 int m, |
|
164 int N, |
|
165 int mm |
|
166 ) |
|
167 { |
|
168 const kiss_twiddle_cpx *tw1,*tw2,*tw3; |
|
169 kiss_fft_cpx scratch[6]; |
|
170 const size_t m2=2*m; |
|
171 const size_t m3=3*m; |
|
172 int i, j; |
|
173 |
|
174 kiss_fft_cpx * Fout_beg = Fout; |
|
175 for (i=0;i<N;i++) |
|
176 { |
|
177 Fout = Fout_beg + i*mm; |
|
178 tw3 = tw2 = tw1 = st->twiddles; |
|
179 for (j=0;j<m;j++) |
|
180 { |
|
181 C_MULC(scratch[0],Fout[m] , *tw1 ); |
|
182 C_MULC(scratch[1],Fout[m2] , *tw2 ); |
|
183 C_MULC(scratch[2],Fout[m3] , *tw3 ); |
|
184 |
|
185 C_SUB( scratch[5] , *Fout, scratch[1] ); |
|
186 C_ADDTO(*Fout, scratch[1]); |
|
187 C_ADD( scratch[3] , scratch[0] , scratch[2] ); |
|
188 C_SUB( scratch[4] , scratch[0] , scratch[2] ); |
|
189 C_SUB( Fout[m2], *Fout, scratch[3] ); |
|
190 tw1 += fstride; |
|
191 tw2 += fstride*2; |
|
192 tw3 += fstride*3; |
|
193 C_ADDTO( *Fout , scratch[3] ); |
|
194 |
|
195 Fout[m].r = scratch[5].r - scratch[4].i; |
|
196 Fout[m].i = scratch[5].i + scratch[4].r; |
|
197 Fout[m3].r = scratch[5].r + scratch[4].i; |
|
198 Fout[m3].i = scratch[5].i - scratch[4].r; |
|
199 ++Fout; |
|
200 } |
|
201 } |
|
202 } |
|
203 |
|
204 #ifndef RADIX_TWO_ONLY |
|
205 |
|
206 static void kf_bfly3( |
|
207 kiss_fft_cpx * Fout, |
|
208 const size_t fstride, |
|
209 const kiss_fft_state *st, |
|
210 int m, |
|
211 int N, |
|
212 int mm |
|
213 ) |
|
214 { |
|
215 int i; |
|
216 size_t k; |
|
217 const size_t m2 = 2*m; |
|
218 const kiss_twiddle_cpx *tw1,*tw2; |
|
219 kiss_fft_cpx scratch[5]; |
|
220 kiss_twiddle_cpx epi3; |
|
221 |
|
222 kiss_fft_cpx * Fout_beg = Fout; |
|
223 epi3 = st->twiddles[fstride*m]; |
|
224 for (i=0;i<N;i++) |
|
225 { |
|
226 Fout = Fout_beg + i*mm; |
|
227 tw1=tw2=st->twiddles; |
|
228 k=m; |
|
229 do { |
|
230 C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); |
|
231 |
|
232 C_MUL(scratch[1],Fout[m] , *tw1); |
|
233 C_MUL(scratch[2],Fout[m2] , *tw2); |
|
234 |
|
235 C_ADD(scratch[3],scratch[1],scratch[2]); |
|
236 C_SUB(scratch[0],scratch[1],scratch[2]); |
|
237 tw1 += fstride; |
|
238 tw2 += fstride*2; |
|
239 |
|
240 Fout[m].r = Fout->r - HALF_OF(scratch[3].r); |
|
241 Fout[m].i = Fout->i - HALF_OF(scratch[3].i); |
|
242 |
|
243 C_MULBYSCALAR( scratch[0] , epi3.i ); |
|
244 |
|
245 C_ADDTO(*Fout,scratch[3]); |
|
246 |
|
247 Fout[m2].r = Fout[m].r + scratch[0].i; |
|
248 Fout[m2].i = Fout[m].i - scratch[0].r; |
|
249 |
|
250 Fout[m].r -= scratch[0].i; |
|
251 Fout[m].i += scratch[0].r; |
|
252 |
|
253 ++Fout; |
|
254 } while(--k); |
|
255 } |
|
256 } |
|
257 |
|
258 static void ki_bfly3( |
|
259 kiss_fft_cpx * Fout, |
|
260 const size_t fstride, |
|
261 const kiss_fft_state *st, |
|
262 int m, |
|
263 int N, |
|
264 int mm |
|
265 ) |
|
266 { |
|
267 int i, k; |
|
268 const size_t m2 = 2*m; |
|
269 const kiss_twiddle_cpx *tw1,*tw2; |
|
270 kiss_fft_cpx scratch[5]; |
|
271 kiss_twiddle_cpx epi3; |
|
272 |
|
273 kiss_fft_cpx * Fout_beg = Fout; |
|
274 epi3 = st->twiddles[fstride*m]; |
|
275 for (i=0;i<N;i++) |
|
276 { |
|
277 Fout = Fout_beg + i*mm; |
|
278 tw1=tw2=st->twiddles; |
|
279 k=m; |
|
280 do{ |
|
281 |
|
282 C_MULC(scratch[1],Fout[m] , *tw1); |
|
283 C_MULC(scratch[2],Fout[m2] , *tw2); |
|
284 |
|
285 C_ADD(scratch[3],scratch[1],scratch[2]); |
|
286 C_SUB(scratch[0],scratch[1],scratch[2]); |
|
287 tw1 += fstride; |
|
288 tw2 += fstride*2; |
|
289 |
|
290 Fout[m].r = Fout->r - HALF_OF(scratch[3].r); |
|
291 Fout[m].i = Fout->i - HALF_OF(scratch[3].i); |
|
292 |
|
293 C_MULBYSCALAR( scratch[0] , -epi3.i ); |
|
294 |
|
295 C_ADDTO(*Fout,scratch[3]); |
|
296 |
|
297 Fout[m2].r = Fout[m].r + scratch[0].i; |
|
298 Fout[m2].i = Fout[m].i - scratch[0].r; |
|
299 |
|
300 Fout[m].r -= scratch[0].i; |
|
301 Fout[m].i += scratch[0].r; |
|
302 |
|
303 ++Fout; |
|
304 }while(--k); |
|
305 } |
|
306 } |
|
307 |
|
308 static void kf_bfly5( |
|
309 kiss_fft_cpx * Fout, |
|
310 const size_t fstride, |
|
311 const kiss_fft_state *st, |
|
312 int m, |
|
313 int N, |
|
314 int mm |
|
315 ) |
|
316 { |
|
317 kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; |
|
318 int i, u; |
|
319 kiss_fft_cpx scratch[13]; |
|
320 const kiss_twiddle_cpx * twiddles = st->twiddles; |
|
321 const kiss_twiddle_cpx *tw; |
|
322 kiss_twiddle_cpx ya,yb; |
|
323 kiss_fft_cpx * Fout_beg = Fout; |
|
324 |
|
325 ya = twiddles[fstride*m]; |
|
326 yb = twiddles[fstride*2*m]; |
|
327 tw=st->twiddles; |
|
328 |
|
329 for (i=0;i<N;i++) |
|
330 { |
|
331 Fout = Fout_beg + i*mm; |
|
332 Fout0=Fout; |
|
333 Fout1=Fout0+m; |
|
334 Fout2=Fout0+2*m; |
|
335 Fout3=Fout0+3*m; |
|
336 Fout4=Fout0+4*m; |
|
337 |
|
338 for ( u=0; u<m; ++u ) { |
|
339 C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5); |
|
340 scratch[0] = *Fout0; |
|
341 |
|
342 C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); |
|
343 C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); |
|
344 C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); |
|
345 C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); |
|
346 |
|
347 C_ADD( scratch[7],scratch[1],scratch[4]); |
|
348 C_SUB( scratch[10],scratch[1],scratch[4]); |
|
349 C_ADD( scratch[8],scratch[2],scratch[3]); |
|
350 C_SUB( scratch[9],scratch[2],scratch[3]); |
|
351 |
|
352 Fout0->r += scratch[7].r + scratch[8].r; |
|
353 Fout0->i += scratch[7].i + scratch[8].i; |
|
354 |
|
355 scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); |
|
356 scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); |
|
357 |
|
358 scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); |
|
359 scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); |
|
360 |
|
361 C_SUB(*Fout1,scratch[5],scratch[6]); |
|
362 C_ADD(*Fout4,scratch[5],scratch[6]); |
|
363 |
|
364 scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); |
|
365 scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); |
|
366 scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); |
|
367 scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); |
|
368 |
|
369 C_ADD(*Fout2,scratch[11],scratch[12]); |
|
370 C_SUB(*Fout3,scratch[11],scratch[12]); |
|
371 |
|
372 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; |
|
373 } |
|
374 } |
|
375 } |
|
376 |
|
377 static void ki_bfly5( |
|
378 kiss_fft_cpx * Fout, |
|
379 const size_t fstride, |
|
380 const kiss_fft_state *st, |
|
381 int m, |
|
382 int N, |
|
383 int mm |
|
384 ) |
|
385 { |
|
386 kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; |
|
387 int i, u; |
|
388 kiss_fft_cpx scratch[13]; |
|
389 const kiss_twiddle_cpx * twiddles = st->twiddles; |
|
390 const kiss_twiddle_cpx *tw; |
|
391 kiss_twiddle_cpx ya,yb; |
|
392 kiss_fft_cpx * Fout_beg = Fout; |
|
393 |
|
394 ya = twiddles[fstride*m]; |
|
395 yb = twiddles[fstride*2*m]; |
|
396 tw=st->twiddles; |
|
397 |
|
398 for (i=0;i<N;i++) |
|
399 { |
|
400 Fout = Fout_beg + i*mm; |
|
401 Fout0=Fout; |
|
402 Fout1=Fout0+m; |
|
403 Fout2=Fout0+2*m; |
|
404 Fout3=Fout0+3*m; |
|
405 Fout4=Fout0+4*m; |
|
406 |
|
407 for ( u=0; u<m; ++u ) { |
|
408 scratch[0] = *Fout0; |
|
409 |
|
410 C_MULC(scratch[1] ,*Fout1, tw[u*fstride]); |
|
411 C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]); |
|
412 C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]); |
|
413 C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]); |
|
414 |
|
415 C_ADD( scratch[7],scratch[1],scratch[4]); |
|
416 C_SUB( scratch[10],scratch[1],scratch[4]); |
|
417 C_ADD( scratch[8],scratch[2],scratch[3]); |
|
418 C_SUB( scratch[9],scratch[2],scratch[3]); |
|
419 |
|
420 Fout0->r += scratch[7].r + scratch[8].r; |
|
421 Fout0->i += scratch[7].i + scratch[8].i; |
|
422 |
|
423 scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); |
|
424 scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); |
|
425 |
|
426 scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i); |
|
427 scratch[6].i = S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i); |
|
428 |
|
429 C_SUB(*Fout1,scratch[5],scratch[6]); |
|
430 C_ADD(*Fout4,scratch[5],scratch[6]); |
|
431 |
|
432 scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); |
|
433 scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); |
|
434 scratch[12].r = S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i); |
|
435 scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i); |
|
436 |
|
437 C_ADD(*Fout2,scratch[11],scratch[12]); |
|
438 C_SUB(*Fout3,scratch[11],scratch[12]); |
|
439 |
|
440 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; |
|
441 } |
|
442 } |
|
443 } |
|
444 |
|
445 #endif |
|
446 |
|
447 |
|
448 #ifdef CUSTOM_MODES |
|
449 |
|
450 static |
|
451 void compute_bitrev_table( |
|
452 int Fout, |
|
453 opus_int16 *f, |
|
454 const size_t fstride, |
|
455 int in_stride, |
|
456 opus_int16 * factors, |
|
457 const kiss_fft_state *st |
|
458 ) |
|
459 { |
|
460 const int p=*factors++; /* the radix */ |
|
461 const int m=*factors++; /* stage's fft length/p */ |
|
462 |
|
463 /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/ |
|
464 if (m==1) |
|
465 { |
|
466 int j; |
|
467 for (j=0;j<p;j++) |
|
468 { |
|
469 *f = Fout+j; |
|
470 f += fstride*in_stride; |
|
471 } |
|
472 } else { |
|
473 int j; |
|
474 for (j=0;j<p;j++) |
|
475 { |
|
476 compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st); |
|
477 f += fstride*in_stride; |
|
478 Fout += m; |
|
479 } |
|
480 } |
|
481 } |
|
482 |
|
483 /* facbuf is populated by p1,m1,p2,m2, ... |
|
484 where |
|
485 p[i] * m[i] = m[i-1] |
|
486 m0 = n */ |
|
487 static |
|
488 int kf_factor(int n,opus_int16 * facbuf) |
|
489 { |
|
490 int p=4; |
|
491 |
|
492 /*factor out powers of 4, powers of 2, then any remaining primes */ |
|
493 do { |
|
494 while (n % p) { |
|
495 switch (p) { |
|
496 case 4: p = 2; break; |
|
497 case 2: p = 3; break; |
|
498 default: p += 2; break; |
|
499 } |
|
500 if (p>32000 || (opus_int32)p*(opus_int32)p > n) |
|
501 p = n; /* no more factors, skip to end */ |
|
502 } |
|
503 n /= p; |
|
504 #ifdef RADIX_TWO_ONLY |
|
505 if (p!=2 && p != 4) |
|
506 #else |
|
507 if (p>5) |
|
508 #endif |
|
509 { |
|
510 return 0; |
|
511 } |
|
512 *facbuf++ = p; |
|
513 *facbuf++ = n; |
|
514 } while (n > 1); |
|
515 return 1; |
|
516 } |
|
517 |
|
518 static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft) |
|
519 { |
|
520 int i; |
|
521 #ifdef FIXED_POINT |
|
522 for (i=0;i<nfft;++i) { |
|
523 opus_val32 phase = -i; |
|
524 kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft)); |
|
525 } |
|
526 #else |
|
527 for (i=0;i<nfft;++i) { |
|
528 const double pi=3.14159265358979323846264338327; |
|
529 double phase = ( -2*pi /nfft ) * i; |
|
530 kf_cexp(twiddles+i, phase ); |
|
531 } |
|
532 #endif |
|
533 } |
|
534 |
|
535 /* |
|
536 * |
|
537 * Allocates all necessary storage space for the fft and ifft. |
|
538 * The return value is a contiguous block of memory. As such, |
|
539 * It can be freed with free(). |
|
540 * */ |
|
541 kiss_fft_state *opus_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem, const kiss_fft_state *base) |
|
542 { |
|
543 kiss_fft_state *st=NULL; |
|
544 size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/ |
|
545 |
|
546 if ( lenmem==NULL ) { |
|
547 st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded ); |
|
548 }else{ |
|
549 if (mem != NULL && *lenmem >= memneeded) |
|
550 st = (kiss_fft_state*)mem; |
|
551 *lenmem = memneeded; |
|
552 } |
|
553 if (st) { |
|
554 opus_int16 *bitrev; |
|
555 kiss_twiddle_cpx *twiddles; |
|
556 |
|
557 st->nfft=nfft; |
|
558 #ifndef FIXED_POINT |
|
559 st->scale = 1.f/nfft; |
|
560 #endif |
|
561 if (base != NULL) |
|
562 { |
|
563 st->twiddles = base->twiddles; |
|
564 st->shift = 0; |
|
565 while (nfft<<st->shift != base->nfft && st->shift < 32) |
|
566 st->shift++; |
|
567 if (st->shift>=32) |
|
568 goto fail; |
|
569 } else { |
|
570 st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft); |
|
571 compute_twiddles(twiddles, nfft); |
|
572 st->shift = -1; |
|
573 } |
|
574 if (!kf_factor(nfft,st->factors)) |
|
575 { |
|
576 goto fail; |
|
577 } |
|
578 |
|
579 /* bitrev */ |
|
580 st->bitrev = bitrev = (opus_int16*)KISS_FFT_MALLOC(sizeof(opus_int16)*nfft); |
|
581 if (st->bitrev==NULL) |
|
582 goto fail; |
|
583 compute_bitrev_table(0, bitrev, 1,1, st->factors,st); |
|
584 } |
|
585 return st; |
|
586 fail: |
|
587 opus_fft_free(st); |
|
588 return NULL; |
|
589 } |
|
590 |
|
591 kiss_fft_state *opus_fft_alloc(int nfft,void * mem,size_t * lenmem ) |
|
592 { |
|
593 return opus_fft_alloc_twiddles(nfft, mem, lenmem, NULL); |
|
594 } |
|
595 |
|
596 void opus_fft_free(const kiss_fft_state *cfg) |
|
597 { |
|
598 if (cfg) |
|
599 { |
|
600 opus_free((opus_int16*)cfg->bitrev); |
|
601 if (cfg->shift < 0) |
|
602 opus_free((kiss_twiddle_cpx*)cfg->twiddles); |
|
603 opus_free((kiss_fft_state*)cfg); |
|
604 } |
|
605 } |
|
606 |
|
607 #endif /* CUSTOM_MODES */ |
|
608 |
|
609 void opus_fft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) |
|
610 { |
|
611 int m2, m; |
|
612 int p; |
|
613 int L; |
|
614 int fstride[MAXFACTORS]; |
|
615 int i; |
|
616 int shift; |
|
617 |
|
618 /* st->shift can be -1 */ |
|
619 shift = st->shift>0 ? st->shift : 0; |
|
620 |
|
621 celt_assert2 (fin != fout, "In-place FFT not supported"); |
|
622 /* Bit-reverse the input */ |
|
623 for (i=0;i<st->nfft;i++) |
|
624 { |
|
625 fout[st->bitrev[i]] = fin[i]; |
|
626 #ifndef FIXED_POINT |
|
627 fout[st->bitrev[i]].r *= st->scale; |
|
628 fout[st->bitrev[i]].i *= st->scale; |
|
629 #endif |
|
630 } |
|
631 |
|
632 fstride[0] = 1; |
|
633 L=0; |
|
634 do { |
|
635 p = st->factors[2*L]; |
|
636 m = st->factors[2*L+1]; |
|
637 fstride[L+1] = fstride[L]*p; |
|
638 L++; |
|
639 } while(m!=1); |
|
640 m = st->factors[2*L-1]; |
|
641 for (i=L-1;i>=0;i--) |
|
642 { |
|
643 if (i!=0) |
|
644 m2 = st->factors[2*i-1]; |
|
645 else |
|
646 m2 = 1; |
|
647 switch (st->factors[2*i]) |
|
648 { |
|
649 case 2: |
|
650 kf_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
651 break; |
|
652 case 4: |
|
653 kf_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
654 break; |
|
655 #ifndef RADIX_TWO_ONLY |
|
656 case 3: |
|
657 kf_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
658 break; |
|
659 case 5: |
|
660 kf_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
661 break; |
|
662 #endif |
|
663 } |
|
664 m = m2; |
|
665 } |
|
666 } |
|
667 |
|
668 void opus_ifft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) |
|
669 { |
|
670 int m2, m; |
|
671 int p; |
|
672 int L; |
|
673 int fstride[MAXFACTORS]; |
|
674 int i; |
|
675 int shift; |
|
676 |
|
677 /* st->shift can be -1 */ |
|
678 shift = st->shift>0 ? st->shift : 0; |
|
679 celt_assert2 (fin != fout, "In-place FFT not supported"); |
|
680 /* Bit-reverse the input */ |
|
681 for (i=0;i<st->nfft;i++) |
|
682 fout[st->bitrev[i]] = fin[i]; |
|
683 |
|
684 fstride[0] = 1; |
|
685 L=0; |
|
686 do { |
|
687 p = st->factors[2*L]; |
|
688 m = st->factors[2*L+1]; |
|
689 fstride[L+1] = fstride[L]*p; |
|
690 L++; |
|
691 } while(m!=1); |
|
692 m = st->factors[2*L-1]; |
|
693 for (i=L-1;i>=0;i--) |
|
694 { |
|
695 if (i!=0) |
|
696 m2 = st->factors[2*i-1]; |
|
697 else |
|
698 m2 = 1; |
|
699 switch (st->factors[2*i]) |
|
700 { |
|
701 case 2: |
|
702 ki_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
703 break; |
|
704 case 4: |
|
705 ki_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
706 break; |
|
707 #ifndef RADIX_TWO_ONLY |
|
708 case 3: |
|
709 ki_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
710 break; |
|
711 case 5: |
|
712 ki_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2); |
|
713 break; |
|
714 #endif |
|
715 } |
|
716 m = m2; |
|
717 } |
|
718 } |
|
719 |