Removed the _t from all the celt*_t types to avoid clashing with POSIX
[opus.git] / libcelt / kiss_fft.c
1 /*
2 Copyright (c) 2003-2004, Mark Borgerding
3 Lots of modifications by Jean-Marc Valin
4 Copyright (c) 2005-2007, Xiph.Org Foundation
5 Copyright (c) 2008,      Xiph.Org Foundation, 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 #include "stack_alloc.h"
29
30 /* The guts header contains all the multiplication and addition macros that are defined for
31    complex numbers.  It also delares the kf_ internal functions.
32 */
33
34 static void kf_bfly2(
35                      kiss_fft_cpx * Fout,
36                      const size_t fstride,
37                      const kiss_fft_cfg st,
38                      int m,
39                      int N,
40                      int mm
41                     )
42 {
43    kiss_fft_cpx * Fout2;
44    kiss_twiddle_cpx * tw1;
45    int i,j;
46    kiss_fft_cpx * Fout_beg = Fout;
47    for (i=0;i<N;i++)
48    {
49       Fout = Fout_beg + i*mm;
50       Fout2 = Fout + m;
51       tw1 = st->twiddles;
52       for(j=0;j<m;j++)
53       {
54          kiss_fft_cpx t;
55          Fout->r = SHR(Fout->r, 1);Fout->i = SHR(Fout->i, 1);
56          Fout2->r = SHR(Fout2->r, 1);Fout2->i = SHR(Fout2->i, 1);
57          C_MUL (t,  *Fout2 , *tw1);
58          tw1 += fstride;
59          C_SUB( *Fout2 ,  *Fout , t );
60          C_ADDTO( *Fout ,  t );
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_twiddle_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_twiddle_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 = PSHR(Fout->r, 2);
125          Fout->i = PSHR(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 = PSHR(Fout[m2].r, 2);
131          Fout[m2].i = PSHR(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_twiddle_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 #ifndef RADIX_TWO_ONLY
193
194 static void kf_bfly3(
195                      kiss_fft_cpx * Fout,
196                      const size_t fstride,
197                      const kiss_fft_cfg st,
198                      size_t m
199                     )
200 {
201    size_t k=m;
202    const size_t m2 = 2*m;
203    kiss_twiddle_cpx *tw1,*tw2;
204    kiss_fft_cpx scratch[5];
205    kiss_twiddle_cpx epi3;
206    epi3 = st->twiddles[fstride*m];
207
208    tw1=tw2=st->twiddles;
209    do{
210       C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
211
212       C_MUL(scratch[1],Fout[m] , *tw1);
213       C_MUL(scratch[2],Fout[m2] , *tw2);
214
215       C_ADD(scratch[3],scratch[1],scratch[2]);
216       C_SUB(scratch[0],scratch[1],scratch[2]);
217       tw1 += fstride;
218       tw2 += fstride*2;
219
220       Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
221       Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
222
223       C_MULBYSCALAR( scratch[0] , epi3.i );
224
225       C_ADDTO(*Fout,scratch[3]);
226
227       Fout[m2].r = Fout[m].r + scratch[0].i;
228       Fout[m2].i = Fout[m].i - scratch[0].r;
229
230       Fout[m].r -= scratch[0].i;
231       Fout[m].i += scratch[0].r;
232
233       ++Fout;
234    }while(--k);
235 }
236
237 static void ki_bfly3(
238                      kiss_fft_cpx * Fout,
239                      const size_t fstride,
240                      const kiss_fft_cfg st,
241                      size_t m
242                     )
243 {
244    size_t k=m;
245    const size_t m2 = 2*m;
246    kiss_twiddle_cpx *tw1,*tw2;
247    kiss_fft_cpx scratch[5];
248    kiss_twiddle_cpx epi3;
249    epi3 = st->twiddles[fstride*m];
250
251    tw1=tw2=st->twiddles;
252    do{
253
254       C_MULC(scratch[1],Fout[m] , *tw1);
255       C_MULC(scratch[2],Fout[m2] , *tw2);
256
257       C_ADD(scratch[3],scratch[1],scratch[2]);
258       C_SUB(scratch[0],scratch[1],scratch[2]);
259       tw1 += fstride;
260       tw2 += fstride*2;
261
262       Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
263       Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
264
265       C_MULBYSCALAR( scratch[0] , -epi3.i );
266
267       C_ADDTO(*Fout,scratch[3]);
268
269       Fout[m2].r = Fout[m].r + scratch[0].i;
270       Fout[m2].i = Fout[m].i - scratch[0].r;
271
272       Fout[m].r -= scratch[0].i;
273       Fout[m].i += scratch[0].r;
274
275       ++Fout;
276    }while(--k);
277 }
278
279
280 static void kf_bfly5(
281                      kiss_fft_cpx * Fout,
282                      const size_t fstride,
283                      const kiss_fft_cfg st,
284                      int m
285                     )
286 {
287    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
288    int u;
289    kiss_fft_cpx scratch[13];
290    kiss_twiddle_cpx * twiddles = st->twiddles;
291    kiss_twiddle_cpx *tw;
292    kiss_twiddle_cpx ya,yb;
293    ya = twiddles[fstride*m];
294    yb = twiddles[fstride*2*m];
295
296    Fout0=Fout;
297    Fout1=Fout0+m;
298    Fout2=Fout0+2*m;
299    Fout3=Fout0+3*m;
300    Fout4=Fout0+4*m;
301
302    tw=st->twiddles;
303    for ( u=0; u<m; ++u ) {
304       C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
305       scratch[0] = *Fout0;
306
307       C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
308       C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
309       C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
310       C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
311
312       C_ADD( scratch[7],scratch[1],scratch[4]);
313       C_SUB( scratch[10],scratch[1],scratch[4]);
314       C_ADD( scratch[8],scratch[2],scratch[3]);
315       C_SUB( scratch[9],scratch[2],scratch[3]);
316
317       Fout0->r += scratch[7].r + scratch[8].r;
318       Fout0->i += scratch[7].i + scratch[8].i;
319
320       scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
321       scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
322
323       scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
324       scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
325
326       C_SUB(*Fout1,scratch[5],scratch[6]);
327       C_ADD(*Fout4,scratch[5],scratch[6]);
328
329       scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
330       scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
331       scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
332       scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
333
334       C_ADD(*Fout2,scratch[11],scratch[12]);
335       C_SUB(*Fout3,scratch[11],scratch[12]);
336
337       ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
338    }
339 }
340
341 static void ki_bfly5(
342                      kiss_fft_cpx * Fout,
343                      const size_t fstride,
344                      const kiss_fft_cfg st,
345                      int m
346                     )
347 {
348    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
349    int u;
350    kiss_fft_cpx scratch[13];
351    kiss_twiddle_cpx * twiddles = st->twiddles;
352    kiss_twiddle_cpx *tw;
353    kiss_twiddle_cpx ya,yb;
354    ya = twiddles[fstride*m];
355    yb = twiddles[fstride*2*m];
356
357    Fout0=Fout;
358    Fout1=Fout0+m;
359    Fout2=Fout0+2*m;
360    Fout3=Fout0+3*m;
361    Fout4=Fout0+4*m;
362
363    tw=st->twiddles;
364    for ( u=0; u<m; ++u ) {
365       scratch[0] = *Fout0;
366
367       C_MULC(scratch[1] ,*Fout1, tw[u*fstride]);
368       C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]);
369       C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]);
370       C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]);
371
372       C_ADD( scratch[7],scratch[1],scratch[4]);
373       C_SUB( scratch[10],scratch[1],scratch[4]);
374       C_ADD( scratch[8],scratch[2],scratch[3]);
375       C_SUB( scratch[9],scratch[2],scratch[3]);
376
377       Fout0->r += scratch[7].r + scratch[8].r;
378       Fout0->i += scratch[7].i + scratch[8].i;
379
380       scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
381       scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
382
383       scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i);
384       scratch[6].i =  S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i);
385
386       C_SUB(*Fout1,scratch[5],scratch[6]);
387       C_ADD(*Fout4,scratch[5],scratch[6]);
388
389       scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
390       scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
391       scratch[12].r =  S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i);
392       scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i);
393
394       C_ADD(*Fout2,scratch[11],scratch[12]);
395       C_SUB(*Fout3,scratch[11],scratch[12]);
396
397       ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
398    }
399 }
400
401 /* perform the butterfly for one stage of a mixed radix FFT */
402 static void kf_bfly_generic(
403                             kiss_fft_cpx * Fout,
404                             const size_t fstride,
405                             const kiss_fft_cfg st,
406                             int m,
407                             int p
408                            )
409 {
410    int u,k,q1,q;
411    kiss_twiddle_cpx * twiddles = st->twiddles;
412    kiss_fft_cpx t;
413    VARDECL(kiss_fft_cpx, scratchbuf);
414    int Norig = st->nfft;
415    ALLOC(scratchbuf, p, kiss_fft_cpx);
416
417    for ( u=0; u<m; ++u ) {
418       k=u;
419       for ( q1=0 ; q1<p ; ++q1 ) {
420          scratchbuf[q1] = Fout[ k  ];
421          C_FIXDIV(scratchbuf[q1],p);
422          k += m;
423       }
424
425       k=u;
426       for ( q1=0 ; q1<p ; ++q1 ) {
427          int twidx=0;
428          Fout[ k ] = scratchbuf[0];
429          for (q=1;q<p;++q ) {
430             twidx += fstride * k;
431             if (twidx>=Norig) twidx-=Norig;
432             C_MUL(t,scratchbuf[q] , twiddles[twidx] );
433             C_ADDTO( Fout[ k ] ,t);
434          }
435          k += m;
436       }
437    }
438 }
439
440 static void ki_bfly_generic(
441                             kiss_fft_cpx * Fout,
442                             const size_t fstride,
443                             const kiss_fft_cfg st,
444                             int m,
445                             int p
446                            )
447 {
448    int u,k,q1,q;
449    kiss_twiddle_cpx * twiddles = st->twiddles;
450    kiss_fft_cpx t;
451    VARDECL(kiss_fft_cpx, scratchbuf);
452    int Norig = st->nfft;
453    ALLOC(scratchbuf, p, kiss_fft_cpx);
454
455    for ( u=0; u<m; ++u ) {
456       k=u;
457       for ( q1=0 ; q1<p ; ++q1 ) {
458          scratchbuf[q1] = Fout[ k  ];
459          k += m;
460       }
461
462       k=u;
463       for ( q1=0 ; q1<p ; ++q1 ) {
464          int twidx=0;
465          Fout[ k ] = scratchbuf[0];
466          for (q=1;q<p;++q ) {
467             twidx += fstride * k;
468             if (twidx>=Norig) twidx-=Norig;
469             C_MULC(t,scratchbuf[q] , twiddles[twidx] );
470             C_ADDTO( Fout[ k ] ,t);
471          }
472          k += m;
473       }
474    }
475 }
476 #endif
477
478 static
479 void compute_bitrev_table(
480          int Fout,
481          int *f,
482          const size_t fstride,
483          int in_stride,
484          int * factors,
485          const kiss_fft_cfg st
486             )
487 {
488    const int p=*factors++; /* the radix  */
489    const int m=*factors++; /* stage's fft length/p */
490    
491     /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
492    if (m==1)
493    {
494       int j;
495       for (j=0;j<p;j++)
496       {
497          *f = Fout+j;
498          f += fstride*in_stride;
499       }
500    } else {
501       int j;
502       for (j=0;j<p;j++)
503       {
504          compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
505          f += fstride*in_stride;
506          Fout += m;
507       }
508    }
509 }
510
511
512 void kf_work(
513         kiss_fft_cpx * Fout,
514         const kiss_fft_cpx * f,
515         const size_t fstride,
516         int in_stride,
517         int * factors,
518         const kiss_fft_cfg st,
519         int N,
520         int s2,
521         int m2
522         )
523 {
524 #ifndef RADIX_TWO_ONLY
525     int i;
526     kiss_fft_cpx * Fout_beg=Fout;
527 #endif
528     const int p=*factors++; /* the radix  */
529     const int m=*factors++; /* stage's fft length/p */
530     /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
531     if (m!=1) 
532         kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
533
534     switch (p) {
535         case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break;
536         case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
537 #ifndef RADIX_TWO_ONLY
538         case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly3(Fout,fstride,st,m);} break; 
539         case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly5(Fout,fstride,st,m);} break; 
540         default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly_generic(Fout,fstride,st,m,p);} break;
541 #else
542        default: celt_fatal("kiss_fft: only powers of two enabled");
543 #endif
544     }    
545 }
546
547
548 void ki_work(
549              kiss_fft_cpx * Fout,
550              const kiss_fft_cpx * f,
551              const size_t fstride,
552              int in_stride,
553              int * factors,
554              const kiss_fft_cfg st,
555              int N,
556              int s2,
557              int m2
558             )
559 {
560 #ifndef RADIX_TWO_ONLY
561    int i;
562    kiss_fft_cpx * Fout_beg=Fout;
563 #endif
564    const int p=*factors++; /* the radix  */
565    const int m=*factors++; /* stage's fft length/p */
566    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
567    if (m!=1) 
568       ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
569
570    switch (p) {
571       case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
572       case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
573 #ifndef RADIX_TWO_ONLY
574       case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly3(Fout,fstride,st,m);} break; 
575       case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly5(Fout,fstride,st,m);} break; 
576       default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly_generic(Fout,fstride,st,m,p);} break;
577 #else
578       default: celt_fatal("kiss_fft: only powers of two enabled");
579 #endif
580    }    
581 }
582
583 /*  facbuf is populated by p1,m1,p2,m2, ...
584     where 
585     p[i] * m[i] = m[i-1]
586     m0 = n                  */
587 static 
588 void kf_factor(int n,int * facbuf)
589 {
590     int p=4;
591
592     /*factor out powers of 4, powers of 2, then any remaining primes */
593     do {
594         while (n % p) {
595             switch (p) {
596                 case 4: p = 2; break;
597                 case 2: p = 3; break;
598                 default: p += 2; break;
599             }
600             if (p>32000 || (celt_int32)p*(celt_int32)p > n)
601                 p = n;          /* no more factors, skip to end */
602         }
603         n /= p;
604         *facbuf++ = p;
605         *facbuf++ = n;
606     } while (n > 1);
607 }
608 /*
609  *
610  * User-callable function to allocate all necessary storage space for the fft.
611  *
612  * The return value is a contiguous block of memory, allocated with malloc.  As such,
613  * It can be freed with free(), rather than a kiss_fft-specific function.
614  * */
615 kiss_fft_cfg kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
616 {
617     kiss_fft_cfg st=NULL;
618     size_t memneeded = sizeof(struct kiss_fft_state)
619           + sizeof(kiss_twiddle_cpx)*(nfft-1) + sizeof(int)*nfft; /* twiddle factors*/
620
621     if ( lenmem==NULL ) {
622         st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
623     }else{
624         if (mem != NULL && *lenmem >= memneeded)
625             st = (kiss_fft_cfg)mem;
626         *lenmem = memneeded;
627     }
628     if (st) {
629         int i;
630         st->nfft=nfft;
631 #ifndef FIXED_POINT
632         st->scale = 1./nfft;
633 #endif
634 #if defined(FIXED_POINT) && (!defined(DOUBLE_PRECISION) || defined(MIXED_PRECISION))
635         for (i=0;i<nfft;++i) {
636             celt_word32 phase = -i;
637             kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft));
638         }
639 #else
640         for (i=0;i<nfft;++i) {
641            const double pi=3.14159265358979323846264338327;
642            double phase = ( -2*pi /nfft ) * i;
643            kf_cexp(st->twiddles+i, phase );
644         }
645 #endif
646         kf_factor(nfft,st->factors);
647         
648         /* bitrev */
649         st->bitrev = (int*)((char*)st + memneeded - sizeof(int)*nfft);
650         compute_bitrev_table(0, st->bitrev, 1,1, st->factors,st);
651     }
652     return st;
653 }
654
655
656
657     
658 void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
659 {
660     if (fin == fout) 
661     {
662        celt_fatal("In-place FFT not supported");
663     } else {
664        /* Bit-reverse the input */
665        int i;
666        for (i=0;i<st->nfft;i++)
667        {
668           fout[st->bitrev[i]] = fin[i];
669 #ifndef FIXED_POINT
670           fout[st->bitrev[i]].r *= st->scale;
671           fout[st->bitrev[i]].i *= st->scale;
672 #endif
673        }
674        kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
675     }
676 }
677
678 void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
679 {
680     kiss_fft_stride(cfg,fin,fout,1);
681 }
682
683 void kiss_ifft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
684 {
685    if (fin == fout) 
686    {
687       celt_fatal("In-place FFT not supported");
688    } else {
689       /* Bit-reverse the input */
690       int i;
691       for (i=0;i<st->nfft;i++)
692          fout[st->bitrev[i]] = fin[i];
693       ki_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
694    }
695 }
696
697 void kiss_ifft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
698 {
699    kiss_ifft_stride(cfg,fin,fout,1);
700 }
701