Applying the forward FFT gain up-front for fixed-point too
[opus.git] / celt / kiss_fft.c
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          C_MUL (t,  *Fout2 , *tw1);
70          tw1 += fstride;
71          C_SUB( *Fout2 ,  *Fout , t );
72          C_ADDTO( *Fout ,  t );
73          ++Fout2;
74          ++Fout;
75       }
76    }
77 }
78
79 static void ki_bfly2(
80                      kiss_fft_cpx * Fout,
81                      const size_t fstride,
82                      const kiss_fft_state *st,
83                      int m,
84                      int N,
85                      int mm
86                     )
87 {
88    kiss_fft_cpx * Fout2;
89    const kiss_twiddle_cpx * tw1;
90    kiss_fft_cpx t;
91    int i,j;
92    kiss_fft_cpx * Fout_beg = Fout;
93    for (i=0;i<N;i++)
94    {
95       Fout = Fout_beg + i*mm;
96       Fout2 = Fout + m;
97       tw1 = st->twiddles;
98       for(j=0;j<m;j++)
99       {
100          C_MULC (t,  *Fout2 , *tw1);
101          tw1 += fstride;
102          C_SUB( *Fout2 ,  *Fout , t );
103          C_ADDTO( *Fout ,  t );
104          ++Fout2;
105          ++Fout;
106       }
107    }
108 }
109
110 static void kf_bfly4(
111                      kiss_fft_cpx * Fout,
112                      const size_t fstride,
113                      const kiss_fft_state *st,
114                      int m,
115                      int N,
116                      int mm
117                     )
118 {
119    const kiss_twiddle_cpx *tw1,*tw2,*tw3;
120    kiss_fft_cpx scratch[6];
121    const size_t m2=2*m;
122    const size_t m3=3*m;
123    int i, j;
124
125    kiss_fft_cpx * Fout_beg = Fout;
126    for (i=0;i<N;i++)
127    {
128       Fout = Fout_beg + i*mm;
129       tw3 = tw2 = tw1 = st->twiddles;
130       for (j=0;j<m;j++)
131       {
132          C_MUL(scratch[0],Fout[m] , *tw1 );
133          C_MUL(scratch[1],Fout[m2] , *tw2 );
134          C_MUL(scratch[2],Fout[m3] , *tw3 );
135
136          C_SUB( scratch[5] , *Fout, scratch[1] );
137          C_ADDTO(*Fout, scratch[1]);
138          C_ADD( scratch[3] , scratch[0] , scratch[2] );
139          C_SUB( scratch[4] , scratch[0] , scratch[2] );
140          C_SUB( Fout[m2], *Fout, scratch[3] );
141          tw1 += fstride;
142          tw2 += fstride*2;
143          tw3 += fstride*3;
144          C_ADDTO( *Fout , scratch[3] );
145
146          Fout[m].r = scratch[5].r + scratch[4].i;
147          Fout[m].i = scratch[5].i - scratch[4].r;
148          Fout[m3].r = scratch[5].r - scratch[4].i;
149          Fout[m3].i = scratch[5].i + scratch[4].r;
150          ++Fout;
151       }
152    }
153 }
154
155 static void ki_bfly4(
156                      kiss_fft_cpx * Fout,
157                      const size_t fstride,
158                      const kiss_fft_state *st,
159                      int m,
160                      int N,
161                      int mm
162                     )
163 {
164    const kiss_twiddle_cpx *tw1,*tw2,*tw3;
165    kiss_fft_cpx scratch[6];
166    const size_t m2=2*m;
167    const size_t m3=3*m;
168    int i, j;
169
170    kiss_fft_cpx * Fout_beg = Fout;
171    for (i=0;i<N;i++)
172    {
173       Fout = Fout_beg + i*mm;
174       tw3 = tw2 = tw1 = st->twiddles;
175       for (j=0;j<m;j++)
176       {
177          C_MULC(scratch[0],Fout[m] , *tw1 );
178          C_MULC(scratch[1],Fout[m2] , *tw2 );
179          C_MULC(scratch[2],Fout[m3] , *tw3 );
180
181          C_SUB( scratch[5] , *Fout, scratch[1] );
182          C_ADDTO(*Fout, scratch[1]);
183          C_ADD( scratch[3] , scratch[0] , scratch[2] );
184          C_SUB( scratch[4] , scratch[0] , scratch[2] );
185          C_SUB( Fout[m2], *Fout, scratch[3] );
186          tw1 += fstride;
187          tw2 += fstride*2;
188          tw3 += fstride*3;
189          C_ADDTO( *Fout , scratch[3] );
190
191          Fout[m].r = scratch[5].r - scratch[4].i;
192          Fout[m].i = scratch[5].i + scratch[4].r;
193          Fout[m3].r = scratch[5].r + scratch[4].i;
194          Fout[m3].i = scratch[5].i - scratch[4].r;
195          ++Fout;
196       }
197    }
198 }
199
200 #ifndef RADIX_TWO_ONLY
201
202 static void kf_bfly3(
203                      kiss_fft_cpx * Fout,
204                      const size_t fstride,
205                      const kiss_fft_state *st,
206                      int m,
207                      int N,
208                      int mm
209                     )
210 {
211    int i;
212    size_t k;
213    const size_t m2 = 2*m;
214    const kiss_twiddle_cpx *tw1,*tw2;
215    kiss_fft_cpx scratch[5];
216    kiss_twiddle_cpx epi3;
217
218    kiss_fft_cpx * Fout_beg = Fout;
219    epi3 = st->twiddles[fstride*m];
220    for (i=0;i<N;i++)
221    {
222       Fout = Fout_beg + i*mm;
223       tw1=tw2=st->twiddles;
224       k=m;
225       do {
226
227          C_MUL(scratch[1],Fout[m] , *tw1);
228          C_MUL(scratch[2],Fout[m2] , *tw2);
229
230          C_ADD(scratch[3],scratch[1],scratch[2]);
231          C_SUB(scratch[0],scratch[1],scratch[2]);
232          tw1 += fstride;
233          tw2 += fstride*2;
234
235          Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
236          Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
237
238          C_MULBYSCALAR( scratch[0] , epi3.i );
239
240          C_ADDTO(*Fout,scratch[3]);
241
242          Fout[m2].r = Fout[m].r + scratch[0].i;
243          Fout[m2].i = Fout[m].i - scratch[0].r;
244
245          Fout[m].r -= scratch[0].i;
246          Fout[m].i += scratch[0].r;
247
248          ++Fout;
249       } while(--k);
250    }
251 }
252
253 static void ki_bfly3(
254                      kiss_fft_cpx * Fout,
255                      const size_t fstride,
256                      const kiss_fft_state *st,
257                      int m,
258                      int N,
259                      int mm
260                     )
261 {
262    int i, k;
263    const size_t m2 = 2*m;
264    const kiss_twiddle_cpx *tw1,*tw2;
265    kiss_fft_cpx scratch[5];
266    kiss_twiddle_cpx epi3;
267
268    kiss_fft_cpx * Fout_beg = Fout;
269    epi3 = st->twiddles[fstride*m];
270    for (i=0;i<N;i++)
271    {
272       Fout = Fout_beg + i*mm;
273       tw1=tw2=st->twiddles;
274       k=m;
275       do{
276
277          C_MULC(scratch[1],Fout[m] , *tw1);
278          C_MULC(scratch[2],Fout[m2] , *tw2);
279
280          C_ADD(scratch[3],scratch[1],scratch[2]);
281          C_SUB(scratch[0],scratch[1],scratch[2]);
282          tw1 += fstride;
283          tw2 += fstride*2;
284
285          Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
286          Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
287
288          C_MULBYSCALAR( scratch[0] , -epi3.i );
289
290          C_ADDTO(*Fout,scratch[3]);
291
292          Fout[m2].r = Fout[m].r + scratch[0].i;
293          Fout[m2].i = Fout[m].i - scratch[0].r;
294
295          Fout[m].r -= scratch[0].i;
296          Fout[m].i += scratch[0].r;
297
298          ++Fout;
299       }while(--k);
300    }
301 }
302
303 static void kf_bfly5(
304                      kiss_fft_cpx * Fout,
305                      const size_t fstride,
306                      const kiss_fft_state *st,
307                      int m,
308                      int N,
309                      int mm
310                     )
311 {
312    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
313    int i, u;
314    kiss_fft_cpx scratch[13];
315    const kiss_twiddle_cpx * twiddles = st->twiddles;
316    const kiss_twiddle_cpx *tw;
317    kiss_twiddle_cpx ya,yb;
318    kiss_fft_cpx * Fout_beg = Fout;
319
320    ya = twiddles[fstride*m];
321    yb = twiddles[fstride*2*m];
322    tw=st->twiddles;
323
324    for (i=0;i<N;i++)
325    {
326       Fout = Fout_beg + i*mm;
327       Fout0=Fout;
328       Fout1=Fout0+m;
329       Fout2=Fout0+2*m;
330       Fout3=Fout0+3*m;
331       Fout4=Fout0+4*m;
332
333       for ( u=0; u<m; ++u ) {
334          scratch[0] = *Fout0;
335
336          C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
337          C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
338          C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
339          C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
340
341          C_ADD( scratch[7],scratch[1],scratch[4]);
342          C_SUB( scratch[10],scratch[1],scratch[4]);
343          C_ADD( scratch[8],scratch[2],scratch[3]);
344          C_SUB( scratch[9],scratch[2],scratch[3]);
345
346          Fout0->r += scratch[7].r + scratch[8].r;
347          Fout0->i += scratch[7].i + scratch[8].i;
348
349          scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
350          scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
351
352          scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
353          scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
354
355          C_SUB(*Fout1,scratch[5],scratch[6]);
356          C_ADD(*Fout4,scratch[5],scratch[6]);
357
358          scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
359          scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
360          scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
361          scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
362
363          C_ADD(*Fout2,scratch[11],scratch[12]);
364          C_SUB(*Fout3,scratch[11],scratch[12]);
365
366          ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
367       }
368    }
369 }
370
371 static void ki_bfly5(
372                      kiss_fft_cpx * Fout,
373                      const size_t fstride,
374                      const kiss_fft_state *st,
375                      int m,
376                      int N,
377                      int mm
378                     )
379 {
380    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
381    int i, u;
382    kiss_fft_cpx scratch[13];
383    const kiss_twiddle_cpx * twiddles = st->twiddles;
384    const kiss_twiddle_cpx *tw;
385    kiss_twiddle_cpx ya,yb;
386    kiss_fft_cpx * Fout_beg = Fout;
387
388    ya = twiddles[fstride*m];
389    yb = twiddles[fstride*2*m];
390    tw=st->twiddles;
391
392    for (i=0;i<N;i++)
393    {
394       Fout = Fout_beg + i*mm;
395       Fout0=Fout;
396       Fout1=Fout0+m;
397       Fout2=Fout0+2*m;
398       Fout3=Fout0+3*m;
399       Fout4=Fout0+4*m;
400
401       for ( u=0; u<m; ++u ) {
402          scratch[0] = *Fout0;
403
404          C_MULC(scratch[1] ,*Fout1, tw[u*fstride]);
405          C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]);
406          C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]);
407          C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]);
408
409          C_ADD( scratch[7],scratch[1],scratch[4]);
410          C_SUB( scratch[10],scratch[1],scratch[4]);
411          C_ADD( scratch[8],scratch[2],scratch[3]);
412          C_SUB( scratch[9],scratch[2],scratch[3]);
413
414          Fout0->r += scratch[7].r + scratch[8].r;
415          Fout0->i += scratch[7].i + scratch[8].i;
416
417          scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
418          scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
419
420          scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i);
421          scratch[6].i =  S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i);
422
423          C_SUB(*Fout1,scratch[5],scratch[6]);
424          C_ADD(*Fout4,scratch[5],scratch[6]);
425
426          scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
427          scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
428          scratch[12].r =  S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i);
429          scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i);
430
431          C_ADD(*Fout2,scratch[11],scratch[12]);
432          C_SUB(*Fout3,scratch[11],scratch[12]);
433
434          ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
435       }
436    }
437 }
438
439 #endif
440
441
442 #ifdef CUSTOM_MODES
443
444 static
445 void compute_bitrev_table(
446          int Fout,
447          opus_int16 *f,
448          const size_t fstride,
449          int in_stride,
450          opus_int16 * factors,
451          const kiss_fft_state *st
452             )
453 {
454    const int p=*factors++; /* the radix  */
455    const int m=*factors++; /* stage's fft length/p */
456
457     /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
458    if (m==1)
459    {
460       int j;
461       for (j=0;j<p;j++)
462       {
463          *f = Fout+j;
464          f += fstride*in_stride;
465       }
466    } else {
467       int j;
468       for (j=0;j<p;j++)
469       {
470          compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
471          f += fstride*in_stride;
472          Fout += m;
473       }
474    }
475 }
476
477 /*  facbuf is populated by p1,m1,p2,m2, ...
478     where
479     p[i] * m[i] = m[i-1]
480     m0 = n                  */
481 static
482 int kf_factor(int n,opus_int16 * facbuf)
483 {
484     int p=4;
485
486     /*factor out powers of 4, powers of 2, then any remaining primes */
487     do {
488         while (n % p) {
489             switch (p) {
490                 case 4: p = 2; break;
491                 case 2: p = 3; break;
492                 default: p += 2; break;
493             }
494             if (p>32000 || (opus_int32)p*(opus_int32)p > n)
495                 p = n;          /* no more factors, skip to end */
496         }
497         n /= p;
498 #ifdef RADIX_TWO_ONLY
499         if (p!=2 && p != 4)
500 #else
501         if (p>5)
502 #endif
503         {
504            return 0;
505         }
506         *facbuf++ = p;
507         *facbuf++ = n;
508     } while (n > 1);
509     return 1;
510 }
511
512 static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft)
513 {
514    int i;
515 #ifdef FIXED_POINT
516    for (i=0;i<nfft;++i) {
517       opus_val32 phase = -i;
518       kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft));
519    }
520 #else
521    for (i=0;i<nfft;++i) {
522       const double pi=3.14159265358979323846264338327;
523       double phase = ( -2*pi /nfft ) * i;
524       kf_cexp(twiddles+i, phase );
525    }
526 #endif
527 }
528
529 /*
530  *
531  * Allocates all necessary storage space for the fft and ifft.
532  * The return value is a contiguous block of memory.  As such,
533  * It can be freed with free().
534  * */
535 kiss_fft_state *opus_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem,  const kiss_fft_state *base)
536 {
537     kiss_fft_state *st=NULL;
538     size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/
539
540     if ( lenmem==NULL ) {
541         st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded );
542     }else{
543         if (mem != NULL && *lenmem >= memneeded)
544             st = (kiss_fft_state*)mem;
545         *lenmem = memneeded;
546     }
547     if (st) {
548         opus_int16 *bitrev;
549         kiss_twiddle_cpx *twiddles;
550
551         st->nfft=nfft;
552 #ifndef FIXED_POINT
553         st->scale = 1.f/nfft;
554 #endif
555         if (base != NULL)
556         {
557            st->twiddles = base->twiddles;
558            st->shift = 0;
559            while (nfft<<st->shift != base->nfft && st->shift < 32)
560               st->shift++;
561            if (st->shift>=32)
562               goto fail;
563         } else {
564            st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft);
565            compute_twiddles(twiddles, nfft);
566            st->shift = -1;
567         }
568         if (!kf_factor(nfft,st->factors))
569         {
570            goto fail;
571         }
572
573         /* bitrev */
574         st->bitrev = bitrev = (opus_int16*)KISS_FFT_MALLOC(sizeof(opus_int16)*nfft);
575         if (st->bitrev==NULL)
576             goto fail;
577         compute_bitrev_table(0, bitrev, 1,1, st->factors,st);
578     }
579     return st;
580 fail:
581     opus_fft_free(st);
582     return NULL;
583 }
584
585 kiss_fft_state *opus_fft_alloc(int nfft,void * mem,size_t * lenmem )
586 {
587    return opus_fft_alloc_twiddles(nfft, mem, lenmem, NULL);
588 }
589
590 void opus_fft_free(const kiss_fft_state *cfg)
591 {
592    if (cfg)
593    {
594       opus_free((opus_int16*)cfg->bitrev);
595       if (cfg->shift < 0)
596          opus_free((kiss_twiddle_cpx*)cfg->twiddles);
597       opus_free((kiss_fft_state*)cfg);
598    }
599 }
600
601 #endif /* CUSTOM_MODES */
602
603 void opus_fft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
604 {
605     int m2, m;
606     int p;
607     int L;
608     int fstride[MAXFACTORS];
609     int i;
610     int shift;
611 #ifdef FIXED_POINT
612     /* FIXME: This should eventually just go in the state. */
613     opus_val16 scale;
614     int scale_shift;
615     scale_shift = celt_ilog2(st->nfft);
616     if (st->nfft == 1<<scale_shift)
617        scale = Q15ONE;
618     else
619        scale = (1073741824+st->nfft/2)/st->nfft>>(15-scale_shift);
620 #endif
621
622     /* st->shift can be -1 */
623     shift = st->shift>0 ? st->shift : 0;
624
625     celt_assert2 (fin != fout, "In-place FFT not supported");
626     /* Bit-reverse the input */
627     for (i=0;i<st->nfft;i++)
628     {
629        kiss_fft_cpx x = fin[i];
630 #ifdef FIXED_POINT
631        fout[st->bitrev[i]].r = SHR32(MULT16_32_Q15(scale, x.r), scale_shift);
632        fout[st->bitrev[i]].i = SHR32(MULT16_32_Q15(scale, x.i), scale_shift);
633 #else
634        fout[st->bitrev[i]].r = st->scale*x.r;
635        fout[st->bitrev[i]].i = st->scale*x.i;
636 #endif
637     }
638
639     fstride[0] = 1;
640     L=0;
641     do {
642        p = st->factors[2*L];
643        m = st->factors[2*L+1];
644        fstride[L+1] = fstride[L]*p;
645        L++;
646     } while(m!=1);
647     m = st->factors[2*L-1];
648     for (i=L-1;i>=0;i--)
649     {
650        if (i!=0)
651           m2 = st->factors[2*i-1];
652        else
653           m2 = 1;
654        switch (st->factors[2*i])
655        {
656        case 2:
657           kf_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2);
658           break;
659        case 4:
660           kf_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2);
661           break;
662  #ifndef RADIX_TWO_ONLY
663        case 3:
664           kf_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2);
665           break;
666        case 5:
667           kf_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2);
668           break;
669  #endif
670        }
671        m = m2;
672     }
673 }
674
675 void opus_ifft_impl(const kiss_fft_state *st,kiss_fft_cpx *fout)
676 {
677    int m2, m;
678    int p;
679    int L;
680    int fstride[MAXFACTORS];
681    int i;
682    int shift;
683
684    /* st->shift can be -1 */
685    shift = st->shift>0 ? st->shift : 0;
686    fstride[0] = 1;
687    L=0;
688    do {
689       p = st->factors[2*L];
690       m = st->factors[2*L+1];
691       fstride[L+1] = fstride[L]*p;
692       L++;
693    } while(m!=1);
694    m = st->factors[2*L-1];
695    for (i=L-1;i>=0;i--)
696    {
697       if (i!=0)
698          m2 = st->factors[2*i-1];
699       else
700          m2 = 1;
701       switch (st->factors[2*i])
702       {
703       case 2:
704          ki_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2);
705          break;
706       case 4:
707          ki_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2);
708          break;
709 #ifndef RADIX_TWO_ONLY
710       case 3:
711          ki_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2);
712          break;
713       case 5:
714          ki_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2);
715          break;
716 #endif
717       }
718       m = m2;
719    }
720 }
721
722 #ifdef TEST_UNIT_DFT_C
723 void opus_ifft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
724 {
725    int i;
726    celt_assert2 (fin != fout, "In-place FFT not supported");
727    /* Bit-reverse the input */
728    for (i=0;i<st->nfft;i++)
729       fout[st->bitrev[i]] = fin[i];
730    opus_ifft_impl(st, fout);
731 }
732 #endif