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_state *st,
  38                      int m,
  39                      int N,
  40                      int mm
  41                     )
  42 {
  43    kiss_fft_cpx * Fout2;
  44    const 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_state *st,
  71                      int m,
  72                      int N,
  73                      int mm
  74                     )
  75 {
  76    kiss_fft_cpx * Fout2;
  77    const 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_state *st,
 102                      int m,
 103                      int N,
 104                      int mm
 105                     )
 106 {
 107    const 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_state *st,
 151                      int m,
 152                      int N,
 153                      int mm
 154                     )
 155 {
 156    const 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_state *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    const 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_state *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    const 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_state *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    const kiss_twiddle_cpx * twiddles = st->twiddles;
 310    const 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_state *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    const kiss_twiddle_cpx * twiddles = st->twiddles;
 379    const 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 void kf_work(
 437         kiss_fft_cpx * Fout,
 438         const kiss_fft_cpx * f,
 439         size_t fstride,
 440         int in_stride,
 441         const celt_int16 * factors,
 442         const kiss_fft_state *st,
 443         int N,
 444         int m2
 445         )
 446 {
 447     const int p=*factors++; /* the radix  */
 448     const int m=*factors++; /* stage's fft length/p */
 449     /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
 450     if (m!=1)
 451         kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, m);
 452
 453     /* Compensate for longer twiddles table (when sharing) */
 454     if (st->shift>0)
 455        fstride <<= st->shift;
 456     switch (p) {
 457         case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break;
 458         case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
 459 #ifndef RADIX_TWO_ONLY
 460         case 3: kf_bfly3(Fout,fstride,st,m, N, m2); break;
 461         case 5: kf_bfly5(Fout,fstride,st,m, N, m2); break;
 462 #endif
 463     }
 464 }
 465
 466
 467 static void ki_work(
 468              kiss_fft_cpx * Fout,
 469              const kiss_fft_cpx * f,
 470              size_t fstride,
 471              int in_stride,
 472              const celt_int16 * factors,
 473              const kiss_fft_state *st,
 474              int N,
 475              int m2
 476             )
 477 {
 478    const int p=*factors++; /* the radix  */
 479    const int m=*factors++; /* stage's fft length/p */
 480    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
 481    if (m!=1)
 482       ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, m);
 483
 484    /* Compensate for longer twiddles table (when sharing) */
 485    if (st->shift>0)
 486       fstride <<= st->shift;
 487    switch (p) {
 488       case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
 489       case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
 490 #ifndef RADIX_TWO_ONLY
 491       case 3: ki_bfly3(Fout,fstride,st,m, N, m2); break;
 492       case 5: ki_bfly5(Fout,fstride,st,m, N, m2); break;
 493 #endif
 494    }
 495 }
 496
 497
 498 #ifdef CUSTOM_MODES
 499
 500 static
 501 void compute_bitrev_table(
 502          int Fout,
 503          celt_int16 *f,
 504          const size_t fstride,
 505          int in_stride,
 506          celt_int16 * factors,
 507          const kiss_fft_state *st
 508             )
 509 {
 510    const int p=*factors++; /* the radix  */
 511    const int m=*factors++; /* stage's fft length/p */
 512
 513     /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
 514    if (m==1)
 515    {
 516       int j;
 517       for (j=0;j<p;j++)
 518       {
 519          *f = Fout+j;
 520          f += fstride*in_stride;
 521       }
 522    } else {
 523       int j;
 524       for (j=0;j<p;j++)
 525       {
 526          compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
 527          f += fstride*in_stride;
 528          Fout += m;
 529       }
 530    }
 531 }
 532
 533
 534 /*  facbuf is populated by p1,m1,p2,m2, ...
 535     where
 536     p[i] * m[i] = m[i-1]
 537     m0 = n                  */
 538 static
 539 int kf_factor(int n,celt_int16 * facbuf)
 540 {
 541     int p=4;
 542
 543     /*factor out powers of 4, powers of 2, then any remaining primes */
 544     do {
 545         while (n % p) {
 546             switch (p) {
 547                 case 4: p = 2; break;
 548                 case 2: p = 3; break;
 549                 default: p += 2; break;
 550             }
 551             if (p>32000 || (celt_int32)p*(celt_int32)p > n)
 552                 p = n;          /* no more factors, skip to end */
 553         }
 554         n /= p;
 555 #ifdef RADIX_TWO_ONLY
 556         if (p!=2 && p != 4)
 557 #else
 558         if (p>5)
 559 #endif
 560         {
 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 #ifdef FIXED_POINT
 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_state *kiss_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem,  const kiss_fft_state *base)
 595 {
 596     kiss_fft_state *st=NULL;
 597     size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/
 598
 599     if ( lenmem==NULL ) {
 600         st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded );
 601     }else{
 602         if (mem != NULL && *lenmem >= memneeded)
 603             st = (kiss_fft_state*)mem;
 604         *lenmem = memneeded;
 605     }
 606     if (st) {
 607         celt_int16 *bitrev;
 608         kiss_twiddle_cpx *twiddles;
 609
 610         st->nfft=nfft;
 611 #ifndef FIXED_POINT
 612         st->scale = 1./nfft;
 613 #endif
 614         if (base != NULL)
 615         {
 616            st->twiddles = base->twiddles;
 617            st->shift = 0;
 618            while (nfft<<st->shift != base->nfft && st->shift < 32)
 619               st->shift++;
 620            /* FIXME: Report error and do proper cleanup */
 621            if (st->shift>=32)
 622               return NULL;
 623         } else {
 624            st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft);
 625            compute_twiddles(twiddles, nfft);
 626            st->shift = -1;
 627         }
 628         if (!kf_factor(nfft,st->factors))
 629         {
 630            kiss_fft_free(st);
 631            return NULL;
 632         }
 633
 634         /* bitrev */
 635         st->bitrev = bitrev = (celt_int16*)KISS_FFT_MALLOC(sizeof(celt_int16)*nfft);
 636         compute_bitrev_table(0, bitrev, 1,1, st->factors,st);
 637     }
 638     return st;
 639 }
 640
 641 kiss_fft_state *kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
 642 {
 643    return kiss_fft_alloc_twiddles(nfft, mem, lenmem, NULL);
 644 }
 645
 646 void kiss_fft_free(const kiss_fft_state *cfg)
 647 {
 648    celt_free((celt_int16*)cfg->bitrev);
 649    if (cfg->shift < 0)
 650       celt_free((kiss_twiddle_cpx*)cfg->twiddles);
 651    celt_free((kiss_fft_state*)cfg);
 652 }
 653
 654 #endif /* CUSTOM_MODES */
 655
 656 static void kiss_fft_stride(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
 657 {
 658     int i;
 659     celt_assert2 (fin != fout, "In-place FFT not supported");
 660     /* Bit-reverse the input */
 661     for (i=0;i<st->nfft;i++)
 662     {
 663        fout[st->bitrev[i]] = fin[i];
 664 #ifndef FIXED_POINT
 665        fout[st->bitrev[i]].r *= st->scale;
 666        fout[st->bitrev[i]].i *= st->scale;
 667 #endif
 668     }
 669     kf_work( fout, fin, 1,in_stride, st->factors,st, 1, 1);
 670 }
 671
 672 void kiss_fft(const kiss_fft_state *cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
 673 {
 674     kiss_fft_stride(cfg,fin,fout,1);
 675 }
 676
 677 static void kiss_ifft_stride(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
 678 {
 679    int i;
 680    celt_assert2 (fin != fout, "In-place FFT not supported");
 681    /* Bit-reverse the input */
 682    for (i=0;i<st->nfft;i++)
 683       fout[st->bitrev[i]] = fin[i];
 684    ki_work( fout, fin, 1,in_stride, st->factors,st, 1, 1);
 685 }
 686
 687 void kiss_ifft(const kiss_fft_state *cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
 688 {
 689    kiss_ifft_stride(cfg,fin,fout,1);
 690 }
 691