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