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