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