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