1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/celt/kiss_fft.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,719 @@
1.4 +/*Copyright (c) 2003-2004, Mark Borgerding
1.5 + Lots of modifications by Jean-Marc Valin
1.6 + Copyright (c) 2005-2007, Xiph.Org Foundation
1.7 + Copyright (c) 2008, Xiph.Org Foundation, CSIRO
1.8 +
1.9 + All rights reserved.
1.10 +
1.11 + Redistribution and use in source and binary forms, with or without
1.12 + modification, are permitted provided that the following conditions are met:
1.13 +
1.14 + * Redistributions of source code must retain the above copyright notice,
1.15 + this list of conditions and the following disclaimer.
1.16 + * Redistributions in binary form must reproduce the above copyright notice,
1.17 + this list of conditions and the following disclaimer in the
1.18 + documentation and/or other materials provided with the distribution.
1.19 +
1.20 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1.21 + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1.22 + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1.23 + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
1.24 + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
1.25 + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
1.26 + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
1.27 + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
1.28 + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
1.29 + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
1.30 + POSSIBILITY OF SUCH DAMAGE.*/
1.31 +
1.32 +/* This code is originally from Mark Borgerding's KISS-FFT but has been
1.33 + heavily modified to better suit Opus */
1.34 +
1.35 +#ifndef SKIP_CONFIG_H
1.36 +# ifdef HAVE_CONFIG_H
1.37 +# include "config.h"
1.38 +# endif
1.39 +#endif
1.40 +
1.41 +#include "_kiss_fft_guts.h"
1.42 +#include "arch.h"
1.43 +#include "os_support.h"
1.44 +#include "mathops.h"
1.45 +#include "stack_alloc.h"
1.46 +
1.47 +/* The guts header contains all the multiplication and addition macros that are defined for
1.48 + complex numbers. It also delares the kf_ internal functions.
1.49 +*/
1.50 +
1.51 +static void kf_bfly2(
1.52 + kiss_fft_cpx * Fout,
1.53 + const size_t fstride,
1.54 + const kiss_fft_state *st,
1.55 + int m,
1.56 + int N,
1.57 + int mm
1.58 + )
1.59 +{
1.60 + kiss_fft_cpx * Fout2;
1.61 + const kiss_twiddle_cpx * tw1;
1.62 + int i,j;
1.63 + kiss_fft_cpx * Fout_beg = Fout;
1.64 + for (i=0;i<N;i++)
1.65 + {
1.66 + Fout = Fout_beg + i*mm;
1.67 + Fout2 = Fout + m;
1.68 + tw1 = st->twiddles;
1.69 + for(j=0;j<m;j++)
1.70 + {
1.71 + kiss_fft_cpx t;
1.72 + Fout->r = SHR32(Fout->r, 1);Fout->i = SHR32(Fout->i, 1);
1.73 + Fout2->r = SHR32(Fout2->r, 1);Fout2->i = SHR32(Fout2->i, 1);
1.74 + C_MUL (t, *Fout2 , *tw1);
1.75 + tw1 += fstride;
1.76 + C_SUB( *Fout2 , *Fout , t );
1.77 + C_ADDTO( *Fout , t );
1.78 + ++Fout2;
1.79 + ++Fout;
1.80 + }
1.81 + }
1.82 +}
1.83 +
1.84 +static void ki_bfly2(
1.85 + kiss_fft_cpx * Fout,
1.86 + const size_t fstride,
1.87 + const kiss_fft_state *st,
1.88 + int m,
1.89 + int N,
1.90 + int mm
1.91 + )
1.92 +{
1.93 + kiss_fft_cpx * Fout2;
1.94 + const kiss_twiddle_cpx * tw1;
1.95 + kiss_fft_cpx t;
1.96 + int i,j;
1.97 + kiss_fft_cpx * Fout_beg = Fout;
1.98 + for (i=0;i<N;i++)
1.99 + {
1.100 + Fout = Fout_beg + i*mm;
1.101 + Fout2 = Fout + m;
1.102 + tw1 = st->twiddles;
1.103 + for(j=0;j<m;j++)
1.104 + {
1.105 + C_MULC (t, *Fout2 , *tw1);
1.106 + tw1 += fstride;
1.107 + C_SUB( *Fout2 , *Fout , t );
1.108 + C_ADDTO( *Fout , t );
1.109 + ++Fout2;
1.110 + ++Fout;
1.111 + }
1.112 + }
1.113 +}
1.114 +
1.115 +static void kf_bfly4(
1.116 + kiss_fft_cpx * Fout,
1.117 + const size_t fstride,
1.118 + const kiss_fft_state *st,
1.119 + int m,
1.120 + int N,
1.121 + int mm
1.122 + )
1.123 +{
1.124 + const kiss_twiddle_cpx *tw1,*tw2,*tw3;
1.125 + kiss_fft_cpx scratch[6];
1.126 + const size_t m2=2*m;
1.127 + const size_t m3=3*m;
1.128 + int i, j;
1.129 +
1.130 + kiss_fft_cpx * Fout_beg = Fout;
1.131 + for (i=0;i<N;i++)
1.132 + {
1.133 + Fout = Fout_beg + i*mm;
1.134 + tw3 = tw2 = tw1 = st->twiddles;
1.135 + for (j=0;j<m;j++)
1.136 + {
1.137 + C_MUL4(scratch[0],Fout[m] , *tw1 );
1.138 + C_MUL4(scratch[1],Fout[m2] , *tw2 );
1.139 + C_MUL4(scratch[2],Fout[m3] , *tw3 );
1.140 +
1.141 + Fout->r = PSHR32(Fout->r, 2);
1.142 + Fout->i = PSHR32(Fout->i, 2);
1.143 + C_SUB( scratch[5] , *Fout, scratch[1] );
1.144 + C_ADDTO(*Fout, scratch[1]);
1.145 + C_ADD( scratch[3] , scratch[0] , scratch[2] );
1.146 + C_SUB( scratch[4] , scratch[0] , scratch[2] );
1.147 + C_SUB( Fout[m2], *Fout, scratch[3] );
1.148 + tw1 += fstride;
1.149 + tw2 += fstride*2;
1.150 + tw3 += fstride*3;
1.151 + C_ADDTO( *Fout , scratch[3] );
1.152 +
1.153 + Fout[m].r = scratch[5].r + scratch[4].i;
1.154 + Fout[m].i = scratch[5].i - scratch[4].r;
1.155 + Fout[m3].r = scratch[5].r - scratch[4].i;
1.156 + Fout[m3].i = scratch[5].i + scratch[4].r;
1.157 + ++Fout;
1.158 + }
1.159 + }
1.160 +}
1.161 +
1.162 +static void ki_bfly4(
1.163 + kiss_fft_cpx * Fout,
1.164 + const size_t fstride,
1.165 + const kiss_fft_state *st,
1.166 + int m,
1.167 + int N,
1.168 + int mm
1.169 + )
1.170 +{
1.171 + const kiss_twiddle_cpx *tw1,*tw2,*tw3;
1.172 + kiss_fft_cpx scratch[6];
1.173 + const size_t m2=2*m;
1.174 + const size_t m3=3*m;
1.175 + int i, j;
1.176 +
1.177 + kiss_fft_cpx * Fout_beg = Fout;
1.178 + for (i=0;i<N;i++)
1.179 + {
1.180 + Fout = Fout_beg + i*mm;
1.181 + tw3 = tw2 = tw1 = st->twiddles;
1.182 + for (j=0;j<m;j++)
1.183 + {
1.184 + C_MULC(scratch[0],Fout[m] , *tw1 );
1.185 + C_MULC(scratch[1],Fout[m2] , *tw2 );
1.186 + C_MULC(scratch[2],Fout[m3] , *tw3 );
1.187 +
1.188 + C_SUB( scratch[5] , *Fout, scratch[1] );
1.189 + C_ADDTO(*Fout, scratch[1]);
1.190 + C_ADD( scratch[3] , scratch[0] , scratch[2] );
1.191 + C_SUB( scratch[4] , scratch[0] , scratch[2] );
1.192 + C_SUB( Fout[m2], *Fout, scratch[3] );
1.193 + tw1 += fstride;
1.194 + tw2 += fstride*2;
1.195 + tw3 += fstride*3;
1.196 + C_ADDTO( *Fout , scratch[3] );
1.197 +
1.198 + Fout[m].r = scratch[5].r - scratch[4].i;
1.199 + Fout[m].i = scratch[5].i + scratch[4].r;
1.200 + Fout[m3].r = scratch[5].r + scratch[4].i;
1.201 + Fout[m3].i = scratch[5].i - scratch[4].r;
1.202 + ++Fout;
1.203 + }
1.204 + }
1.205 +}
1.206 +
1.207 +#ifndef RADIX_TWO_ONLY
1.208 +
1.209 +static void kf_bfly3(
1.210 + kiss_fft_cpx * Fout,
1.211 + const size_t fstride,
1.212 + const kiss_fft_state *st,
1.213 + int m,
1.214 + int N,
1.215 + int mm
1.216 + )
1.217 +{
1.218 + int i;
1.219 + size_t k;
1.220 + const size_t m2 = 2*m;
1.221 + const kiss_twiddle_cpx *tw1,*tw2;
1.222 + kiss_fft_cpx scratch[5];
1.223 + kiss_twiddle_cpx epi3;
1.224 +
1.225 + kiss_fft_cpx * Fout_beg = Fout;
1.226 + epi3 = st->twiddles[fstride*m];
1.227 + for (i=0;i<N;i++)
1.228 + {
1.229 + Fout = Fout_beg + i*mm;
1.230 + tw1=tw2=st->twiddles;
1.231 + k=m;
1.232 + do {
1.233 + C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
1.234 +
1.235 + C_MUL(scratch[1],Fout[m] , *tw1);
1.236 + C_MUL(scratch[2],Fout[m2] , *tw2);
1.237 +
1.238 + C_ADD(scratch[3],scratch[1],scratch[2]);
1.239 + C_SUB(scratch[0],scratch[1],scratch[2]);
1.240 + tw1 += fstride;
1.241 + tw2 += fstride*2;
1.242 +
1.243 + Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
1.244 + Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
1.245 +
1.246 + C_MULBYSCALAR( scratch[0] , epi3.i );
1.247 +
1.248 + C_ADDTO(*Fout,scratch[3]);
1.249 +
1.250 + Fout[m2].r = Fout[m].r + scratch[0].i;
1.251 + Fout[m2].i = Fout[m].i - scratch[0].r;
1.252 +
1.253 + Fout[m].r -= scratch[0].i;
1.254 + Fout[m].i += scratch[0].r;
1.255 +
1.256 + ++Fout;
1.257 + } while(--k);
1.258 + }
1.259 +}
1.260 +
1.261 +static void ki_bfly3(
1.262 + kiss_fft_cpx * Fout,
1.263 + const size_t fstride,
1.264 + const kiss_fft_state *st,
1.265 + int m,
1.266 + int N,
1.267 + int mm
1.268 + )
1.269 +{
1.270 + int i, k;
1.271 + const size_t m2 = 2*m;
1.272 + const kiss_twiddle_cpx *tw1,*tw2;
1.273 + kiss_fft_cpx scratch[5];
1.274 + kiss_twiddle_cpx epi3;
1.275 +
1.276 + kiss_fft_cpx * Fout_beg = Fout;
1.277 + epi3 = st->twiddles[fstride*m];
1.278 + for (i=0;i<N;i++)
1.279 + {
1.280 + Fout = Fout_beg + i*mm;
1.281 + tw1=tw2=st->twiddles;
1.282 + k=m;
1.283 + do{
1.284 +
1.285 + C_MULC(scratch[1],Fout[m] , *tw1);
1.286 + C_MULC(scratch[2],Fout[m2] , *tw2);
1.287 +
1.288 + C_ADD(scratch[3],scratch[1],scratch[2]);
1.289 + C_SUB(scratch[0],scratch[1],scratch[2]);
1.290 + tw1 += fstride;
1.291 + tw2 += fstride*2;
1.292 +
1.293 + Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
1.294 + Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
1.295 +
1.296 + C_MULBYSCALAR( scratch[0] , -epi3.i );
1.297 +
1.298 + C_ADDTO(*Fout,scratch[3]);
1.299 +
1.300 + Fout[m2].r = Fout[m].r + scratch[0].i;
1.301 + Fout[m2].i = Fout[m].i - scratch[0].r;
1.302 +
1.303 + Fout[m].r -= scratch[0].i;
1.304 + Fout[m].i += scratch[0].r;
1.305 +
1.306 + ++Fout;
1.307 + }while(--k);
1.308 + }
1.309 +}
1.310 +
1.311 +static void kf_bfly5(
1.312 + kiss_fft_cpx * Fout,
1.313 + const size_t fstride,
1.314 + const kiss_fft_state *st,
1.315 + int m,
1.316 + int N,
1.317 + int mm
1.318 + )
1.319 +{
1.320 + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
1.321 + int i, u;
1.322 + kiss_fft_cpx scratch[13];
1.323 + const kiss_twiddle_cpx * twiddles = st->twiddles;
1.324 + const kiss_twiddle_cpx *tw;
1.325 + kiss_twiddle_cpx ya,yb;
1.326 + kiss_fft_cpx * Fout_beg = Fout;
1.327 +
1.328 + ya = twiddles[fstride*m];
1.329 + yb = twiddles[fstride*2*m];
1.330 + tw=st->twiddles;
1.331 +
1.332 + for (i=0;i<N;i++)
1.333 + {
1.334 + Fout = Fout_beg + i*mm;
1.335 + Fout0=Fout;
1.336 + Fout1=Fout0+m;
1.337 + Fout2=Fout0+2*m;
1.338 + Fout3=Fout0+3*m;
1.339 + Fout4=Fout0+4*m;
1.340 +
1.341 + for ( u=0; u<m; ++u ) {
1.342 + C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
1.343 + scratch[0] = *Fout0;
1.344 +
1.345 + C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
1.346 + C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
1.347 + C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
1.348 + C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
1.349 +
1.350 + C_ADD( scratch[7],scratch[1],scratch[4]);
1.351 + C_SUB( scratch[10],scratch[1],scratch[4]);
1.352 + C_ADD( scratch[8],scratch[2],scratch[3]);
1.353 + C_SUB( scratch[9],scratch[2],scratch[3]);
1.354 +
1.355 + Fout0->r += scratch[7].r + scratch[8].r;
1.356 + Fout0->i += scratch[7].i + scratch[8].i;
1.357 +
1.358 + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
1.359 + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
1.360 +
1.361 + scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
1.362 + scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
1.363 +
1.364 + C_SUB(*Fout1,scratch[5],scratch[6]);
1.365 + C_ADD(*Fout4,scratch[5],scratch[6]);
1.366 +
1.367 + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
1.368 + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
1.369 + scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
1.370 + scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
1.371 +
1.372 + C_ADD(*Fout2,scratch[11],scratch[12]);
1.373 + C_SUB(*Fout3,scratch[11],scratch[12]);
1.374 +
1.375 + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
1.376 + }
1.377 + }
1.378 +}
1.379 +
1.380 +static void ki_bfly5(
1.381 + kiss_fft_cpx * Fout,
1.382 + const size_t fstride,
1.383 + const kiss_fft_state *st,
1.384 + int m,
1.385 + int N,
1.386 + int mm
1.387 + )
1.388 +{
1.389 + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
1.390 + int i, u;
1.391 + kiss_fft_cpx scratch[13];
1.392 + const kiss_twiddle_cpx * twiddles = st->twiddles;
1.393 + const kiss_twiddle_cpx *tw;
1.394 + kiss_twiddle_cpx ya,yb;
1.395 + kiss_fft_cpx * Fout_beg = Fout;
1.396 +
1.397 + ya = twiddles[fstride*m];
1.398 + yb = twiddles[fstride*2*m];
1.399 + tw=st->twiddles;
1.400 +
1.401 + for (i=0;i<N;i++)
1.402 + {
1.403 + Fout = Fout_beg + i*mm;
1.404 + Fout0=Fout;
1.405 + Fout1=Fout0+m;
1.406 + Fout2=Fout0+2*m;
1.407 + Fout3=Fout0+3*m;
1.408 + Fout4=Fout0+4*m;
1.409 +
1.410 + for ( u=0; u<m; ++u ) {
1.411 + scratch[0] = *Fout0;
1.412 +
1.413 + C_MULC(scratch[1] ,*Fout1, tw[u*fstride]);
1.414 + C_MULC(scratch[2] ,*Fout2, tw[2*u*fstride]);
1.415 + C_MULC(scratch[3] ,*Fout3, tw[3*u*fstride]);
1.416 + C_MULC(scratch[4] ,*Fout4, tw[4*u*fstride]);
1.417 +
1.418 + C_ADD( scratch[7],scratch[1],scratch[4]);
1.419 + C_SUB( scratch[10],scratch[1],scratch[4]);
1.420 + C_ADD( scratch[8],scratch[2],scratch[3]);
1.421 + C_SUB( scratch[9],scratch[2],scratch[3]);
1.422 +
1.423 + Fout0->r += scratch[7].r + scratch[8].r;
1.424 + Fout0->i += scratch[7].i + scratch[8].i;
1.425 +
1.426 + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
1.427 + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
1.428 +
1.429 + scratch[6].r = -S_MUL(scratch[10].i,ya.i) - S_MUL(scratch[9].i,yb.i);
1.430 + scratch[6].i = S_MUL(scratch[10].r,ya.i) + S_MUL(scratch[9].r,yb.i);
1.431 +
1.432 + C_SUB(*Fout1,scratch[5],scratch[6]);
1.433 + C_ADD(*Fout4,scratch[5],scratch[6]);
1.434 +
1.435 + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
1.436 + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
1.437 + scratch[12].r = S_MUL(scratch[10].i,yb.i) - S_MUL(scratch[9].i,ya.i);
1.438 + scratch[12].i = -S_MUL(scratch[10].r,yb.i) + S_MUL(scratch[9].r,ya.i);
1.439 +
1.440 + C_ADD(*Fout2,scratch[11],scratch[12]);
1.441 + C_SUB(*Fout3,scratch[11],scratch[12]);
1.442 +
1.443 + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
1.444 + }
1.445 + }
1.446 +}
1.447 +
1.448 +#endif
1.449 +
1.450 +
1.451 +#ifdef CUSTOM_MODES
1.452 +
1.453 +static
1.454 +void compute_bitrev_table(
1.455 + int Fout,
1.456 + opus_int16 *f,
1.457 + const size_t fstride,
1.458 + int in_stride,
1.459 + opus_int16 * factors,
1.460 + const kiss_fft_state *st
1.461 + )
1.462 +{
1.463 + const int p=*factors++; /* the radix */
1.464 + const int m=*factors++; /* stage's fft length/p */
1.465 +
1.466 + /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/
1.467 + if (m==1)
1.468 + {
1.469 + int j;
1.470 + for (j=0;j<p;j++)
1.471 + {
1.472 + *f = Fout+j;
1.473 + f += fstride*in_stride;
1.474 + }
1.475 + } else {
1.476 + int j;
1.477 + for (j=0;j<p;j++)
1.478 + {
1.479 + compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st);
1.480 + f += fstride*in_stride;
1.481 + Fout += m;
1.482 + }
1.483 + }
1.484 +}
1.485 +
1.486 +/* facbuf is populated by p1,m1,p2,m2, ...
1.487 + where
1.488 + p[i] * m[i] = m[i-1]
1.489 + m0 = n */
1.490 +static
1.491 +int kf_factor(int n,opus_int16 * facbuf)
1.492 +{
1.493 + int p=4;
1.494 +
1.495 + /*factor out powers of 4, powers of 2, then any remaining primes */
1.496 + do {
1.497 + while (n % p) {
1.498 + switch (p) {
1.499 + case 4: p = 2; break;
1.500 + case 2: p = 3; break;
1.501 + default: p += 2; break;
1.502 + }
1.503 + if (p>32000 || (opus_int32)p*(opus_int32)p > n)
1.504 + p = n; /* no more factors, skip to end */
1.505 + }
1.506 + n /= p;
1.507 +#ifdef RADIX_TWO_ONLY
1.508 + if (p!=2 && p != 4)
1.509 +#else
1.510 + if (p>5)
1.511 +#endif
1.512 + {
1.513 + return 0;
1.514 + }
1.515 + *facbuf++ = p;
1.516 + *facbuf++ = n;
1.517 + } while (n > 1);
1.518 + return 1;
1.519 +}
1.520 +
1.521 +static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft)
1.522 +{
1.523 + int i;
1.524 +#ifdef FIXED_POINT
1.525 + for (i=0;i<nfft;++i) {
1.526 + opus_val32 phase = -i;
1.527 + kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft));
1.528 + }
1.529 +#else
1.530 + for (i=0;i<nfft;++i) {
1.531 + const double pi=3.14159265358979323846264338327;
1.532 + double phase = ( -2*pi /nfft ) * i;
1.533 + kf_cexp(twiddles+i, phase );
1.534 + }
1.535 +#endif
1.536 +}
1.537 +
1.538 +/*
1.539 + *
1.540 + * Allocates all necessary storage space for the fft and ifft.
1.541 + * The return value is a contiguous block of memory. As such,
1.542 + * It can be freed with free().
1.543 + * */
1.544 +kiss_fft_state *opus_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem, const kiss_fft_state *base)
1.545 +{
1.546 + kiss_fft_state *st=NULL;
1.547 + size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/
1.548 +
1.549 + if ( lenmem==NULL ) {
1.550 + st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded );
1.551 + }else{
1.552 + if (mem != NULL && *lenmem >= memneeded)
1.553 + st = (kiss_fft_state*)mem;
1.554 + *lenmem = memneeded;
1.555 + }
1.556 + if (st) {
1.557 + opus_int16 *bitrev;
1.558 + kiss_twiddle_cpx *twiddles;
1.559 +
1.560 + st->nfft=nfft;
1.561 +#ifndef FIXED_POINT
1.562 + st->scale = 1.f/nfft;
1.563 +#endif
1.564 + if (base != NULL)
1.565 + {
1.566 + st->twiddles = base->twiddles;
1.567 + st->shift = 0;
1.568 + while (nfft<<st->shift != base->nfft && st->shift < 32)
1.569 + st->shift++;
1.570 + if (st->shift>=32)
1.571 + goto fail;
1.572 + } else {
1.573 + st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft);
1.574 + compute_twiddles(twiddles, nfft);
1.575 + st->shift = -1;
1.576 + }
1.577 + if (!kf_factor(nfft,st->factors))
1.578 + {
1.579 + goto fail;
1.580 + }
1.581 +
1.582 + /* bitrev */
1.583 + st->bitrev = bitrev = (opus_int16*)KISS_FFT_MALLOC(sizeof(opus_int16)*nfft);
1.584 + if (st->bitrev==NULL)
1.585 + goto fail;
1.586 + compute_bitrev_table(0, bitrev, 1,1, st->factors,st);
1.587 + }
1.588 + return st;
1.589 +fail:
1.590 + opus_fft_free(st);
1.591 + return NULL;
1.592 +}
1.593 +
1.594 +kiss_fft_state *opus_fft_alloc(int nfft,void * mem,size_t * lenmem )
1.595 +{
1.596 + return opus_fft_alloc_twiddles(nfft, mem, lenmem, NULL);
1.597 +}
1.598 +
1.599 +void opus_fft_free(const kiss_fft_state *cfg)
1.600 +{
1.601 + if (cfg)
1.602 + {
1.603 + opus_free((opus_int16*)cfg->bitrev);
1.604 + if (cfg->shift < 0)
1.605 + opus_free((kiss_twiddle_cpx*)cfg->twiddles);
1.606 + opus_free((kiss_fft_state*)cfg);
1.607 + }
1.608 +}
1.609 +
1.610 +#endif /* CUSTOM_MODES */
1.611 +
1.612 +void opus_fft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
1.613 +{
1.614 + int m2, m;
1.615 + int p;
1.616 + int L;
1.617 + int fstride[MAXFACTORS];
1.618 + int i;
1.619 + int shift;
1.620 +
1.621 + /* st->shift can be -1 */
1.622 + shift = st->shift>0 ? st->shift : 0;
1.623 +
1.624 + celt_assert2 (fin != fout, "In-place FFT not supported");
1.625 + /* Bit-reverse the input */
1.626 + for (i=0;i<st->nfft;i++)
1.627 + {
1.628 + fout[st->bitrev[i]] = fin[i];
1.629 +#ifndef FIXED_POINT
1.630 + fout[st->bitrev[i]].r *= st->scale;
1.631 + fout[st->bitrev[i]].i *= st->scale;
1.632 +#endif
1.633 + }
1.634 +
1.635 + fstride[0] = 1;
1.636 + L=0;
1.637 + do {
1.638 + p = st->factors[2*L];
1.639 + m = st->factors[2*L+1];
1.640 + fstride[L+1] = fstride[L]*p;
1.641 + L++;
1.642 + } while(m!=1);
1.643 + m = st->factors[2*L-1];
1.644 + for (i=L-1;i>=0;i--)
1.645 + {
1.646 + if (i!=0)
1.647 + m2 = st->factors[2*i-1];
1.648 + else
1.649 + m2 = 1;
1.650 + switch (st->factors[2*i])
1.651 + {
1.652 + case 2:
1.653 + kf_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.654 + break;
1.655 + case 4:
1.656 + kf_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.657 + break;
1.658 + #ifndef RADIX_TWO_ONLY
1.659 + case 3:
1.660 + kf_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.661 + break;
1.662 + case 5:
1.663 + kf_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.664 + break;
1.665 + #endif
1.666 + }
1.667 + m = m2;
1.668 + }
1.669 +}
1.670 +
1.671 +void opus_ifft(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
1.672 +{
1.673 + int m2, m;
1.674 + int p;
1.675 + int L;
1.676 + int fstride[MAXFACTORS];
1.677 + int i;
1.678 + int shift;
1.679 +
1.680 + /* st->shift can be -1 */
1.681 + shift = st->shift>0 ? st->shift : 0;
1.682 + celt_assert2 (fin != fout, "In-place FFT not supported");
1.683 + /* Bit-reverse the input */
1.684 + for (i=0;i<st->nfft;i++)
1.685 + fout[st->bitrev[i]] = fin[i];
1.686 +
1.687 + fstride[0] = 1;
1.688 + L=0;
1.689 + do {
1.690 + p = st->factors[2*L];
1.691 + m = st->factors[2*L+1];
1.692 + fstride[L+1] = fstride[L]*p;
1.693 + L++;
1.694 + } while(m!=1);
1.695 + m = st->factors[2*L-1];
1.696 + for (i=L-1;i>=0;i--)
1.697 + {
1.698 + if (i!=0)
1.699 + m2 = st->factors[2*i-1];
1.700 + else
1.701 + m2 = 1;
1.702 + switch (st->factors[2*i])
1.703 + {
1.704 + case 2:
1.705 + ki_bfly2(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.706 + break;
1.707 + case 4:
1.708 + ki_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.709 + break;
1.710 +#ifndef RADIX_TWO_ONLY
1.711 + case 3:
1.712 + ki_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.713 + break;
1.714 + case 5:
1.715 + ki_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2);
1.716 + break;
1.717 +#endif
1.718 + }
1.719 + m = m2;
1.720 + }
1.721 +}
1.722 +
