fixed-point: initial support for using the fixed-point MDCT (rest is still all
[opus.git] / libcelt / _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 #define MIN(a,b) ((a)<(b) ? (a):(b))
16 #define MAX(a,b) ((a)>(b) ? (a):(b))
17
18 /* kiss_fft.h
19    defines kiss_fft_scalar as either short or a float type
20    and defines
21    typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
22 #include "kiss_fft.h"
23
24 #define MAXFACTORS 32
25 /* e.g. an fft of length 128 has 4 factors 
26  as far as kissfft is concerned
27  4*4*4*2
28  */
29
30 struct kiss_fft_state{
31     int nfft;
32 #ifndef FIXED_POINT
33     kiss_fft_scalar scale;
34 #endif
35     int factors[2*MAXFACTORS];
36     int *bitrev;
37     kiss_twiddle_cpx twiddles[1];
38 };
39
40 /*
41   Explanation of macros dealing with complex math:
42
43    C_MUL(m,a,b)         : m = a*b
44    C_FIXDIV( c , div )  : if a fixed point impl., c /= div. noop otherwise
45    C_SUB( res, a,b)     : res = a - b
46    C_SUBFROM( res , a)  : res -= a
47    C_ADDTO( res , a)    : res += a
48  * */
49 #ifdef FIXED_POINT
50 #include "arch.h"
51
52 #ifdef DOUBLE_PRECISION
53
54 # define FRACBITS 31
55 # define SAMPPROD celt_int64_t 
56 #define SAMP_MAX 2147483647
57 #ifdef MIXED_PRECISION
58 #define TWID_MAX 32767
59 #define TRIG_UPSCALE 1
60 #else
61 #define TRIG_UPSCALE 65536
62 #define TWID_MAX 2147483647
63 #endif
64 #else /* DOUBLE_PRECISION */
65
66 # define FRACBITS 15
67 # define SAMPPROD celt_int32_t 
68 #define SAMP_MAX 32767
69 #define TRIG_UPSCALE 1
70
71 #endif /* !DOUBLE_PRECISION */
72
73 #define SAMP_MIN -SAMP_MAX
74
75 #if defined(CHECK_OVERFLOW)
76 #  define CHECK_OVERFLOW_OP(a,op,b)  \
77         if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
78                 fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) );  }
79 #endif
80
81 #   define smul(a,b) ( (SAMPPROD)(a)*(b) )
82 #   define sround( x )  (kiss_fft_scalar)( ( (x) + ((SAMPPROD)1<<(FRACBITS-1)) ) >> FRACBITS )
83
84 #ifdef MIXED_PRECISION
85
86 #undef MULT16_32_Q15
87 #define MULT16_16SU(a,b) ((celt_word32_t)(celt_word16_t)(a)*(celt_word32_t)(celt_uint16_t)(b))
88 /*#define MULT16_32_Q15(a,b) ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15))*/
89 #define MULT16_32_Q15(a,b) ADD32(SHL(MULT16_16((a),SHR((b),16)),1), SHR(MULT16_16SU((a),((b)&0x0000ffff)),15))
90
91 #   define S_MUL(a,b) MULT16_32_Q15(b, a)
92
93 #   define C_MUL(m,a,b) \
94       do{ (m).r = S_MUL((a).r,(b).r) - S_MUL((a).i,(b).i); \
95           (m).i = S_MUL((a).r,(b).i) + S_MUL((a).i,(b).r); }while(0)
96
97 #   define C_MULC(m,a,b) \
98       do{ (m).r = S_MUL((a).r,(b).r) + S_MUL((a).i,(b).i); \
99           (m).i = S_MUL((a).i,(b).r) - S_MUL((a).r,(b).i); }while(0)
100
101 #   define C_MUL4(m,a,b) \
102       do{ (m).r = SHR(S_MUL((a).r,(b).r) - S_MUL((a).i,(b).i),2); \
103           (m).i = SHR(S_MUL((a).r,(b).i) + S_MUL((a).i,(b).r),2); }while(0)
104
105 #   define C_MULBYSCALAR( c, s ) \
106       do{ (c).r =  S_MUL( (c).r , s ) ;\
107           (c).i =  S_MUL( (c).i , s ) ; }while(0)
108
109 #   define DIVSCALAR(x,k) \
110         (x) = S_MUL(  x, (TWID_MAX-((k)>>1))/(k)+1 )
111
112 #   define C_FIXDIV(c,div) \
113         do {    DIVSCALAR( (c).r , div);  \
114                 DIVSCALAR( (c).i  , div); }while (0)
115
116 #else /* MIXED_PRECISION */
117 #   define sround4( x )  (kiss_fft_scalar)( ( (x) + ((SAMPPROD)1<<(FRACBITS-1)) ) >> (FRACBITS+2) )
118
119 #   define S_MUL(a,b) sround( smul(a,b) )
120
121 #   define C_MUL(m,a,b) \
122       do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
123           (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
124 #   define C_MULC(m,a,b) \
125       do{ (m).r = sround( smul((a).r,(b).r) + smul((a).i,(b).i) ); \
126           (m).i = sround( smul((a).i,(b).r) - smul((a).r,(b).i) ); }while(0)
127
128 #   define C_MUL4(m,a,b) \
129                do{ (m).r = sround4( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
130                (m).i = sround4( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
131
132 #   define C_MULBYSCALAR( c, s ) \
133                do{ (c).r =  sround( smul( (c).r , s ) ) ;\
134                (c).i =  sround( smul( (c).i , s ) ) ; }while(0)
135
136 #   define DIVSCALAR(x,k) \
137         (x) = sround( smul(  x, SAMP_MAX/k ) )
138
139 #   define C_FIXDIV(c,div) \
140         do {    DIVSCALAR( (c).r , div);  \
141                 DIVSCALAR( (c).i  , div); }while (0)
142
143 #endif /* !MIXED_PRECISION */
144
145
146                
147 #define L1 32767
148 #define L2 -7651
149 #define L3 8277
150 #define L4 -626
151
152 static inline celt_word16_t _celt_cos_pi_2(celt_word16_t x)
153 {
154    celt_word16_t x2;
155    
156    x2 = MULT16_16_P15(x,x);
157    return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2
158                                                                                 ))))))));
159 }
160
161 static inline celt_word16_t celt_cos_norm(celt_word32_t x)
162 {
163    x = x&0x0001ffff;
164    if (x>SHL32(EXTEND32(1), 16))
165       x = SUB32(SHL32(EXTEND32(1), 17),x);
166    if (x&0x00007fff)
167    {
168       if (x<SHL32(EXTEND32(1), 15))
169       {
170          return _celt_cos_pi_2(EXTRACT16(x));
171       } else {
172          return NEG32(_celt_cos_pi_2(EXTRACT16(65536-x)));
173       }
174    } else {
175       if (x&0x0000ffff)
176          return 0;
177       else if (x&0x0001ffff)
178          return -32767;
179       else
180          return 32767;
181    }
182 }
183
184 #else  /* not FIXED_POINT*/
185
186 #   define S_MUL(a,b) ( (a)*(b) )
187 #define C_MUL(m,a,b) \
188     do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
189         (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
190 #define C_MULC(m,a,b) \
191     do{ (m).r = (a).r*(b).r + (a).i*(b).i;\
192         (m).i = (a).i*(b).r - (a).r*(b).i; }while(0)
193
194 #define C_MUL4(m,a,b) C_MUL(m,a,b)
195
196 #   define C_FIXDIV(c,div) /* NOOP */
197 #   define C_MULBYSCALAR( c, s ) \
198     do{ (c).r *= (s);\
199         (c).i *= (s); }while(0)
200 #endif
201
202 #ifndef CHECK_OVERFLOW_OP
203 #  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
204 #endif
205
206 #define  C_ADD( res, a,b)\
207     do { \
208             CHECK_OVERFLOW_OP((a).r,+,(b).r)\
209             CHECK_OVERFLOW_OP((a).i,+,(b).i)\
210             (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
211     }while(0)
212 #define  C_SUB( res, a,b)\
213     do { \
214             CHECK_OVERFLOW_OP((a).r,-,(b).r)\
215             CHECK_OVERFLOW_OP((a).i,-,(b).i)\
216             (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
217     }while(0)
218 #define C_ADDTO( res , a)\
219     do { \
220             CHECK_OVERFLOW_OP((res).r,+,(a).r)\
221             CHECK_OVERFLOW_OP((res).i,+,(a).i)\
222             (res).r += (a).r;  (res).i += (a).i;\
223     }while(0)
224
225 #define C_SUBFROM( res , a)\
226     do {\
227             CHECK_OVERFLOW_OP((res).r,-,(a).r)\
228             CHECK_OVERFLOW_OP((res).i,-,(a).i)\
229             (res).r -= (a).r;  (res).i -= (a).i; \
230     }while(0)
231
232
233 #ifdef FIXED_POINT
234 /*#  define KISS_FFT_COS(phase)  TRIG_UPSCALE*floor(MIN(32767,MAX(-32767,.5+32768 * cos (phase))))
235 #  define KISS_FFT_SIN(phase)  TRIG_UPSCALE*floor(MIN(32767,MAX(-32767,.5+32768 * sin (phase))))*/
236 #  define KISS_FFT_COS(phase)  floor(.5+TWID_MAX*cos (phase))
237 #  define KISS_FFT_SIN(phase)  floor(.5+TWID_MAX*sin (phase))
238 #  define HALF_OF(x) ((x)>>1)
239 #elif defined(USE_SIMD)
240 #  define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
241 #  define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
242 #  define HALF_OF(x) ((x)*_mm_set1_ps(.5))
243 #else
244 #  define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
245 #  define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
246 #  define HALF_OF(x) ((x)*.5)
247 #endif
248
249 #define  kf_cexp(x,phase) \
250         do{ \
251                 (x)->r = KISS_FFT_COS(phase);\
252                 (x)->i = KISS_FFT_SIN(phase);\
253         }while(0)
254    
255 #define  kf_cexp2(x,phase) \
256    do{ \
257       (x)->r = TRIG_UPSCALE*celt_cos_norm((phase));\
258       (x)->i = TRIG_UPSCALE*celt_cos_norm((phase)-32768);\
259 }while(0)
260
261 /* a debugging function */
262 #define pcpx(c)\
263     fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )