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