Blame libcelt/kiss_fft.c

Packit 664db3
/*
Packit 664db3
Copyright (c) 2003-2004, Mark Borgerding
Packit 664db3
Lots of modifications by JMV
Packit 664db3
Copyright (c) 2005-2007, Jean-Marc Valin
Packit 664db3
Copyright (c) 2008,      Jean-Marc Valin, CSIRO
Packit 664db3
Packit 664db3
All rights reserved.
Packit 664db3
Packit 664db3
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
Packit 664db3
Packit 664db3
    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Packit 664db3
    * 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.
Packit 664db3
    * 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.
Packit 664db3
Packit 664db3
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.
Packit 664db3
*/
Packit 664db3
Packit 664db3
#ifndef SKIP_CONFIG_H
Packit 664db3
#  ifdef HAVE_CONFIG_H
Packit 664db3
#    include "config.h"
Packit 664db3
#  endif
Packit 664db3
#endif
Packit 664db3
Packit 664db3
#include "_kiss_fft_guts.h"
Packit 664db3
#include "arch.h"
Packit 664db3
#include "os_support.h"
Packit 664db3
#include "mathops.h"
Packit 664db3
#include "stack_alloc.h"
Packit 664db3
Packit 664db3
/* The guts header contains all the multiplication and addition macros that are defined for
Packit 664db3
   complex numbers.  It also delares the kf_ internal functions.
Packit 664db3
*/
Packit 664db3
Packit 664db3
static void kf_bfly2(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m,
Packit 664db3
                     int N,
Packit 664db3
                     int mm
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_fft_cpx * Fout2;
Packit 664db3
   kiss_twiddle_cpx * tw1;
Packit 664db3
   int i,j;
Packit 664db3
   kiss_fft_cpx * Fout_beg = Fout;
Packit 664db3
   for (i=0;i
Packit 664db3
   {
Packit 664db3
      Fout = Fout_beg + i*mm;
Packit 664db3
      Fout2 = Fout + m;
Packit 664db3
      tw1 = st->twiddles;
Packit 664db3
      for(j=0;j
Packit 664db3
      {
Packit 664db3
         kiss_fft_cpx t;
Packit 664db3
         Fout->r = SHR(Fout->r, 1);Fout->i = SHR(Fout->i, 1);
Packit 664db3
         Fout2->r = SHR(Fout2->r, 1);Fout2->i = SHR(Fout2->i, 1);
Packit 664db3
         C_MUL (t,  *Fout2 , *tw1);
Packit 664db3
         tw1 += fstride;
Packit 664db3
         C_SUB( *Fout2 ,  *Fout , t );
Packit 664db3
         C_ADDTO( *Fout ,  t );
Packit 664db3
         ++Fout2;
Packit 664db3
         ++Fout;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
static void ki_bfly2(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m,
Packit 664db3
                     int N,
Packit 664db3
                     int mm
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_fft_cpx * Fout2;
Packit 664db3
   kiss_twiddle_cpx * tw1;
Packit 664db3
   kiss_fft_cpx t;
Packit 664db3
   int i,j;
Packit 664db3
   kiss_fft_cpx * Fout_beg = Fout;
Packit 664db3
   for (i=0;i
Packit 664db3
   {
Packit 664db3
      Fout = Fout_beg + i*mm;
Packit 664db3
      Fout2 = Fout + m;
Packit 664db3
      tw1 = st->twiddles;
Packit 664db3
      for(j=0;j
Packit 664db3
      {
Packit 664db3
         C_MULC (t,  *Fout2 , *tw1);
Packit 664db3
         tw1 += fstride;
Packit 664db3
         C_SUB( *Fout2 ,  *Fout , t );
Packit 664db3
         C_ADDTO( *Fout ,  t );
Packit 664db3
         ++Fout2;
Packit 664db3
         ++Fout;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
static void kf_bfly4(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m,
Packit 664db3
                     int N,
Packit 664db3
                     int mm
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_twiddle_cpx *tw1,*tw2,*tw3;
Packit 664db3
   kiss_fft_cpx scratch[6];
Packit 664db3
   const size_t m2=2*m;
Packit 664db3
   const size_t m3=3*m;
Packit 664db3
   int i, j;
Packit 664db3
Packit 664db3
   kiss_fft_cpx * Fout_beg = Fout;
Packit 664db3
   for (i=0;i
Packit 664db3
   {
Packit 664db3
      Fout = Fout_beg + i*mm;
Packit 664db3
      tw3 = tw2 = tw1 = st->twiddles;
Packit 664db3
      for (j=0;j
Packit 664db3
      {
Packit 664db3
         C_MUL4(scratch[0],Fout[m] , *tw1 );
Packit 664db3
         C_MUL4(scratch[1],Fout[m2] , *tw2 );
Packit 664db3
         C_MUL4(scratch[2],Fout[m3] , *tw3 );
Packit 664db3
             
Packit 664db3
         Fout->r = PSHR(Fout->r, 2);
Packit 664db3
         Fout->i = PSHR(Fout->i, 2);
Packit 664db3
         C_SUB( scratch[5] , *Fout, scratch[1] );
Packit 664db3
         C_ADDTO(*Fout, scratch[1]);
Packit 664db3
         C_ADD( scratch[3] , scratch[0] , scratch[2] );
Packit 664db3
         C_SUB( scratch[4] , scratch[0] , scratch[2] );
Packit 664db3
         Fout[m2].r = PSHR(Fout[m2].r, 2);
Packit 664db3
         Fout[m2].i = PSHR(Fout[m2].i, 2);
Packit 664db3
         C_SUB( Fout[m2], *Fout, scratch[3] );
Packit 664db3
         tw1 += fstride;
Packit 664db3
         tw2 += fstride*2;
Packit 664db3
         tw3 += fstride*3;
Packit 664db3
         C_ADDTO( *Fout , scratch[3] );
Packit 664db3
             
Packit 664db3
         Fout[m].r = scratch[5].r + scratch[4].i;
Packit 664db3
         Fout[m].i = scratch[5].i - scratch[4].r;
Packit 664db3
         Fout[m3].r = scratch[5].r - scratch[4].i;
Packit 664db3
         Fout[m3].i = scratch[5].i + scratch[4].r;
Packit 664db3
         ++Fout;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
static void ki_bfly4(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m,
Packit 664db3
                     int N,
Packit 664db3
                     int mm
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_twiddle_cpx *tw1,*tw2,*tw3;
Packit 664db3
   kiss_fft_cpx scratch[6];
Packit 664db3
   const size_t m2=2*m;
Packit 664db3
   const size_t m3=3*m;
Packit 664db3
   int i, j;
Packit 664db3
Packit 664db3
   kiss_fft_cpx * Fout_beg = Fout;
Packit 664db3
   for (i=0;i
Packit 664db3
   {
Packit 664db3
      Fout = Fout_beg + i*mm;
Packit 664db3
      tw3 = tw2 = tw1 = st->twiddles;
Packit 664db3
      for (j=0;j
Packit 664db3
      {
Packit 664db3
         C_MULC(scratch[0],Fout[m] , *tw1 );
Packit 664db3
         C_MULC(scratch[1],Fout[m2] , *tw2 );
Packit 664db3
         C_MULC(scratch[2],Fout[m3] , *tw3 );
Packit 664db3
             
Packit 664db3
         C_SUB( scratch[5] , *Fout, scratch[1] );
Packit 664db3
         C_ADDTO(*Fout, scratch[1]);
Packit 664db3
         C_ADD( scratch[3] , scratch[0] , scratch[2] );
Packit 664db3
         C_SUB( scratch[4] , scratch[0] , scratch[2] );
Packit 664db3
         C_SUB( Fout[m2], *Fout, scratch[3] );
Packit 664db3
         tw1 += fstride;
Packit 664db3
         tw2 += fstride*2;
Packit 664db3
         tw3 += fstride*3;
Packit 664db3
         C_ADDTO( *Fout , scratch[3] );
Packit 664db3
             
Packit 664db3
         Fout[m].r = scratch[5].r - scratch[4].i;
Packit 664db3
         Fout[m].i = scratch[5].i + scratch[4].r;
Packit 664db3
         Fout[m3].r = scratch[5].r + scratch[4].i;
Packit 664db3
         Fout[m3].i = scratch[5].i - scratch[4].r;
Packit 664db3
         ++Fout;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
#ifndef RADIX_TWO_ONLY
Packit 664db3
Packit 664db3
static void kf_bfly3(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     size_t m
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   size_t k=m;
Packit 664db3
   const size_t m2 = 2*m;
Packit 664db3
   kiss_twiddle_cpx *tw1,*tw2;
Packit 664db3
   kiss_fft_cpx scratch[5];
Packit 664db3
   kiss_twiddle_cpx epi3;
Packit 664db3
   epi3 = st->twiddles[fstride*m];
Packit 664db3
Packit 664db3
   tw1=tw2=st->twiddles;
Packit 664db3
   do{
Packit 664db3
      C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
Packit 664db3
Packit 664db3
      C_MUL(scratch[1],Fout[m] , *tw1);
Packit 664db3
      C_MUL(scratch[2],Fout[m2] , *tw2);
Packit 664db3
Packit 664db3
      C_ADD(scratch[3],scratch[1],scratch[2]);
Packit 664db3
      C_SUB(scratch[0],scratch[1],scratch[2]);
Packit 664db3
      tw1 += fstride;
Packit 664db3
      tw2 += fstride*2;
Packit 664db3
Packit 664db3
      Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
Packit 664db3
      Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
Packit 664db3
Packit 664db3
      C_MULBYSCALAR( scratch[0] , epi3.i );
Packit 664db3
Packit 664db3
      C_ADDTO(*Fout,scratch[3]);
Packit 664db3
Packit 664db3
      Fout[m2].r = Fout[m].r + scratch[0].i;
Packit 664db3
      Fout[m2].i = Fout[m].i - scratch[0].r;
Packit 664db3
Packit 664db3
      Fout[m].r -= scratch[0].i;
Packit 664db3
      Fout[m].i += scratch[0].r;
Packit 664db3
Packit 664db3
      ++Fout;
Packit 664db3
   }while(--k);
Packit 664db3
}
Packit 664db3
Packit 664db3
static void ki_bfly3(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     size_t m
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   size_t k=m;
Packit 664db3
   const size_t m2 = 2*m;
Packit 664db3
   kiss_twiddle_cpx *tw1,*tw2;
Packit 664db3
   kiss_fft_cpx scratch[5];
Packit 664db3
   kiss_twiddle_cpx epi3;
Packit 664db3
   epi3 = st->twiddles[fstride*m];
Packit 664db3
Packit 664db3
   tw1=tw2=st->twiddles;
Packit 664db3
   do{
Packit 664db3
Packit 664db3
      C_MULC(scratch[1],Fout[m] , *tw1);
Packit 664db3
      C_MULC(scratch[2],Fout[m2] , *tw2);
Packit 664db3
Packit 664db3
      C_ADD(scratch[3],scratch[1],scratch[2]);
Packit 664db3
      C_SUB(scratch[0],scratch[1],scratch[2]);
Packit 664db3
      tw1 += fstride;
Packit 664db3
      tw2 += fstride*2;
Packit 664db3
Packit 664db3
      Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
Packit 664db3
      Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
Packit 664db3
Packit 664db3
      C_MULBYSCALAR( scratch[0] , -epi3.i );
Packit 664db3
Packit 664db3
      C_ADDTO(*Fout,scratch[3]);
Packit 664db3
Packit 664db3
      Fout[m2].r = Fout[m].r + scratch[0].i;
Packit 664db3
      Fout[m2].i = Fout[m].i - scratch[0].r;
Packit 664db3
Packit 664db3
      Fout[m].r -= scratch[0].i;
Packit 664db3
      Fout[m].i += scratch[0].r;
Packit 664db3
Packit 664db3
      ++Fout;
Packit 664db3
   }while(--k);
Packit 664db3
}
Packit 664db3
Packit 664db3
Packit 664db3
static void kf_bfly5(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
Packit 664db3
   int u;
Packit 664db3
   kiss_fft_cpx scratch[13];
Packit 664db3
   kiss_twiddle_cpx * twiddles = st->twiddles;
Packit 664db3
   kiss_twiddle_cpx *tw;
Packit 664db3
   kiss_twiddle_cpx ya,yb;
Packit 664db3
   ya = twiddles[fstride*m];
Packit 664db3
   yb = twiddles[fstride*2*m];
Packit 664db3
Packit 664db3
   Fout0=Fout;
Packit 664db3
   Fout1=Fout0+m;
Packit 664db3
   Fout2=Fout0+2*m;
Packit 664db3
   Fout3=Fout0+3*m;
Packit 664db3
   Fout4=Fout0+4*m;
Packit 664db3
Packit 664db3
   tw=st->twiddles;
Packit 664db3
   for ( u=0; u
Packit 664db3
      C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
Packit 664db3
      scratch[0] = *Fout0;
Packit 664db3
Packit 664db3
      C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
Packit 664db3
      C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
Packit 664db3
      C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
Packit 664db3
      C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
Packit 664db3
Packit 664db3
      C_ADD( scratch[7],scratch[1],scratch[4]);
Packit 664db3
      C_SUB( scratch[10],scratch[1],scratch[4]);
Packit 664db3
      C_ADD( scratch[8],scratch[2],scratch[3]);
Packit 664db3
      C_SUB( scratch[9],scratch[2],scratch[3]);
Packit 664db3
Packit 664db3
      Fout0->r += scratch[7].r + scratch[8].r;
Packit 664db3
      Fout0->i += scratch[7].i + scratch[8].i;
Packit 664db3
Packit 664db3
      scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
Packit 664db3
      scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
Packit 664db3
Packit 664db3
      scratch[6].r =  S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
Packit 664db3
      scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
Packit 664db3
Packit 664db3
      C_SUB(*Fout1,scratch[5],scratch[6]);
Packit 664db3
      C_ADD(*Fout4,scratch[5],scratch[6]);
Packit 664db3
Packit 664db3
      scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
Packit 664db3
      scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
Packit 664db3
      scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
Packit 664db3
      scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
Packit 664db3
Packit 664db3
      C_ADD(*Fout2,scratch[11],scratch[12]);
Packit 664db3
      C_SUB(*Fout3,scratch[11],scratch[12]);
Packit 664db3
Packit 664db3
      ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
static void ki_bfly5(
Packit 664db3
                     kiss_fft_cpx * Fout,
Packit 664db3
                     const size_t fstride,
Packit 664db3
                     const kiss_fft_cfg st,
Packit 664db3
                     int m
Packit 664db3
                    )
Packit 664db3
{
Packit 664db3
   kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
Packit 664db3
   int u;
Packit 664db3
   kiss_fft_cpx scratch[13];
Packit 664db3
   kiss_twiddle_cpx * twiddles = st->twiddles;
Packit 664db3
   kiss_twiddle_cpx *tw;
Packit 664db3
   kiss_twiddle_cpx ya,yb;
Packit 664db3
   ya = twiddles[fstride*m];
Packit 664db3
   yb = twiddles[fstride*2*m];
Packit 664db3
Packit 664db3
   Fout0=Fout;
Packit 664db3
   Fout1=Fout0+m;
Packit 664db3
   Fout2=Fout0+2*m;
Packit 664db3
   Fout3=Fout0+3*m;
Packit 664db3
   Fout4=Fout0+4*m;
Packit 664db3
Packit 664db3
   tw=st->twiddles;
Packit 664db3
   for ( u=0; u
Packit 664db3
      scratch[0] = *Fout0;
Packit 664db3
Packit 664db3
      C_MULC(scratch[1] ,*Fout1, tw[u*fstride]);
Packit 664db3
      C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]);
Packit 664db3
      C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]);
Packit 664db3
      C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]);
Packit 664db3
Packit 664db3
      C_ADD( scratch[7],scratch[1],scratch[4]);
Packit 664db3
      C_SUB( scratch[10],scratch[1],scratch[4]);
Packit 664db3
      C_ADD( scratch[8],scratch[2],scratch[3]);
Packit 664db3
      C_SUB( scratch[9],scratch[2],scratch[3]);
Packit 664db3
Packit 664db3
      Fout0->r += scratch[7].r + scratch[8].r;
Packit 664db3
      Fout0->i += scratch[7].i + scratch[8].i;
Packit 664db3
Packit 664db3
      scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
Packit 664db3
      scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
Packit 664db3
Packit 664db3
      scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i);
Packit 664db3
      scratch[6].i =  S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i);
Packit 664db3
Packit 664db3
      C_SUB(*Fout1,scratch[5],scratch[6]);
Packit 664db3
      C_ADD(*Fout4,scratch[5],scratch[6]);
Packit 664db3
Packit 664db3
      scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
Packit 664db3
      scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
Packit 664db3
      scratch[12].r =  S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i);
Packit 664db3
      scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i);
Packit 664db3
Packit 664db3
      C_ADD(*Fout2,scratch[11],scratch[12]);
Packit 664db3
      C_SUB(*Fout3,scratch[11],scratch[12]);
Packit 664db3
Packit 664db3
      ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
/* perform the butterfly for one stage of a mixed radix FFT */
Packit 664db3
static void kf_bfly_generic(
Packit 664db3
                            kiss_fft_cpx * Fout,
Packit 664db3
                            const size_t fstride,
Packit 664db3
                            const kiss_fft_cfg st,
Packit 664db3
                            int m,
Packit 664db3
                            int p
Packit 664db3
                           )
Packit 664db3
{
Packit 664db3
   int u,k,q1,q;
Packit 664db3
   kiss_twiddle_cpx * twiddles = st->twiddles;
Packit 664db3
   kiss_fft_cpx t;
Packit 664db3
   VARDECL(kiss_fft_cpx, scratchbuf);
Packit 664db3
   int Norig = st->nfft;
Packit 664db3
   ALLOC(scratchbuf, p, kiss_fft_cpx);
Packit 664db3
Packit 664db3
   for ( u=0; u
Packit 664db3
      k=u;
Packit 664db3
      for ( q1=0 ; q1
Packit 664db3
         scratchbuf[q1] = Fout[ k  ];
Packit 664db3
         C_FIXDIV(scratchbuf[q1],p);
Packit 664db3
         k += m;
Packit 664db3
      }
Packit 664db3
Packit 664db3
      k=u;
Packit 664db3
      for ( q1=0 ; q1
Packit 664db3
         int twidx=0;
Packit 664db3
         Fout[ k ] = scratchbuf[0];
Packit 664db3
         for (q=1;q
Packit 664db3
            twidx += fstride * k;
Packit 664db3
            if (twidx>=Norig) twidx-=Norig;
Packit 664db3
            C_MUL(t,scratchbuf[q] , twiddles[twidx] );
Packit 664db3
            C_ADDTO( Fout[ k ] ,t);
Packit 664db3
         }
Packit 664db3
         k += m;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
static void ki_bfly_generic(
Packit 664db3
                            kiss_fft_cpx * Fout,
Packit 664db3
                            const size_t fstride,
Packit 664db3
                            const kiss_fft_cfg st,
Packit 664db3
                            int m,
Packit 664db3
                            int p
Packit 664db3
                           )
Packit 664db3
{
Packit 664db3
   int u,k,q1,q;
Packit 664db3
   kiss_twiddle_cpx * twiddles = st->twiddles;
Packit 664db3
   kiss_fft_cpx t;
Packit 664db3
   VARDECL(kiss_fft_cpx, scratchbuf);
Packit 664db3
   int Norig = st->nfft;
Packit 664db3
   ALLOC(scratchbuf, p, kiss_fft_cpx);
Packit 664db3
Packit 664db3
   for ( u=0; u
Packit 664db3
      k=u;
Packit 664db3
      for ( q1=0 ; q1
Packit 664db3
         scratchbuf[q1] = Fout[ k  ];
Packit 664db3
         k += m;
Packit 664db3
      }
Packit 664db3
Packit 664db3
      k=u;
Packit 664db3
      for ( q1=0 ; q1
Packit 664db3
         int twidx=0;
Packit 664db3
         Fout[ k ] = scratchbuf[0];
Packit 664db3
         for (q=1;q
Packit 664db3
            twidx += fstride * k;
Packit 664db3
            if (twidx>=Norig) twidx-=Norig;
Packit 664db3
            C_MULC(t,scratchbuf[q] , twiddles[twidx] );
Packit 664db3
            C_ADDTO( Fout[ k ] ,t);
Packit 664db3
         }
Packit 664db3
         k += m;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
#endif
Packit 664db3
Packit 664db3
static
Packit 664db3
void compute_bitrev_table(
Packit 664db3
         int Fout,
Packit 664db3
         int *f,
Packit 664db3
         const size_t fstride,
Packit 664db3
         int in_stride,
Packit 664db3
         int * factors,
Packit 664db3
         const kiss_fft_cfg st
Packit 664db3
            )
Packit 664db3
{
Packit 664db3
   const int p=*factors++; /* the radix  */
Packit 664db3
   const int m=*factors++; /* stage's fft length/p */
Packit 664db3
   
Packit 664db3
    /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
Packit 664db3
   if (m==1)
Packit 664db3
   {
Packit 664db3
      int j;
Packit 664db3
      for (j=0;j
Packit 664db3
      {
Packit 664db3
         *f = Fout+j;
Packit 664db3
         f += fstride*in_stride;
Packit 664db3
      }
Packit 664db3
   } else {
Packit 664db3
      int j;
Packit 664db3
      for (j=0;j
Packit 664db3
      {
Packit 664db3
         compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
Packit 664db3
         f += fstride*in_stride;
Packit 664db3
         Fout += m;
Packit 664db3
      }
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
Packit 664db3
void kf_work(
Packit 664db3
        kiss_fft_cpx * Fout,
Packit 664db3
        const kiss_fft_cpx * f,
Packit 664db3
        const size_t fstride,
Packit 664db3
        int in_stride,
Packit 664db3
        int * factors,
Packit 664db3
        const kiss_fft_cfg st,
Packit 664db3
        int N,
Packit 664db3
        int s2,
Packit 664db3
        int m2
Packit 664db3
        )
Packit 664db3
{
Packit 664db3
#ifndef RADIX_TWO_ONLY
Packit 664db3
    int i;
Packit 664db3
    kiss_fft_cpx * Fout_beg=Fout;
Packit 664db3
#endif
Packit 664db3
    const int p=*factors++; /* the radix  */
Packit 664db3
    const int m=*factors++; /* stage's fft length/p */
Packit 664db3
    /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
Packit 664db3
    if (m!=1) 
Packit 664db3
        kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
Packit 664db3
Packit 664db3
    switch (p) {
Packit 664db3
        case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break;
Packit 664db3
        case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break;
Packit 664db3
#ifndef RADIX_TWO_ONLY
Packit 664db3
        case 3: for (i=0;i
Packit 664db3
        case 5: for (i=0;i
Packit 664db3
        default: for (i=0;i
Packit 664db3
#else
Packit 664db3
       default: celt_fatal("kiss_fft: only powers of two enabled");
Packit 664db3
#endif
Packit 664db3
    }    
Packit 664db3
}
Packit 664db3
Packit 664db3
Packit 664db3
void ki_work(
Packit 664db3
             kiss_fft_cpx * Fout,
Packit 664db3
             const kiss_fft_cpx * f,
Packit 664db3
             const size_t fstride,
Packit 664db3
             int in_stride,
Packit 664db3
             int * factors,
Packit 664db3
             const kiss_fft_cfg st,
Packit 664db3
             int N,
Packit 664db3
             int s2,
Packit 664db3
             int m2
Packit 664db3
            )
Packit 664db3
{
Packit 664db3
#ifndef RADIX_TWO_ONLY
Packit 664db3
   int i;
Packit 664db3
   kiss_fft_cpx * Fout_beg=Fout;
Packit 664db3
#endif
Packit 664db3
   const int p=*factors++; /* the radix  */
Packit 664db3
   const int m=*factors++; /* stage's fft length/p */
Packit 664db3
   /*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/
Packit 664db3
   if (m!=1) 
Packit 664db3
      ki_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m);
Packit 664db3
Packit 664db3
   switch (p) {
Packit 664db3
      case 2: ki_bfly2(Fout,fstride,st,m, N, m2); break;
Packit 664db3
      case 4: ki_bfly4(Fout,fstride,st,m, N, m2); break;
Packit 664db3
#ifndef RADIX_TWO_ONLY
Packit 664db3
      case 3: for (i=0;i
Packit 664db3
      case 5: for (i=0;i
Packit 664db3
      default: for (i=0;i
Packit 664db3
#else
Packit 664db3
      default: celt_fatal("kiss_fft: only powers of two enabled");
Packit 664db3
#endif
Packit 664db3
   }    
Packit 664db3
}
Packit 664db3
Packit 664db3
/*  facbuf is populated by p1,m1,p2,m2, ...
Packit 664db3
    where 
Packit 664db3
    p[i] * m[i] = m[i-1]
Packit 664db3
    m0 = n                  */
Packit 664db3
static 
Packit 664db3
void kf_factor(int n,int * facbuf)
Packit 664db3
{
Packit 664db3
    int p=4;
Packit 664db3
Packit 664db3
    /*factor out powers of 4, powers of 2, then any remaining primes */
Packit 664db3
    do {
Packit 664db3
        while (n % p) {
Packit 664db3
            switch (p) {
Packit 664db3
                case 4: p = 2; break;
Packit 664db3
                case 2: p = 3; break;
Packit 664db3
                default: p += 2; break;
Packit 664db3
            }
Packit 664db3
            if (p>32000 || (celt_int32_t)p*(celt_int32_t)p > n)
Packit 664db3
                p = n;          /* no more factors, skip to end */
Packit 664db3
        }
Packit 664db3
        n /= p;
Packit 664db3
        *facbuf++ = p;
Packit 664db3
        *facbuf++ = n;
Packit 664db3
    } while (n > 1);
Packit 664db3
}
Packit 664db3
/*
Packit 664db3
 *
Packit 664db3
 * User-callable function to allocate all necessary storage space for the fft.
Packit 664db3
 *
Packit 664db3
 * The return value is a contiguous block of memory, allocated with malloc.  As such,
Packit 664db3
 * It can be freed with free(), rather than a kiss_fft-specific function.
Packit 664db3
 * */
Packit 664db3
kiss_fft_cfg kiss_fft_alloc(int nfft,void * mem,size_t * lenmem )
Packit 664db3
{
Packit 664db3
    kiss_fft_cfg st=NULL;
Packit 664db3
    size_t memneeded = sizeof(struct kiss_fft_state)
Packit 664db3
          + sizeof(kiss_twiddle_cpx)*(nfft-1) + sizeof(int)*nfft; /* twiddle factors*/
Packit 664db3
Packit 664db3
    if ( lenmem==NULL ) {
Packit 664db3
        st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
Packit 664db3
    }else{
Packit 664db3
        if (mem != NULL && *lenmem >= memneeded)
Packit 664db3
            st = (kiss_fft_cfg)mem;
Packit 664db3
        *lenmem = memneeded;
Packit 664db3
    }
Packit 664db3
    if (st) {
Packit 664db3
        int i;
Packit 664db3
        st->nfft=nfft;
Packit 664db3
#ifndef FIXED_POINT
Packit 664db3
        st->scale = 1./nfft;
Packit 664db3
#endif
Packit 664db3
#if defined(FIXED_POINT) && (!defined(DOUBLE_PRECISION) || defined(MIXED_PRECISION))
Packit 664db3
        for (i=0;i
Packit 664db3
            celt_word32_t phase = -i;
Packit 664db3
            kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft));
Packit 664db3
        }
Packit 664db3
#else
Packit 664db3
        for (i=0;i
Packit 664db3
           const double pi=3.14159265358979323846264338327;
Packit 664db3
           double phase = ( -2*pi /nfft ) * i;
Packit 664db3
           kf_cexp(st->twiddles+i, phase );
Packit 664db3
        }
Packit 664db3
#endif
Packit 664db3
        kf_factor(nfft,st->factors);
Packit 664db3
        
Packit 664db3
        /* bitrev */
Packit 664db3
        st->bitrev = (int*)((char*)st + memneeded - sizeof(int)*nfft);
Packit 664db3
        compute_bitrev_table(0, st->bitrev, 1,1, st->factors,st);
Packit 664db3
    }
Packit 664db3
    return st;
Packit 664db3
}
Packit 664db3
Packit 664db3
Packit 664db3
Packit 664db3
    
Packit 664db3
void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
Packit 664db3
{
Packit 664db3
    if (fin == fout) 
Packit 664db3
    {
Packit 664db3
       celt_fatal("In-place FFT not supported");
Packit 664db3
    } else {
Packit 664db3
       /* Bit-reverse the input */
Packit 664db3
       int i;
Packit 664db3
       for (i=0;i<st->nfft;i++)
Packit 664db3
       {
Packit 664db3
          fout[st->bitrev[i]] = fin[i];
Packit 664db3
#ifndef FIXED_POINT
Packit 664db3
          fout[st->bitrev[i]].r *= st->scale;
Packit 664db3
          fout[st->bitrev[i]].i *= st->scale;
Packit 664db3
#endif
Packit 664db3
       }
Packit 664db3
       kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
Packit 664db3
    }
Packit 664db3
}
Packit 664db3
Packit 664db3
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
Packit 664db3
{
Packit 664db3
    kiss_fft_stride(cfg,fin,fout,1);
Packit 664db3
}
Packit 664db3
Packit 664db3
void kiss_ifft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
Packit 664db3
{
Packit 664db3
   if (fin == fout) 
Packit 664db3
   {
Packit 664db3
      celt_fatal("In-place FFT not supported");
Packit 664db3
   } else {
Packit 664db3
      /* Bit-reverse the input */
Packit 664db3
      int i;
Packit 664db3
      for (i=0;i<st->nfft;i++)
Packit 664db3
         fout[st->bitrev[i]] = fin[i];
Packit 664db3
      ki_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1);
Packit 664db3
   }
Packit 664db3
}
Packit 664db3
Packit 664db3
void kiss_ifft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
Packit 664db3
{
Packit 664db3
   kiss_ifft_stride(cfg,fin,fout,1);
Packit 664db3
}
Packit 664db3