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