Updated Doxygen comments, removed an incorrectly placed LGPL header (we own
[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_fft_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              /* Almost the same as the code path below, except that we divide the input by two
52          (while keeping the best accuracy possible) */
53          celt_word32_t tr, ti;
54          tr = SHR32(SUB32(MULT16_16(Fout2->r , tw1->r),MULT16_16(Fout2->i , tw1->i)), 1);
55          ti = SHR32(ADD32(MULT16_16(Fout2->i , tw1->r),MULT16_16(Fout2->r , tw1->i)), 1);
56          tw1 += fstride;
57          Fout2->r = PSHR32(SUB32(SHL32(EXTEND32(Fout->r), 14), tr), 15);
58          Fout2->i = PSHR32(SUB32(SHL32(EXTEND32(Fout->i), 14), ti), 15);
59          Fout->r = PSHR32(ADD32(SHL32(EXTEND32(Fout->r), 14), tr), 15);
60          Fout->i = PSHR32(ADD32(SHL32(EXTEND32(Fout->i), 14), ti), 15);
61          ++Fout2;
62          ++Fout;
63       }
64    }
65 }
66
67 static void ki_bfly2(
68                      kiss_fft_cpx * Fout,
69                      const size_t fstride,
70                      const kiss_fft_cfg st,
71                      int m,
72                      int N,
73                      int mm
74                     )
75 {
76    kiss_fft_cpx * Fout2;
77    kiss_fft_cpx * tw1;
78    kiss_fft_cpx t;
79    int i,j;
80    kiss_fft_cpx * Fout_beg = Fout;
81    for (i=0;i<N;i++)
82    {
83       Fout = Fout_beg + i*mm;
84       Fout2 = Fout + m;
85       tw1 = st->twiddles;
86       for(j=0;j<m;j++)
87       {
88          C_MULC (t,  *Fout2 , *tw1);
89          tw1 += fstride;
90          C_SUB( *Fout2 ,  *Fout , t );
91          C_ADDTO( *Fout ,  t );
92          ++Fout2;
93          ++Fout;
94       }
95    }
96 }
97
98 static void kf_bfly4(
99                      kiss_fft_cpx * Fout,
100                      const size_t fstride,
101                      const kiss_fft_cfg st,
102                      int m,
103                      int N,
104                      int mm
105                     )
106 {
107    kiss_fft_cpx *tw1,*tw2,*tw3;
108    kiss_fft_cpx scratch[6];
109    const size_t m2=2*m;
110    const size_t m3=3*m;
111    int i, j;
112
113    kiss_fft_cpx * Fout_beg = Fout;
114    for (i=0;i<N;i++)
115    {
116       Fout = Fout_beg + i*mm;
117       tw3 = tw2 = tw1 = st->twiddles;
118       for (j=0;j<m;j++)
119       {
120          C_MUL4(scratch[0],Fout[m] , *tw1 );
121          C_MUL4(scratch[1],Fout[m2] , *tw2 );
122          C_MUL4(scratch[2],Fout[m3] , *tw3 );
123              
124          Fout->r = PSHR16(Fout->r, 2);
125          Fout->i = PSHR16(Fout->i, 2);
126          C_SUB( scratch[5] , *Fout, scratch[1] );
127          C_ADDTO(*Fout, scratch[1]);
128          C_ADD( scratch[3] , scratch[0] , scratch[2] );
129          C_SUB( scratch[4] , scratch[0] , scratch[2] );
130          Fout[m2].r = PSHR16(Fout[m2].r, 2);
131          Fout[m2].i = PSHR16(Fout[m2].i, 2);
132          C_SUB( Fout[m2], *Fout, scratch[3] );
133          tw1 += fstride;
134          tw2 += fstride*2;
135          tw3 += fstride*3;
136          C_ADDTO( *Fout , scratch[3] );
137              
138          Fout[m].r = scratch[5].r + scratch[4].i;
139          Fout[m].i = scratch[5].i - scratch[4].r;
140          Fout[m3].r = scratch[5].r - scratch[4].i;
141          Fout[m3].i = scratch[5].i + scratch[4].r;
142          ++Fout;
143       }
144    }
145 }
146
147 static void ki_bfly4(
148                      kiss_fft_cpx * Fout,
149                      const size_t fstride,
150                      const kiss_fft_cfg st,
151                      int m,
152                      int N,
153                      int mm
154                     )
155 {
156    kiss_fft_cpx *tw1,*tw2,*tw3;
157    kiss_fft_cpx scratch[6];
158    const size_t m2=2*m;
159    const size_t m3=3*m;
160    int i, j;
161
162    kiss_fft_cpx * Fout_beg = Fout;
163    for (i=0;i<N;i++)
164    {
165       Fout = Fout_beg + i*mm;
166       tw3 = tw2 = tw1 = st->twiddles;
167       for (j=0;j<m;j++)
168       {
169          C_MULC(scratch[0],Fout[m] , *tw1 );
170          C_MULC(scratch[1],Fout[m2] , *tw2 );
171          C_MULC(scratch[2],Fout[m3] , *tw3 );
172              
173          C_SUB( scratch[5] , *Fout, scratch[1] );
174          C_ADDTO(*Fout, scratch[1]);
175          C_ADD( scratch[3] , scratch[0] , scratch[2] );
176          C_SUB( scratch[4] , scratch[0] , scratch[2] );
177          C_SUB( Fout[m2], *Fout, scratch[3] );
178          tw1 += fstride;
179          tw2 += fstride*2;
180          tw3 += fstride*3;
181          C_ADDTO( *Fout , scratch[3] );
182              
183          Fout[m].r = scratch[5].r - scratch[4].i;
184          Fout[m].i = scratch[5].i + scratch[4].r;
185          Fout[m3].r = scratch[5].r + scratch[4].i;
186          Fout[m3].i = scratch[5].i - scratch[4].r;
187          ++Fout;
188       }
189    }
190 }
191
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_fft_cpx *tw1,*tw2;
203    kiss_fft_cpx scratch[5];
204    kiss_fft_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_fft_cpx *tw1,*tw2;
246    kiss_fft_cpx scratch[5];
247    kiss_fft_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_fft_cpx * twiddles = st->twiddles;
290    kiss_fft_cpx *tw;
291    kiss_fft_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_fft_cpx * twiddles = st->twiddles;
351    kiss_fft_cpx *tw;
352    kiss_fft_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_fft_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_fft_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
482 static
483 void compute_bitrev_table(
484          int * Fout,
485          int f,
486          const size_t fstride,
487          int in_stride,
488          int * factors,
489          const kiss_fft_cfg st
490             )
491 {
492    const int p=*factors++; /* the radix  */
493    const int m=*factors++; /* stage's fft length/p */
494    
495     /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
496    if (m==1)
497    {
498       int j;
499       for (j=0;j<p;j++)
500       {
501          Fout[j] = f;
502          f += fstride*in_stride;
503       }
504    } else {
505       int j;
506       for (j=0;j<p;j++)
507       {
508          compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
509          f += fstride*in_stride;
510          Fout += m;
511       }
512    }
513 }
514
515 static
516 void kf_work(
517         kiss_fft_cpx * Fout,
518         const kiss_fft_cpx * f,
519         const size_t fstride,
520         int in_stride,
521         int * factors,
522         const kiss_fft_cfg st,
523         int N,
524         int s2,
525         int m2
526         )
527 {
528     int i;
529     kiss_fft_cpx * Fout_beg=Fout;
530     const int p=*factors++; /* the radix  */
531     const int m=*factors++; /* stage's fft length/p */
532     /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
533     if (m!=1) 
534         kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
535
536     switch (p) {
537         case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break;
538         case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly3(Fout,fstride,st,m);} break; 
539         case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
540         case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly5(Fout,fstride,st,m);} break; 
541         default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly_generic(Fout,fstride,st,m,p);} break;
542     }    
543 }
544
545 static
546       void ki_work(
547                    kiss_fft_cpx * Fout,
548                    const kiss_fft_cpx * f,
549                    const size_t fstride,
550                    int in_stride,
551                    int * factors,
552                    const kiss_fft_cfg st,
553                    int N,
554                    int s2,
555                    int m2
556                   )
557 {
558    int i;
559    kiss_fft_cpx * Fout_beg=Fout;
560    const int p=*factors++; /* the radix  */
561    const int m=*factors++; /* stage's fft length/p */
562    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
563    if (m!=1) 
564       ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
565
566    switch (p) {
567       case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
568       case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly3(Fout,fstride,st,m);} break; 
569       case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
570       case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly5(Fout,fstride,st,m);} break; 
571       default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly_generic(Fout,fstride,st,m,p);} break;
572    }    
573 }
574
575 /*  facbuf is populated by p1,m1,p2,m2, ...
576     where 
577     p[i] * m[i] = m[i-1]
578     m0 = n                  */
579 static 
580 void kf_factor(int n,int * facbuf)
581 {
582     int p=4;
583
584     /*factor out powers of 4, powers of 2, then any remaining primes */
585     do {
586         while (n % p) {
587             switch (p) {
588                 case 4: p = 2; break;
589                 case 2: p = 3; break;
590                 default: p += 2; break;
591             }
592             if (p>32000 || (celt_int32_t)p*(celt_int32_t)p > n)
593                 p = n;          /* no more factors, skip to end */
594         }
595         n /= p;
596         *facbuf++ = p;
597         *facbuf++ = n;
598     } while (n > 1);
599 }
600 /*
601  *
602  * User-callable function to allocate all necessary storage space for the fft.
603  *
604  * The return value is a contiguous block of memory, allocated with malloc.  As such,
605  * It can be freed with free(), rather than a kiss_fft-specific function.
606  * */
607 kiss_fft_cfg kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
608 {
609     kiss_fft_cfg st=NULL;
610     size_t memneeded = sizeof(struct kiss_fft_state)
611           + sizeof(kiss_fft_cpx)*(nfft-1) + sizeof(int)*nfft; /* twiddle factors*/
612
613     if ( lenmem==NULL ) {
614         st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
615     }else{
616         if (mem != NULL && *lenmem >= memneeded)
617             st = (kiss_fft_cfg)mem;
618         *lenmem = memneeded;
619     }
620     if (st) {
621         int i;
622         st->nfft=nfft;
623 #ifdef FIXED_POINT
624         for (i=0;i<nfft;++i) {
625             celt_word32_t phase = -i;
626             kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft));
627         }
628 #else
629         for (i=0;i<nfft;++i) {
630            const double pi=3.14159265358979323846264338327;
631            double phase = ( -2*pi /nfft ) * i;
632            kf_cexp(st->twiddles+i, phase );
633         }
634 #endif
635         kf_factor(nfft,st->factors);
636         
637         /* bitrev */
638         st->bitrev = (int*)((char*)st + memneeded - sizeof(int)*nfft);
639         compute_bitrev_table(st->bitrev, 0, 1,1, st->factors,st);
640     }
641     return st;
642 }
643
644
645
646     
647 void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
648 {
649     if (fin == fout) 
650     {
651        celt_fatal("In-place FFT not supported");
652     } else {
653        /* Bit-reverse the input */
654        int i;
655        for (i=0;i<st->nfft;i++)
656           fout[i] = fin[st->bitrev[i]];
657        kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
658     }
659 }
660
661 void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
662 {
663     kiss_fft_stride(cfg,fin,fout,1);
664 }
665
666 void kiss_ifft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
667 {
668    if (fin == fout) 
669    {
670       celt_fatal("In-place FFT not supported");
671    } else {
672       /* Bit-reverse the input */
673       int i;
674       for (i=0;i<st->nfft;i++)
675          fout[i] = fin[st->bitrev[i]];
676       ki_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
677    }
678 }
679
680 void kiss_ifft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
681 {
682    kiss_ifft_stride(cfg,fin,fout,1);
683 }
684