apps/plugins/fft/_kiss_fft_guts.h

   1 /*
   2 Copyright (c) 2003-2004, Mark Borgerding
   3
   4 All rights reserved.
   5
   6 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
   7
   8     * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
   9     * 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.
  10     * 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.
  11
  12 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.
  13 */
  14
  15 /* kiss_fft.h
  16    defines kiss_fft_scalar as either short or a float type
  17    and defines
  18    typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
  19 #include "kiss_fft.h"
  20 #include "math.h"
  21 #include <limits.h>
  22
  23 #define MAXFACTORS 32
  24 /* e.g. an fft of length 128 has 4 factors
  25  as far as kissfft is concerned
  26  4*4*4*2
  27  */
  28
  29 struct kiss_fft_state{
  30     int nfft;
  31     int inverse;
  32     int factors[2*MAXFACTORS];
  33     kiss_fft_cpx twiddles[1];
  34 };
  35
  36 /*
  37   Explanation of macros dealing with complex math:
  38
  39    C_MUL(m,a,b)         : m = a*b
  40    C_FIXDIV( c , div )  : if a fixed point impl., c /= div. noop otherwise
  41    C_SUB( res, a,b)     : res = a - b
  42    C_SUBFROM( res , a)  : res -= a
  43    C_ADDTO( res , a)    : res += a
  44  * */
  45 #ifdef FIXED_POINT
  46 #if (FIXED_POINT==32)
  47 # define FRACBITS 31
  48 # define SAMPPROD int64_t
  49 #define SAMP_MAX 2147483647
  50 #else
  51 # define FRACBITS 15
  52 # define SAMPPROD int32_t
  53 #define SAMP_MAX 32767
  54 #endif
  55
  56 #define SAMP_MIN -SAMP_MAX
  57
  58 #if defined(CHECK_OVERFLOW)
  59 #  define CHECK_OVERFLOW_OP(a,op,b)  \
  60         if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
  61                 fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) );  }
  62 #endif
  63
  64
  65 #   define smul(a,b) ( (SAMPPROD)(a)*(b) )
  66 #   define sround( x )  (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )
  67
  68 #   define S_MUL(a,b) sround( smul(a,b) )
  69
  70 #   define C_MUL(m,a,b) \
  71       do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
  72           (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
  73
  74 #   define DIVSCALAR(x,k) \
  75         (x) = sround( smul(  x, SAMP_MAX/k ) )
  76
  77 #   define C_FIXDIV(c,div) \
  78         do {    DIVSCALAR( (c).r , div);  \
  79                 DIVSCALAR( (c).i  , div); }while (0)
  80
  81 #   define C_MULBYSCALAR( c, s ) \
  82     do{ (c).r =  sround( smul( (c).r , s ) ) ;\
  83         (c).i =  sround( smul( (c).i , s ) ) ; }while(0)
  84
  85 #else  /* not FIXED_POINT*/
  86
  87 #   define S_MUL(a,b) ( (a)*(b) )
  88 #define C_MUL(m,a,b) \
  89     do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
  90         (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
  91 #   define C_FIXDIV(c,div) /* NOOP */
  92 #   define C_MULBYSCALAR( c, s ) \
  93     do{ (c).r *= (s);\
  94         (c).i *= (s); }while(0)
  95 #endif
  96
  97 #ifndef CHECK_OVERFLOW_OP
  98 #  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
  99 #endif
 100
 101 #define  C_ADD( res, a,b)\
 102     do { \
 103             CHECK_OVERFLOW_OP((a).r,+,(b).r)\
 104             CHECK_OVERFLOW_OP((a).i,+,(b).i)\
 105             (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
 106     }while(0)
 107 #define  C_SUB( res, a,b)\
 108     do { \
 109             CHECK_OVERFLOW_OP((a).r,-,(b).r)\
 110             CHECK_OVERFLOW_OP((a).i,-,(b).i)\
 111             (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
 112     }while(0)
 113 #define C_ADDTO( res , a)\
 114     do { \
 115             CHECK_OVERFLOW_OP((res).r,+,(a).r)\
 116             CHECK_OVERFLOW_OP((res).i,+,(a).i)\
 117             (res).r += (a).r;  (res).i += (a).i;\
 118     }while(0)
 119
 120 #define C_SUBFROM( res , a)\
 121     do {\
 122             CHECK_OVERFLOW_OP((res).r,-,(a).r)\
 123             CHECK_OVERFLOW_OP((res).i,-,(a).i)\
 124             (res).r -= (a).r;  (res).i -= (a).i; \
 125     }while(0)
 126
 127
 128 #ifdef FIXED_POINT
 129 #  define HALF_OF(x) ((x)>>1)
 130 #else
 131 #  define HALF_OF(x) ((x)*.5)
 132 #endif
 133
 134 #define  kf_cexp(x, k, n) \
 135         do{ \
 136                 int32_t div = Q_DIV( (k) << 16, (n) << 16, 16 ); \
 137                 long cos, sin = fp_sincos(div << 16, &cos); \
 138                 (x)->r = ( Q_MUL(SAMP_MAX << 16, cos >> 15, 16) ) >> 16; \
 139                 (x)->i = ( Q_MUL(SAMP_MAX << 16, -1*(sin >> 15), 16) ) >> 16; \
 140         }while(0)
 141
 142 #define  kf_cexp_round(x, k, n) \
 143         do{ \
 144                 int32_t div = Q_DIV( (k) << 16, (n) << 16, 16 ) + (1 << 15); \
 145                 long cos, sin = fp_sincos(div << 16, &cos); \
 146                 (x)->r = ( Q_MUL(SAMP_MAX << 16, cos >> 15, 16) ) >> 16; \
 147                 (x)->i = ( Q_MUL(SAMP_MAX << 16, -1*(sin >> 15), 16) ) >> 16; \
 148         }while(0)
 149
 150 /* a debugging function */
 151 #define pcpx(c)\
 152     fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )