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