Making a bunch of things static
[opus.git] / libcelt / kiss_fft.c
1 /*
2 Copyright (c) 2003-2004, Mark Borgerding
3 Lots of modifications by Jean-Marc Valin
4 Copyright (c) 2005-2007, Xiph.Org Foundation
5 Copyright (c) 2008,      Xiph.Org Foundation, CSIRO
6
7 All rights reserved.
8
9 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
10
11     * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
12     * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
13     * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.
14
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
16 */
17
18 #ifndef SKIP_CONFIG_H
19 #  ifdef HAVE_CONFIG_H
20 #    include "config.h"
21 #  endif
22 #endif
23
24 #include "_kiss_fft_guts.h"
25 #include "arch.h"
26 #include "os_support.h"
27 #include "mathops.h"
28 #include "stack_alloc.h"
29
30 /* The guts header contains all the multiplication and addition macros that are defined for
31    complex numbers.  It also delares the kf_ internal functions.
32 */
33
34 static void kf_bfly2(
35                      kiss_fft_cpx * Fout,
36                      const size_t fstride,
37                      const kiss_fft_cfg st,
38                      int m,
39                      int N,
40                      int mm
41                     )
42 {
43    kiss_fft_cpx * Fout2;
44    kiss_twiddle_cpx * tw1;
45    int i,j;
46    kiss_fft_cpx * Fout_beg = Fout;
47    for (i=0;i<N;i++)
48    {
49       Fout = Fout_beg + i*mm;
50       Fout2 = Fout + m;
51       tw1 = st->twiddles;
52       for(j=0;j<m;j++)
53       {
54          kiss_fft_cpx t;
55          Fout->r = SHR(Fout->r, 1);Fout->i = SHR(Fout->i, 1);
56          Fout2->r = SHR(Fout2->r, 1);Fout2->i = SHR(Fout2->i, 1);
57          C_MUL (t,  *Fout2 , *tw1);
58          tw1 += fstride;
59          C_SUB( *Fout2 ,  *Fout , t );
60          C_ADDTO( *Fout ,  t );
61          ++Fout2;
62          ++Fout;
63       }
64    }
65 }
66
67 static void ki_bfly2(
68                      kiss_fft_cpx * Fout,
69                      const size_t fstride,
70                      const kiss_fft_cfg st,
71                      int m,
72                      int N,
73                      int mm
74                     )
75 {
76    kiss_fft_cpx * Fout2;
77    kiss_twiddle_cpx * tw1;
78    kiss_fft_cpx t;
79    int i,j;
80    kiss_fft_cpx * Fout_beg = Fout;
81    for (i=0;i<N;i++)
82    {
83       Fout = Fout_beg + i*mm;
84       Fout2 = Fout + m;
85       tw1 = st->twiddles;
86       for(j=0;j<m;j++)
87       {
88          C_MULC (t,  *Fout2 , *tw1);
89          tw1 += fstride;
90          C_SUB( *Fout2 ,  *Fout , t );
91          C_ADDTO( *Fout ,  t );
92          ++Fout2;
93          ++Fout;
94       }
95    }
96 }
97
98 static void kf_bfly4(
99                      kiss_fft_cpx * Fout,
100                      const size_t fstride,
101                      const kiss_fft_cfg st,
102                      int m,
103                      int N,
104                      int mm
105                     )
106 {
107    kiss_twiddle_cpx *tw1,*tw2,*tw3;
108    kiss_fft_cpx scratch[6];
109    const size_t m2=2*m;
110    const size_t m3=3*m;
111    int i, j;
112
113    kiss_fft_cpx * Fout_beg = Fout;
114    for (i=0;i<N;i++)
115    {
116       Fout = Fout_beg + i*mm;
117       tw3 = tw2 = tw1 = st->twiddles;
118       for (j=0;j<m;j++)
119       {
120          C_MUL4(scratch[0],Fout[m] , *tw1 );
121          C_MUL4(scratch[1],Fout[m2] , *tw2 );
122          C_MUL4(scratch[2],Fout[m3] , *tw3 );
123              
124          Fout->r = PSHR(Fout->r, 2);
125          Fout->i = PSHR(Fout->i, 2);
126          C_SUB( scratch[5] , *Fout, scratch[1] );
127          C_ADDTO(*Fout, scratch[1]);
128          C_ADD( scratch[3] , scratch[0] , scratch[2] );
129          C_SUB( scratch[4] , scratch[0] , scratch[2] );
130          Fout[m2].r = PSHR(Fout[m2].r, 2);
131          Fout[m2].i = PSHR(Fout[m2].i, 2);
132          C_SUB( Fout[m2], *Fout, scratch[3] );
133          tw1 += fstride;
134          tw2 += fstride*2;
135          tw3 += fstride*3;
136          C_ADDTO( *Fout , scratch[3] );
137              
138          Fout[m].r = scratch[5].r + scratch[4].i;
139          Fout[m].i = scratch[5].i - scratch[4].r;
140          Fout[m3].r = scratch[5].r - scratch[4].i;
141          Fout[m3].i = scratch[5].i + scratch[4].r;
142          ++Fout;
143       }
144    }
145 }
146
147 static void ki_bfly4(
148                      kiss_fft_cpx * Fout,
149                      const size_t fstride,
150                      const kiss_fft_cfg st,
151                      int m,
152                      int N,
153                      int mm
154                     )
155 {
156    kiss_twiddle_cpx *tw1,*tw2,*tw3;
157    kiss_fft_cpx scratch[6];
158    const size_t m2=2*m;
159    const size_t m3=3*m;
160    int i, j;
161
162    kiss_fft_cpx * Fout_beg = Fout;
163    for (i=0;i<N;i++)
164    {
165       Fout = Fout_beg + i*mm;
166       tw3 = tw2 = tw1 = st->twiddles;
167       for (j=0;j<m;j++)
168       {
169          C_MULC(scratch[0],Fout[m] , *tw1 );
170          C_MULC(scratch[1],Fout[m2] , *tw2 );
171          C_MULC(scratch[2],Fout[m3] , *tw3 );
172              
173          C_SUB( scratch[5] , *Fout, scratch[1] );
174          C_ADDTO(*Fout, scratch[1]);
175          C_ADD( scratch[3] , scratch[0] , scratch[2] );
176          C_SUB( scratch[4] , scratch[0] , scratch[2] );
177          C_SUB( Fout[m2], *Fout, scratch[3] );
178          tw1 += fstride;
179          tw2 += fstride*2;
180          tw3 += fstride*3;
181          C_ADDTO( *Fout , scratch[3] );
182              
183          Fout[m].r = scratch[5].r - scratch[4].i;
184          Fout[m].i = scratch[5].i + scratch[4].r;
185          Fout[m3].r = scratch[5].r + scratch[4].i;
186          Fout[m3].i = scratch[5].i - scratch[4].r;
187          ++Fout;
188       }
189    }
190 }
191
192 #ifndef RADIX_TWO_ONLY
193
194 static void kf_bfly3(
195                      kiss_fft_cpx * Fout,
196                      const size_t fstride,
197                      const kiss_fft_cfg st,
198                      int m,
199                      int N,
200                      int mm
201                     )
202 {
203    int i;
204    size_t k;
205    const size_t m2 = 2*m;
206    kiss_twiddle_cpx *tw1,*tw2;
207    kiss_fft_cpx scratch[5];
208    kiss_twiddle_cpx epi3;
209
210    kiss_fft_cpx * Fout_beg = Fout;
211    epi3 = st->twiddles[fstride*m];
212    for (i=0;i<N;i++)
213    {
214       Fout = Fout_beg + i*mm;
215       tw1=tw2=st->twiddles;
216       k=m;
217       do {
218          C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
219
220          C_MUL(scratch[1],Fout[m] , *tw1);
221          C_MUL(scratch[2],Fout[m2] , *tw2);
222
223          C_ADD(scratch[3],scratch[1],scratch[2]);
224          C_SUB(scratch[0],scratch[1],scratch[2]);
225          tw1 += fstride;
226          tw2 += fstride*2;
227
228          Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
229          Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
230
231          C_MULBYSCALAR( scratch[0] , epi3.i );
232
233          C_ADDTO(*Fout,scratch[3]);
234
235          Fout[m2].r = Fout[m].r + scratch[0].i;
236          Fout[m2].i = Fout[m].i - scratch[0].r;
237
238          Fout[m].r -= scratch[0].i;
239          Fout[m].i += scratch[0].r;
240
241          ++Fout;
242       } while(--k);
243    }
244 }
245
246 static void ki_bfly3(
247                      kiss_fft_cpx * Fout,
248                      const size_t fstride,
249                      const kiss_fft_cfg st,
250                      size_t m,
251                      int N,
252                      int mm
253                     )
254 {
255    size_t i, k;
256    const size_t m2 = 2*m;
257    kiss_twiddle_cpx *tw1,*tw2;
258    kiss_fft_cpx scratch[5];
259    kiss_twiddle_cpx epi3;
260
261    kiss_fft_cpx * Fout_beg = Fout;
262    epi3 = st->twiddles[fstride*m];
263    for (i=0;i<N;i++)
264    {
265       Fout = Fout_beg + i*mm;
266       tw1=tw2=st->twiddles;
267       k=m;
268       do{
269
270          C_MULC(scratch[1],Fout[m] , *tw1);
271          C_MULC(scratch[2],Fout[m2] , *tw2);
272
273          C_ADD(scratch[3],scratch[1],scratch[2]);
274          C_SUB(scratch[0],scratch[1],scratch[2]);
275          tw1 += fstride;
276          tw2 += fstride*2;
277
278          Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
279          Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
280
281          C_MULBYSCALAR( scratch[0] , -epi3.i );
282
283          C_ADDTO(*Fout,scratch[3]);
284
285          Fout[m2].r = Fout[m].r + scratch[0].i;
286          Fout[m2].i = Fout[m].i - scratch[0].r;
287
288          Fout[m].r -= scratch[0].i;
289          Fout[m].i += scratch[0].r;
290
291          ++Fout;
292       }while(--k);
293    }
294 }
295
296
297 static void kf_bfly5(
298                      kiss_fft_cpx * Fout,
299                      const size_t fstride,
300                      const kiss_fft_cfg st,
301                      int m,
302                      int N,
303                      int mm
304                     )
305 {
306    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
307    int i, u;
308    kiss_fft_cpx scratch[13];
309    kiss_twiddle_cpx * twiddles = st->twiddles;
310    kiss_twiddle_cpx *tw;
311    kiss_twiddle_cpx ya,yb;
312    kiss_fft_cpx * Fout_beg = Fout;
313
314    ya = twiddles[fstride*m];
315    yb = twiddles[fstride*2*m];
316    tw=st->twiddles;
317
318    for (i=0;i<N;i++)
319    {
320       Fout = Fout_beg + i*mm;
321       Fout0=Fout;
322       Fout1=Fout0+m;
323       Fout2=Fout0+2*m;
324       Fout3=Fout0+3*m;
325       Fout4=Fout0+4*m;
326
327       for ( u=0; u<m; ++u ) {
328          C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
329          scratch[0] = *Fout0;
330
331          C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
332          C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
333          C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
334          C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
335
336          C_ADD( scratch[7],scratch[1],scratch[4]);
337          C_SUB( scratch[10],scratch[1],scratch[4]);
338          C_ADD( scratch[8],scratch[2],scratch[3]);
339          C_SUB( scratch[9],scratch[2],scratch[3]);
340
341          Fout0->r += scratch[7].r + scratch[8].r;
342          Fout0->i += scratch[7].i + scratch[8].i;
343
344          scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
345          scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
346
347          scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
348          scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
349
350          C_SUB(*Fout1,scratch[5],scratch[6]);
351          C_ADD(*Fout4,scratch[5],scratch[6]);
352
353          scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
354          scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
355          scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
356          scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
357
358          C_ADD(*Fout2,scratch[11],scratch[12]);
359          C_SUB(*Fout3,scratch[11],scratch[12]);
360
361          ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
362       }
363    }
364 }
365
366 static void ki_bfly5(
367                      kiss_fft_cpx * Fout,
368                      const size_t fstride,
369                      const kiss_fft_cfg st,
370                      int m,
371                      int N,
372                      int mm
373                     )
374 {
375    kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
376    int i, u;
377    kiss_fft_cpx scratch[13];
378    kiss_twiddle_cpx * twiddles = st->twiddles;
379    kiss_twiddle_cpx *tw;
380    kiss_twiddle_cpx ya,yb;
381    kiss_fft_cpx * Fout_beg = Fout;
382
383    ya = twiddles[fstride*m];
384    yb = twiddles[fstride*2*m];
385    tw=st->twiddles;
386
387    for (i=0;i<N;i++)
388    {
389       Fout = Fout_beg + i*mm;
390       Fout0=Fout;
391       Fout1=Fout0+m;
392       Fout2=Fout0+2*m;
393       Fout3=Fout0+3*m;
394       Fout4=Fout0+4*m;
395
396       for ( u=0; u<m; ++u ) {
397          scratch[0] = *Fout0;
398
399          C_MULC(scratch[1] ,*Fout1, tw[u*fstride]);
400          C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]);
401          C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]);
402          C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]);
403
404          C_ADD( scratch[7],scratch[1],scratch[4]);
405          C_SUB( scratch[10],scratch[1],scratch[4]);
406          C_ADD( scratch[8],scratch[2],scratch[3]);
407          C_SUB( scratch[9],scratch[2],scratch[3]);
408
409          Fout0->r += scratch[7].r + scratch[8].r;
410          Fout0->i += scratch[7].i + scratch[8].i;
411
412          scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
413          scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
414
415          scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i);
416          scratch[6].i =  S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i);
417
418          C_SUB(*Fout1,scratch[5],scratch[6]);
419          C_ADD(*Fout4,scratch[5],scratch[6]);
420
421          scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
422          scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
423          scratch[12].r =  S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i);
424          scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i);
425
426          C_ADD(*Fout2,scratch[11],scratch[12]);
427          C_SUB(*Fout3,scratch[11],scratch[12]);
428
429          ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
430       }
431    }
432 }
433
434 #endif
435
436 static
437 void compute_bitrev_table(
438          int Fout,
439          celt_int16 *f,
440          const size_t fstride,
441          int in_stride,
442          int * factors,
443          const kiss_fft_cfg st
444             )
445 {
446    const int p=*factors++; /* the radix  */
447    const int m=*factors++; /* stage's fft length/p */
448    
449     /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
450    if (m==1)
451    {
452       int j;
453       for (j=0;j<p;j++)
454       {
455          *f = Fout+j;
456          f += fstride*in_stride;
457       }
458    } else {
459       int j;
460       for (j=0;j<p;j++)
461       {
462          compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
463          f += fstride*in_stride;
464          Fout += m;
465       }
466    }
467 }
468
469
470 static void kf_work(
471         kiss_fft_cpx * Fout,
472         const kiss_fft_cpx * f,
473         size_t fstride,
474         int in_stride,
475         int * factors,
476         const kiss_fft_cfg st,
477         int N,
478         int s2,
479         int m2
480         )
481 {
482     const int p=*factors++; /* the radix  */
483     const int m=*factors++; /* stage's fft length/p */
484     /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
485     if (m!=1) 
486         kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
487
488     /* Compensate for longer twiddles table (when sharing) */
489     if (st->shift>0)
490        fstride <<= st->shift;
491     switch (p) {
492         case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break;
493         case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
494 #ifndef RADIX_TWO_ONLY
495         case 3: kf_bfly3(Fout,fstride,st,m, N, m2); break;
496         case 5: kf_bfly5(Fout,fstride,st,m, N, m2); break;
497 #else
498        default: celt_fatal("kiss_fft: only powers of two enabled");
499 #endif
500     }    
501 }
502
503
504 static void ki_work(
505              kiss_fft_cpx * Fout,
506              const kiss_fft_cpx * f,
507              size_t fstride,
508              int in_stride,
509              int * factors,
510              const kiss_fft_cfg st,
511              int N,
512              int s2,
513              int m2
514             )
515 {
516    const int p=*factors++; /* the radix  */
517    const int m=*factors++; /* stage's fft length/p */
518    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
519    if (m!=1) 
520       ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
521
522    /* Compensate for longer twiddles table (when sharing) */
523    if (st->shift>0)
524       fstride <<= st->shift;
525    switch (p) {
526       case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
527       case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
528 #ifndef RADIX_TWO_ONLY
529       case 3: ki_bfly3(Fout,fstride,st,m, N, m2); break;
530       case 5: ki_bfly5(Fout,fstride,st,m, N, m2); break;
531 #else
532       default: celt_fatal("kiss_fft: only powers of two enabled");
533 #endif
534    }    
535 }
536
537 /*  facbuf is populated by p1,m1,p2,m2, ...
538     where 
539     p[i] * m[i] = m[i-1]
540     m0 = n                  */
541 static 
542 int kf_factor(int n,int * facbuf)
543 {
544     int p=4;
545
546     /*factor out powers of 4, powers of 2, then any remaining primes */
547     do {
548         while (n % p) {
549             switch (p) {
550                 case 4: p = 2; break;
551                 case 2: p = 3; break;
552                 default: p += 2; break;
553             }
554             if (p>32000 || (celt_int32)p*(celt_int32)p > n)
555                 p = n;          /* no more factors, skip to end */
556         }
557         n /= p;
558         if (p>5)
559         {
560            celt_warning("Only powers of 2, 3 and 5 are supported");
561            return 0;
562         }
563         *facbuf++ = p;
564         *facbuf++ = n;
565     } while (n > 1);
566     return 1;
567 }
568
569 static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft)
570 {
571    int i;
572 #if defined(FIXED_POINT) && (!defined(DOUBLE_PRECISION) || defined(MIXED_PRECISION))
573    for (i=0;i<nfft;++i) {
574       celt_word32 phase = -i;
575       kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft));
576    }
577 #else
578    for (i=0;i<nfft;++i) {
579       const double pi=3.14159265358979323846264338327;
580       double phase = ( -2*pi /nfft ) * i;
581       kf_cexp(twiddles+i, phase );
582    }
583 #endif
584 }
585
586
587 /*
588  *
589  * User-callable function to allocate all necessary storage space for the fft.
590  *
591  * The return value is a contiguous block of memory, allocated with malloc.  As such,
592  * It can be freed with free(), rather than a kiss_fft-specific function.
593  * */
594 kiss_fft_cfg kiss_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem,  kiss_fft_cfg base)
595 {
596     kiss_fft_cfg st=NULL;
597     size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/
598
599     if ( lenmem==NULL ) {
600         st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
601     }else{
602         if (mem != NULL && *lenmem >= memneeded)
603             st = (kiss_fft_cfg)mem;
604         *lenmem = memneeded;
605     }
606     if (st) {
607         st->nfft=nfft;
608 #ifndef FIXED_POINT
609         st->scale = 1./nfft;
610 #endif
611         if (base != NULL)
612         {
613            st->twiddles = base->twiddles;
614            st->shift = 0;
615            while (nfft<<st->shift != base->nfft && st->shift < 32)
616               st->shift++;
617            /* FIXME: Report error and do proper cleanup */
618            if (st->shift>=32)
619               return NULL;
620         } else {
621            st->twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft);
622            compute_twiddles(st->twiddles, nfft);
623            st->shift = -1;
624         }
625         if (!kf_factor(nfft,st->factors))
626         {
627            kiss_fft_free(st);
628            return NULL;
629         }
630         
631         /* bitrev */
632         st->bitrev = (celt_int16*)KISS_FFT_MALLOC(sizeof(celt_int16)*nfft);
633         compute_bitrev_table(0, st->bitrev, 1,1, st->factors,st);
634     }
635     return st;
636 }
637
638 kiss_fft_cfg kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
639 {
640    return kiss_fft_alloc_twiddles(nfft, mem, lenmem, NULL);
641 }
642
643     
644 static void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
645 {
646     if (fin == fout) 
647     {
648        celt_fatal("In-place FFT not supported");
649     } else {
650        /* Bit-reverse the input */
651        int i;
652        for (i=0;i<st->nfft;i++)
653        {
654           fout[st->bitrev[i]] = fin[i];
655 #ifndef FIXED_POINT
656           fout[st->bitrev[i]].r *= st->scale;
657           fout[st->bitrev[i]].i *= st->scale;
658 #endif
659        }
660        kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
661     }
662 }
663
664 void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
665 {
666     kiss_fft_stride(cfg,fin,fout,1);
667 }
668
669 static void kiss_ifft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
670 {
671    if (fin == fout) 
672    {
673       celt_fatal("In-place FFT not supported");
674    } else {
675       /* Bit-reverse the input */
676       int i;
677       for (i=0;i<st->nfft;i++)
678          fout[st->bitrev[i]] = fin[i];
679       ki_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
680    }
681 }
682
683 void kiss_ifft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
684 {
685    kiss_ifft_stride(cfg,fin,fout,1);
686 }
687
688 void kiss_fft_free(kiss_fft_cfg cfg)
689 {
690    celt_free(cfg->bitrev);
691    if (cfg->shift < 0)
692       celt_free(cfg->twiddles);
693    celt_free(cfg);
694 }