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