1 /*
2  *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3  *
4  *  Use of this source code is governed by a BSD-style license
5  *  that can be found in the LICENSE file in the root of the source
6  *  tree. An additional intellectual property rights grant can be found
7  *  in the file PATENTS.  All contributing project authors may
8  *  be found in the AUTHORS file in the root of the source tree.
9  */
10 
11 /*
12  * arith_routines.h
13  *
14  * This file contains functions for arithmatically encoding and
15  * decoding DFT coefficients.
16  *
17  */
18 
19 
20 #include "arith_routines.h"
21 
22 
23 
24 static const int32_t kHistEdgesQ15[51] = {
25   -327680, -314573, -301466, -288359, -275252, -262144, -249037, -235930, -222823, -209716,
26   -196608, -183501, -170394, -157287, -144180, -131072, -117965, -104858, -91751, -78644,
27   -65536, -52429, -39322, -26215, -13108,  0,  13107,  26214,  39321,  52428,
28   65536,  78643,  91750,  104857,  117964,  131072,  144179,  157286,  170393,  183500,
29   196608,  209715,  222822,  235929,  249036,  262144,  275251,  288358,  301465,  314572,
30   327680};
31 
32 
33 static const int kCdfSlopeQ0[51] = {  /* Q0 */
34   5,  5,  5,  5,  5,  5,  5,  5,  5,  5,
35   5,  5,  13,  23,  47,  87,  154,  315,  700,  1088,
36   2471,  6064,  14221,  21463,  36634,  36924,  19750,  13270,  5806,  2312,
37   1095,  660,  316,  145,  86,  41,  32,  5,  5,  5,
38   5,  5,  5,  5,  5,  5,  5,  5,  5,  2, 0};
39 
40 
41 static const int kCdfQ16[51] = {  /* Q16 */
42   0,  2,  4,  6,  8,  10,  12,  14,  16,  18,
43   20,  22,  24,  29,  38,  57,  92,  153,  279,  559,
44   994,  1983,  4408,  10097,  18682,  33336,  48105,  56005,  61313,  63636,
45   64560,  64998,  65262,  65389,  65447,  65481,  65497,  65510,  65512,  65514,
46   65516,  65518,  65520,  65522,  65524,  65526,  65528,  65530,  65532,  65534,
47   65535};
48 
49 
50 
51 /* function to be converted to fixed point */
piecewise(int32_t xinQ15)52 static __inline uint32_t piecewise(int32_t xinQ15) {
53 
54   int32_t ind, qtmp1, qtmp2, qtmp3;
55   uint32_t tmpUW32;
56 
57 
58   qtmp2 = xinQ15;
59 
60   if (qtmp2 < kHistEdgesQ15[0]) {
61     qtmp2 = kHistEdgesQ15[0];
62   }
63   if (qtmp2 > kHistEdgesQ15[50]) {
64     qtmp2 = kHistEdgesQ15[50];
65   }
66 
67   qtmp1 = qtmp2 - kHistEdgesQ15[0];       /* Q15 - Q15 = Q15        */
68   ind = (qtmp1 * 5) >> 16;              /* 2^16 / 5 = 0.4 in Q15  */
69   /* Q15 -> Q0              */
70   qtmp1 = qtmp2 - kHistEdgesQ15[ind];     /* Q15 - Q15 = Q15        */
71   qtmp2 = kCdfSlopeQ0[ind] * qtmp1;      /* Q0 * Q15 = Q15         */
72   qtmp3 = qtmp2>>15;                    /* Q15 -> Q0              */
73 
74   tmpUW32 = kCdfQ16[ind] + qtmp3;    /* Q0 + Q0 = Q0           */
75   return tmpUW32;
76 }
77 
78 
79 
WebRtcIsac_EncLogisticMulti2(Bitstr * streamdata,int16_t * dataQ7,const uint16_t * envQ8,const int N,const int16_t isSWB12kHz)80 int WebRtcIsac_EncLogisticMulti2(
81     Bitstr *streamdata,      /* in-/output struct containing bitstream */
82     int16_t *dataQ7,    /* input: data vector */
83     const uint16_t *envQ8, /* input: side info vector defining the width of the pdf */
84     const int N,       /* input: data vector length / 2 */
85     const int16_t isSWB12kHz)
86 {
87   uint32_t W_lower, W_upper;
88   uint32_t W_upper_LSB, W_upper_MSB;
89   uint8_t *stream_ptr;
90   uint8_t *maxStreamPtr;
91   uint8_t *stream_ptr_carry;
92   uint32_t cdf_lo, cdf_hi;
93   int k;
94 
95   /* point to beginning of stream buffer */
96   stream_ptr = streamdata->stream + streamdata->stream_index;
97   W_upper = streamdata->W_upper;
98 
99   maxStreamPtr = streamdata->stream + STREAM_SIZE_MAX_60 - 1;
100   for (k = 0; k < N; k++)
101   {
102     /* compute cdf_lower and cdf_upper by evaluating the piecewise linear cdf */
103     cdf_lo = piecewise((*dataQ7 - 64) * *envQ8);
104     cdf_hi = piecewise((*dataQ7 + 64) * *envQ8);
105 
106     /* test and clip if probability gets too small */
107     while (cdf_lo+1 >= cdf_hi) {
108       /* clip */
109       if (*dataQ7 > 0) {
110         *dataQ7 -= 128;
111         cdf_hi = cdf_lo;
112         cdf_lo = piecewise((*dataQ7 - 64) * *envQ8);
113       } else {
114         *dataQ7 += 128;
115         cdf_lo = cdf_hi;
116         cdf_hi = piecewise((*dataQ7 + 64) * *envQ8);
117       }
118     }
119 
120     dataQ7++;
121     // increment only once per 4 iterations for SWB-16kHz or WB
122     // increment only once per 2 iterations for SWB-12kHz
123     envQ8 += (isSWB12kHz)? (k & 1):((k & 1) & (k >> 1));
124 
125 
126     /* update interval */
127     W_upper_LSB = W_upper & 0x0000FFFF;
128     W_upper_MSB = W_upper >> 16;
129     W_lower = W_upper_MSB * cdf_lo;
130     W_lower += (W_upper_LSB * cdf_lo) >> 16;
131     W_upper = W_upper_MSB * cdf_hi;
132     W_upper += (W_upper_LSB * cdf_hi) >> 16;
133 
134     /* shift interval such that it begins at zero */
135     W_upper -= ++W_lower;
136 
137     /* add integer to bitstream */
138     streamdata->streamval += W_lower;
139 
140     /* handle carry */
141     if (streamdata->streamval < W_lower)
142     {
143       /* propagate carry */
144       stream_ptr_carry = stream_ptr;
145       while (!(++(*--stream_ptr_carry)));
146     }
147 
148     /* renormalize interval, store most significant byte of streamval and update streamval */
149     while ( !(W_upper & 0xFF000000) )      /* W_upper < 2^24 */
150     {
151       W_upper <<= 8;
152       *stream_ptr++ = (uint8_t) (streamdata->streamval >> 24);
153 
154       if(stream_ptr > maxStreamPtr)
155       {
156         return -ISAC_DISALLOWED_BITSTREAM_LENGTH;
157       }
158       streamdata->streamval <<= 8;
159     }
160   }
161 
162   /* calculate new stream_index */
163   streamdata->stream_index = (int)(stream_ptr - streamdata->stream);
164   streamdata->W_upper = W_upper;
165 
166   return 0;
167 }
168 
169 
170 
WebRtcIsac_DecLogisticMulti2(int16_t * dataQ7,Bitstr * streamdata,const uint16_t * envQ8,const int16_t * ditherQ7,const int N,const int16_t isSWB12kHz)171 int WebRtcIsac_DecLogisticMulti2(
172     int16_t *dataQ7,       /* output: data vector */
173     Bitstr *streamdata,      /* in-/output struct containing bitstream */
174     const uint16_t *envQ8, /* input: side info vector defining the width of the pdf */
175     const int16_t *ditherQ7,/* input: dither vector */
176     const int N,         /* input: data vector length */
177     const int16_t isSWB12kHz)
178 {
179   uint32_t    W_lower, W_upper;
180   uint32_t    W_tmp;
181   uint32_t    W_upper_LSB, W_upper_MSB;
182   uint32_t    streamval;
183   const uint8_t *stream_ptr;
184   uint32_t    cdf_tmp;
185   int16_t     candQ7;
186   int             k;
187 
188   stream_ptr = streamdata->stream + streamdata->stream_index;
189   W_upper = streamdata->W_upper;
190   if (streamdata->stream_index == 0)   /* first time decoder is called for this stream */
191   {
192     /* read first word from bytestream */
193     streamval = *stream_ptr << 24;
194     streamval |= *++stream_ptr << 16;
195     streamval |= *++stream_ptr << 8;
196     streamval |= *++stream_ptr;
197   } else {
198     streamval = streamdata->streamval;
199   }
200 
201 
202   for (k = 0; k < N; k++)
203   {
204     /* find the integer *data for which streamval lies in [W_lower+1, W_upper] */
205     W_upper_LSB = W_upper & 0x0000FFFF;
206     W_upper_MSB = W_upper >> 16;
207 
208     /* find first candidate by inverting the logistic cdf */
209     candQ7 = - *ditherQ7 + 64;
210     cdf_tmp = piecewise(candQ7 * *envQ8);
211 
212     W_tmp = W_upper_MSB * cdf_tmp;
213     W_tmp += (W_upper_LSB * cdf_tmp) >> 16;
214     if (streamval > W_tmp)
215     {
216       W_lower = W_tmp;
217       candQ7 += 128;
218       cdf_tmp = piecewise(candQ7 * *envQ8);
219 
220       W_tmp = W_upper_MSB * cdf_tmp;
221       W_tmp += (W_upper_LSB * cdf_tmp) >> 16;
222       while (streamval > W_tmp)
223       {
224         W_lower = W_tmp;
225         candQ7 += 128;
226         cdf_tmp = piecewise(candQ7 * *envQ8);
227 
228         W_tmp = W_upper_MSB * cdf_tmp;
229         W_tmp += (W_upper_LSB * cdf_tmp) >> 16;
230 
231         /* error check */
232         if (W_lower == W_tmp) return -1;
233       }
234       W_upper = W_tmp;
235 
236       /* another sample decoded */
237       *dataQ7 = candQ7 - 64;
238     }
239     else
240     {
241       W_upper = W_tmp;
242       candQ7 -= 128;
243       cdf_tmp = piecewise(candQ7 * *envQ8);
244 
245       W_tmp = W_upper_MSB * cdf_tmp;
246       W_tmp += (W_upper_LSB * cdf_tmp) >> 16;
247       while ( !(streamval > W_tmp) )
248       {
249         W_upper = W_tmp;
250         candQ7 -= 128;
251         cdf_tmp = piecewise(candQ7 * *envQ8);
252 
253         W_tmp = W_upper_MSB * cdf_tmp;
254         W_tmp += (W_upper_LSB * cdf_tmp) >> 16;
255 
256         /* error check */
257         if (W_upper == W_tmp) return -1;
258       }
259       W_lower = W_tmp;
260 
261       /* another sample decoded */
262       *dataQ7 = candQ7 + 64;
263     }
264     ditherQ7++;
265     dataQ7++;
266     // increment only once per 4 iterations for SWB-16kHz or WB
267     // increment only once per 2 iterations for SWB-12kHz
268     envQ8 += (isSWB12kHz)? (k & 1):((k & 1) & (k >> 1));
269 
270     /* shift interval to start at zero */
271     W_upper -= ++W_lower;
272 
273     /* add integer to bitstream */
274     streamval -= W_lower;
275 
276     /* renormalize interval and update streamval */
277     while ( !(W_upper & 0xFF000000) )    /* W_upper < 2^24 */
278     {
279       /* read next byte from stream */
280       streamval = (streamval << 8) | *++stream_ptr;
281       W_upper <<= 8;
282     }
283   }
284 
285   streamdata->stream_index = (int)(stream_ptr - streamdata->stream);
286   streamdata->W_upper = W_upper;
287   streamdata->streamval = streamval;
288 
289   /* find number of bytes in original stream (determined by current interval width) */
290   if ( W_upper > 0x01FFFFFF )
291     return streamdata->stream_index - 2;
292   else
293     return streamdata->stream_index - 1;
294 }
295