fixed leaked ritrev table
[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    kiss_fft_cpx t;
43    int i,j;
44    kiss_fft_cpx * Fout_beg = Fout;
45    for (i=0;i<N;i++)
46    {
47       Fout = Fout_beg + i*mm;
48       Fout2 = Fout + m;
49       tw1 = st->twiddles;
50       for(j=0;j<m;j++)
51       {
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          ++Fout2;
63          ++Fout;
64       }
65    }
66 }
67
68 static void ki_bfly2(
69                      kiss_fft_cpx * Fout,
70                      const size_t fstride,
71                      const kiss_fft_cfg st,
72                      int m,
73                      int N,
74                      int mm
75                     )
76 {
77    kiss_fft_cpx * Fout2;
78    kiss_fft_cpx * tw1;
79    kiss_fft_cpx t;
80    int i,j;
81    kiss_fft_cpx * Fout_beg = Fout;
82    for (i=0;i<N;i++)
83    {
84       Fout = Fout_beg + i*mm;
85       Fout2 = Fout + m;
86       tw1 = st->twiddles;
87       for(j=0;j<m;j++)
88       {
89          C_MULC (t,  *Fout2 , *tw1);
90          tw1 += fstride;
91          C_SUB( *Fout2 ,  *Fout , t );
92          C_ADDTO( *Fout ,  t );
93          ++Fout2;
94          ++Fout;
95       }
96    }
97 }
98
99 static void kf_bfly4(
100                      kiss_fft_cpx * Fout,
101                      const size_t fstride,
102                      const kiss_fft_cfg st,
103                      int m,
104                      int N,
105                      int mm
106                     )
107 {
108    kiss_fft_cpx *tw1,*tw2,*tw3;
109    kiss_fft_cpx scratch[6];
110    const size_t m2=2*m;
111    const size_t m3=3*m;
112    int i, j;
113
114    kiss_fft_cpx * Fout_beg = Fout;
115    for (i=0;i<N;i++)
116    {
117       Fout = Fout_beg + i*mm;
118       tw3 = tw2 = tw1 = st->twiddles;
119       for (j=0;j<m;j++)
120       {
121          C_MUL4(scratch[0],Fout[m] , *tw1 );
122          C_MUL4(scratch[1],Fout[m2] , *tw2 );
123          C_MUL4(scratch[2],Fout[m3] , *tw3 );
124              
125          Fout->r = PSHR16(Fout->r, 2);
126          Fout->i = PSHR16(Fout->i, 2);
127          C_SUB( scratch[5] , *Fout, scratch[1] );
128          C_ADDTO(*Fout, scratch[1]);
129          C_ADD( scratch[3] , scratch[0] , scratch[2] );
130          C_SUB( scratch[4] , scratch[0] , scratch[2] );
131          Fout[m2].r = PSHR16(Fout[m2].r, 2);
132          Fout[m2].i = PSHR16(Fout[m2].i, 2);
133          C_SUB( Fout[m2], *Fout, scratch[3] );
134          tw1 += fstride;
135          tw2 += fstride*2;
136          tw3 += fstride*3;
137          C_ADDTO( *Fout , scratch[3] );
138              
139          Fout[m].r = scratch[5].r + scratch[4].i;
140          Fout[m].i = scratch[5].i - scratch[4].r;
141          Fout[m3].r = scratch[5].r - scratch[4].i;
142          Fout[m3].i = scratch[5].i + scratch[4].r;
143          ++Fout;
144       }
145    }
146 }
147
148 static void ki_bfly4(
149                      kiss_fft_cpx * Fout,
150                      const size_t fstride,
151                      const kiss_fft_cfg st,
152                      int m,
153                      int N,
154                      int mm
155                     )
156 {
157    kiss_fft_cpx *tw1,*tw2,*tw3;
158    kiss_fft_cpx scratch[6];
159    const size_t m2=2*m;
160    const size_t m3=3*m;
161    int i, j;
162
163    kiss_fft_cpx * Fout_beg = Fout;
164    for (i=0;i<N;i++)
165    {
166       Fout = Fout_beg + i*mm;
167       tw3 = tw2 = tw1 = st->twiddles;
168       for (j=0;j<m;j++)
169       {
170          C_MULC(scratch[0],Fout[m] , *tw1 );
171          C_MULC(scratch[1],Fout[m2] , *tw2 );
172          C_MULC(scratch[2],Fout[m3] , *tw3 );
173              
174          C_SUB( scratch[5] , *Fout, scratch[1] );
175          C_ADDTO(*Fout, scratch[1]);
176          C_ADD( scratch[3] , scratch[0] , scratch[2] );
177          C_SUB( scratch[4] , scratch[0] , scratch[2] );
178          C_SUB( Fout[m2], *Fout, scratch[3] );
179          tw1 += fstride;
180          tw2 += fstride*2;
181          tw3 += fstride*3;
182          C_ADDTO( *Fout , scratch[3] );
183              
184          Fout[m].r = scratch[5].r - scratch[4].i;
185          Fout[m].i = scratch[5].i + scratch[4].r;
186          Fout[m3].r = scratch[5].r + scratch[4].i;
187          Fout[m3].i = scratch[5].i - scratch[4].r;
188          ++Fout;
189       }
190    }
191 }
192
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_fft_cpx *tw1,*tw2;
204    kiss_fft_cpx scratch[5];
205    kiss_fft_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_fft_cpx *tw1,*tw2;
247    kiss_fft_cpx scratch[5];
248    kiss_fft_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_fft_cpx * twiddles = st->twiddles;
291    kiss_fft_cpx *tw;
292    kiss_fft_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_fft_cpx * twiddles = st->twiddles;
352    kiss_fft_cpx *tw;
353    kiss_fft_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_fft_cpx * twiddles = st->twiddles;
412    kiss_fft_cpx t;
413    kiss_fft_cpx scratchbuf[17];
414    int Norig = st->nfft;
415
416    /*CHECKBUF(scratchbuf,nscratchbuf,p);*/
417    if (p>17)
418       celt_fatal("KissFFT: max radix supported is 17");
419     
420    for ( u=0; u<m; ++u ) {
421       k=u;
422       for ( q1=0 ; q1<p ; ++q1 ) {
423          scratchbuf[q1] = Fout[ k  ];
424          C_FIXDIV(scratchbuf[q1],p);
425          k += m;
426       }
427
428       k=u;
429       for ( q1=0 ; q1<p ; ++q1 ) {
430          int twidx=0;
431          Fout[ k ] = scratchbuf[0];
432          for (q=1;q<p;++q ) {
433             twidx += fstride * k;
434             if (twidx>=Norig) twidx-=Norig;
435             C_MUL(t,scratchbuf[q] , twiddles[twidx] );
436             C_ADDTO( Fout[ k ] ,t);
437          }
438          k += m;
439       }
440    }
441 }
442
443 static void ki_bfly_generic(
444                             kiss_fft_cpx * Fout,
445                             const size_t fstride,
446                             const kiss_fft_cfg st,
447                             int m,
448                             int p
449                            )
450 {
451    int u,k,q1,q;
452    kiss_fft_cpx * twiddles = st->twiddles;
453    kiss_fft_cpx t;
454    kiss_fft_cpx scratchbuf[17];
455    int Norig = st->nfft;
456
457    /*CHECKBUF(scratchbuf,nscratchbuf,p);*/
458    if (p>17)
459       celt_fatal("KissFFT: max radix supported is 17");
460     
461    for ( u=0; u<m; ++u ) {
462       k=u;
463       for ( q1=0 ; q1<p ; ++q1 ) {
464          scratchbuf[q1] = Fout[ k  ];
465          k += m;
466       }
467
468       k=u;
469       for ( q1=0 ; q1<p ; ++q1 ) {
470          int twidx=0;
471          Fout[ k ] = scratchbuf[0];
472          for (q=1;q<p;++q ) {
473             twidx += fstride * k;
474             if (twidx>=Norig) twidx-=Norig;
475             C_MULC(t,scratchbuf[q] , twiddles[twidx] );
476             C_ADDTO( Fout[ k ] ,t);
477          }
478          k += m;
479       }
480    }
481 }
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          Fout[j] = f;
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 static
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 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly3(Fout,fstride,st,m);} break; 
540         case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
541         case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly5(Fout,fstride,st,m);} break; 
542         default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly_generic(Fout,fstride,st,m,p);} break;
543     }    
544 }
545
546 static
547       void ki_work(
548                    kiss_fft_cpx * Fout,
549                    const kiss_fft_cpx * f,
550                    const size_t fstride,
551                    int in_stride,
552                    int * factors,
553                    const kiss_fft_cfg st,
554                    int N,
555                    int s2,
556                    int m2
557                   )
558 {
559    int i;
560    kiss_fft_cpx * Fout_beg=Fout;
561    const int p=*factors++; /* the radix  */
562    const int m=*factors++; /* stage's fft length/p */
563    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
564    if (m!=1) 
565       ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
566
567    switch (p) {
568       case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
569       case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly3(Fout,fstride,st,m);} break; 
570       case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
571       case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly5(Fout,fstride,st,m);} break; 
572       default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; ki_bfly_generic(Fout,fstride,st,m,p);} break;
573    }    
574 }
575
576 /*  facbuf is populated by p1,m1,p2,m2, ...
577     where 
578     p[i] * m[i] = m[i-1]
579     m0 = n                  */
580 static 
581 void kf_factor(int n,int * facbuf)
582 {
583     int p=4;
584
585     /*factor out powers of 4, powers of 2, then any remaining primes */
586     do {
587         while (n % p) {
588             switch (p) {
589                 case 4: p = 2; break;
590                 case 2: p = 3; break;
591                 default: p += 2; break;
592             }
593             if (p>32000 || (celt_int32_t)p*(celt_int32_t)p > n)
594                 p = n;          /* no more factors, skip to end */
595         }
596         n /= p;
597         *facbuf++ = p;
598         *facbuf++ = n;
599     } while (n > 1);
600 }
601 /*
602  *
603  * User-callable function to allocate all necessary storage space for the fft.
604  *
605  * The return value is a contiguous block of memory, allocated with malloc.  As such,
606  * It can be freed with free(), rather than a kiss_fft-specific function.
607  * */
608 kiss_fft_cfg kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
609 {
610     kiss_fft_cfg st=NULL;
611     size_t memneeded = sizeof(struct kiss_fft_state)
612           + sizeof(kiss_fft_cpx)*(nfft-1) + sizeof(int)*nfft; /* twiddle factors*/
613
614     if ( lenmem==NULL ) {
615         st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
616     }else{
617         if (mem != NULL && *lenmem >= memneeded)
618             st = (kiss_fft_cfg)mem;
619         *lenmem = memneeded;
620     }
621     if (st) {
622         int i;
623         st->nfft=nfft;
624 #ifdef FIXED_POINT
625         for (i=0;i<nfft;++i) {
626             celt_word32_t phase = -i;
627             kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft));
628         }
629 #else
630         for (i=0;i<nfft;++i) {
631            const double pi=3.14159265358979323846264338327;
632            double phase = ( -2*pi /nfft ) * i;
633            kf_cexp(st->twiddles+i, phase );
634         }
635 #endif
636         kf_factor(nfft,st->factors);
637         
638         /* bitrev */
639         st->bitrev = (int*)((char*)st + memneeded - sizeof(int)*nfft);
640         compute_bitrev_table(st->bitrev, 0, 1,1, st->factors,st);
641     }
642     return st;
643 }
644
645
646
647     
648 void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
649 {
650     if (fin == fout) 
651     {
652        celt_fatal("In-place FFT not supported");
653     } else {
654        /* Bit-reverse the input */
655        int i;
656        for (i=0;i<st->nfft;i++)
657           fout[i] = fin[st->bitrev[i]];
658        kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
659     }
660 }
661
662 void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
663 {
664     kiss_fft_stride(cfg,fin,fout,1);
665 }
666
667 void kiss_ifft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
668 {
669    if (fin == fout) 
670    {
671       celt_fatal("In-place FFT not supported");
672    } else {
673       /* Bit-reverse the input */
674       int i;
675       for (i=0;i<st->nfft;i++)
676          fout[i] = fin[st->bitrev[i]];
677       ki_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
678    }
679 }
680
681 void kiss_ifft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
682 {
683    kiss_ifft_stride(cfg,fin,fout,1);
684 }
685