fixed-point: converted intra prediction and folding, unb0rked mixed-precision
[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 #define EXT32(a) (a)
65
66 #else /* DOUBLE_PRECISION */
67
68 # define FRACBITS 15
69 # define SAMPPROD celt_int32_t 
70 #define SAMP_MAX 32767
71 #define TRIG_UPSCALE 1
72 #define EXT32(a) EXTEND32(a)
73
74 #endif /* !DOUBLE_PRECISION */
75
76 #define SAMP_MIN -SAMP_MAX
77
78 #if defined(CHECK_OVERFLOW)
79 #  define CHECK_OVERFLOW_OP(a,op,b)  \
80         if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
81                 fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) );  }
82 #endif
83
84 #   define smul(a,b) ( (SAMPPROD)(a)*(b) )
85 #   define sround( x )  (kiss_fft_scalar)( ( (x) + ((SAMPPROD)1<<(FRACBITS-1)) ) >> FRACBITS )
86
87 #ifdef MIXED_PRECISION
88
89 #undef MULT16_32_Q15
90 #define MULT16_16SU(a,b) ((celt_word32_t)(celt_word16_t)(a)*(celt_word32_t)(celt_uint16_t)(b))
91 /*#define MULT16_32_Q15(a,b) ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15))*/
92 #define MULT16_32_Q15(a,b) ADD32(SHL(MULT16_16((a),SHR((b),16)),1), SHR(MULT16_16SU((a),((b)&0x0000ffff)),15))
93
94 #   define S_MUL(a,b) MULT16_32_Q15(b, a)
95
96 #   define C_MUL(m,a,b) \
97       do{ (m).r = S_MUL((a).r,(b).r) - S_MUL((a).i,(b).i); \
98           (m).i = S_MUL((a).r,(b).i) + S_MUL((a).i,(b).r); }while(0)
99
100 #   define C_MULC(m,a,b) \
101       do{ (m).r = S_MUL((a).r,(b).r) + S_MUL((a).i,(b).i); \
102           (m).i = S_MUL((a).i,(b).r) - S_MUL((a).r,(b).i); }while(0)
103
104 #   define C_MUL4(m,a,b) \
105       do{ (m).r = SHR(S_MUL((a).r,(b).r) - S_MUL((a).i,(b).i),2); \
106           (m).i = SHR(S_MUL((a).r,(b).i) + S_MUL((a).i,(b).r),2); }while(0)
107
108 #   define C_MULBYSCALAR( c, s ) \
109       do{ (c).r =  S_MUL( (c).r , s ) ;\
110           (c).i =  S_MUL( (c).i , s ) ; }while(0)
111
112 #   define DIVSCALAR(x,k) \
113         (x) = S_MUL(  x, (TWID_MAX-((k)>>1))/(k)+1 )
114
115 #   define C_FIXDIV(c,div) \
116         do {    DIVSCALAR( (c).r , div);  \
117                 DIVSCALAR( (c).i  , div); }while (0)
118
119 #else /* MIXED_PRECISION */
120 #   define sround4( x )  (kiss_fft_scalar)( ( (x) + ((SAMPPROD)1<<(FRACBITS-1)) ) >> (FRACBITS+2) )
121
122 #   define S_MUL(a,b) sround( smul(a,b) )
123
124 #   define C_MUL(m,a,b) \
125       do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
126           (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
127 #   define C_MULC(m,a,b) \
128       do{ (m).r = sround( smul((a).r,(b).r) + smul((a).i,(b).i) ); \
129           (m).i = sround( smul((a).i,(b).r) - smul((a).r,(b).i) ); }while(0)
130
131 #   define C_MUL4(m,a,b) \
132                do{ (m).r = sround4( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
133                (m).i = sround4( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
134
135 #   define C_MULBYSCALAR( c, s ) \
136                do{ (c).r =  sround( smul( (c).r , s ) ) ;\
137                (c).i =  sround( smul( (c).i , s ) ) ; }while(0)
138
139 #   define DIVSCALAR(x,k) \
140         (x) = sround( smul(  x, SAMP_MAX/k ) )
141
142 #   define C_FIXDIV(c,div) \
143         do {    DIVSCALAR( (c).r , div);  \
144                 DIVSCALAR( (c).i  , div); }while (0)
145
146 #endif /* !MIXED_PRECISION */
147
148
149
150 #else  /* not FIXED_POINT*/
151
152 #define EXT32(a) (a)
153
154 #   define S_MUL(a,b) ( (a)*(b) )
155 #define C_MUL(m,a,b) \
156     do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
157         (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
158 #define C_MULC(m,a,b) \
159     do{ (m).r = (a).r*(b).r + (a).i*(b).i;\
160         (m).i = (a).i*(b).r - (a).r*(b).i; }while(0)
161
162 #define C_MUL4(m,a,b) C_MUL(m,a,b)
163
164 #   define C_FIXDIV(c,div) /* NOOP */
165 #   define C_MULBYSCALAR( c, s ) \
166     do{ (c).r *= (s);\
167         (c).i *= (s); }while(0)
168 #endif
169
170 #ifndef CHECK_OVERFLOW_OP
171 #  define CHECK_OVERFLOW_OP(a,op,b) /* noop */
172 #endif
173
174 #define  C_ADD( res, a,b)\
175     do { \
176             CHECK_OVERFLOW_OP((a).r,+,(b).r)\
177             CHECK_OVERFLOW_OP((a).i,+,(b).i)\
178             (res).r=(a).r+(b).r;  (res).i=(a).i+(b).i; \
179     }while(0)
180 #define  C_SUB( res, a,b)\
181     do { \
182             CHECK_OVERFLOW_OP((a).r,-,(b).r)\
183             CHECK_OVERFLOW_OP((a).i,-,(b).i)\
184             (res).r=(a).r-(b).r;  (res).i=(a).i-(b).i; \
185     }while(0)
186 #define C_ADDTO( res , a)\
187     do { \
188             CHECK_OVERFLOW_OP((res).r,+,(a).r)\
189             CHECK_OVERFLOW_OP((res).i,+,(a).i)\
190             (res).r += (a).r;  (res).i += (a).i;\
191     }while(0)
192
193 #define C_SUBFROM( res , a)\
194     do {\
195             CHECK_OVERFLOW_OP((res).r,-,(a).r)\
196             CHECK_OVERFLOW_OP((res).i,-,(a).i)\
197             (res).r -= (a).r;  (res).i -= (a).i; \
198     }while(0)
199
200
201 #ifdef FIXED_POINT
202 /*#  define KISS_FFT_COS(phase)  TRIG_UPSCALE*floor(MIN(32767,MAX(-32767,.5+32768 * cos (phase))))
203 #  define KISS_FFT_SIN(phase)  TRIG_UPSCALE*floor(MIN(32767,MAX(-32767,.5+32768 * sin (phase))))*/
204 #  define KISS_FFT_COS(phase)  floor(.5+TWID_MAX*cos (phase))
205 #  define KISS_FFT_SIN(phase)  floor(.5+TWID_MAX*sin (phase))
206 #  define HALF_OF(x) ((x)>>1)
207 #elif defined(USE_SIMD)
208 #  define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
209 #  define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
210 #  define HALF_OF(x) ((x)*_mm_set1_ps(.5))
211 #else
212 #  define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
213 #  define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
214 #  define HALF_OF(x) ((x)*.5)
215 #endif
216
217 #define  kf_cexp(x,phase) \
218         do{ \
219                 (x)->r = KISS_FFT_COS(phase);\
220                 (x)->i = KISS_FFT_SIN(phase);\
221         }while(0)
222    
223 #define  kf_cexp2(x,phase) \
224    do{ \
225       (x)->r = TRIG_UPSCALE*celt_cos_norm((phase));\
226       (x)->i = TRIG_UPSCALE*celt_cos_norm((phase)-32768);\
227 }while(0)
228
229 /* a debugging function */
230 #define pcpx(c)\
231     fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )