1 /*
2  * Siren Encoder/Decoder library
3  *
4  *   @author: Youness Alaoui <kakaroto@kakaroto.homelinux.net>
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Library General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Library General Public License for more details.
15  *
16  * You should have received a copy of the GNU Library General Public
17  * License along with this library; if not, write to the
18  * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
19  * Boston, MA 02110-1301, USA.
20  */
21 
22 
23 #include "siren7.h"
24 #include "huffman_consts.h"
25 
26 
27 static short current_word = 0;
28 static int bit_idx = 0;
29 static int *bitstream_ptr = NULL;
30 
31 int
next_bit(void)32 next_bit (void)
33 {
34   if (bitstream_ptr == NULL)
35     return -1;
36 
37   if (bit_idx == 0) {
38     current_word = *bitstream_ptr++;
39     bit_idx = 16;
40   }
41 
42   return (current_word >> --bit_idx) & 1;
43 }
44 
45 void
set_bitstream(int * stream)46 set_bitstream (int *stream)
47 {
48   bitstream_ptr = stream;
49   current_word = *bitstream_ptr;
50   bit_idx = 0;
51 }
52 
53 
54 int
compute_region_powers(int number_of_regions,float * coefs,int * drp_num_bits,int * drp_code_bits,int * absolute_region_power_index,int esf_adjustment)55 compute_region_powers (int number_of_regions, float *coefs, int *drp_num_bits,
56     int *drp_code_bits, int *absolute_region_power_index, int esf_adjustment)
57 {
58   float region_power = 0;
59   int num_bits;
60   int idx;
61   int max_idx, min_idx;
62   int region, i;
63 
64   for (region = 0; region < number_of_regions; region++) {
65     region_power = 0.0f;
66     for (i = 0; i < region_size; i++) {
67       region_power +=
68           coefs[(region * region_size) + i] * coefs[(region * region_size) + i];
69     }
70     region_power *= region_size_inverse;
71 
72     min_idx = 0;
73     max_idx = 64;
74     for (i = 0; i < 6; i++) {
75       idx = (min_idx + max_idx) / 2;
76       if (region_power_table_boundary[idx - 1] <= region_power) {
77         min_idx = idx;
78       } else {
79         max_idx = idx;
80       }
81     }
82     absolute_region_power_index[region] = min_idx - 24;
83 
84   }
85 
86   for (region = number_of_regions - 2; region >= 0; region--) {
87     if (absolute_region_power_index[region] <
88         absolute_region_power_index[region + 1] - 11)
89       absolute_region_power_index[region] =
90           absolute_region_power_index[region + 1] - 11;
91   }
92 
93   if (absolute_region_power_index[0] < (1 - esf_adjustment))
94     absolute_region_power_index[0] = (1 - esf_adjustment);
95 
96   if (absolute_region_power_index[0] > (31 - esf_adjustment))
97     absolute_region_power_index[0] = (31 - esf_adjustment);
98 
99   drp_num_bits[0] = 5;
100   drp_code_bits[0] = absolute_region_power_index[0] + esf_adjustment;
101 
102 
103   for (region = 1; region < number_of_regions; region++) {
104     if (absolute_region_power_index[region] < (-8 - esf_adjustment))
105       absolute_region_power_index[region] = (-8 - esf_adjustment);
106     if (absolute_region_power_index[region] > (31 - esf_adjustment))
107       absolute_region_power_index[region] = (31 - esf_adjustment);
108   }
109 
110   num_bits = 5;
111 
112   for (region = 0; region < number_of_regions - 1; region++) {
113     idx =
114         absolute_region_power_index[region + 1] -
115         absolute_region_power_index[region] + 12;
116     if (idx < 0)
117       idx = 0;
118 
119     absolute_region_power_index[region + 1] =
120         absolute_region_power_index[region] + idx - 12;
121     drp_num_bits[region + 1] = differential_region_power_bits[region][idx];
122     drp_code_bits[region + 1] = differential_region_power_codes[region][idx];
123     num_bits += drp_num_bits[region + 1];
124   }
125 
126   return num_bits;
127 }
128 
129 
130 int
decode_envelope(int number_of_regions,float * decoder_standard_deviation,int * absolute_region_power_index,int esf_adjustment)131 decode_envelope (int number_of_regions, float *decoder_standard_deviation,
132     int *absolute_region_power_index, int esf_adjustment)
133 {
134   int index;
135   int i;
136   int envelope_bits = 0;
137 
138   index = 0;
139   for (i = 0; i < 5; i++)
140     index = (index << 1) | next_bit ();
141   envelope_bits = 5;
142 
143   absolute_region_power_index[0] = index - esf_adjustment;
144   decoder_standard_deviation[0] =
145       standard_deviation[absolute_region_power_index[0] + 24];
146 
147   for (i = 1; i < number_of_regions; i++) {
148     index = 0;
149     do {
150       index = differential_decoder_tree[i - 1][index][next_bit ()];
151       envelope_bits++;
152     } while (index > 0);
153 
154     absolute_region_power_index[i] =
155         absolute_region_power_index[i - 1] - index - 12;
156     if (absolute_region_power_index[i] < -24)
157       absolute_region_power_index[i] = -24;
158     else if (absolute_region_power_index[i] > 39)
159       absolute_region_power_index[i] = 39;
160     decoder_standard_deviation[i] =
161         standard_deviation[absolute_region_power_index[i] + 24];
162   }
163 
164   return envelope_bits;
165 }
166 
167 
168 
169 static int
huffman_vector(int category,int power_idx,float * mlts,int * out)170 huffman_vector (int category, int power_idx, float *mlts, int *out)
171 {
172   int i, j;
173   float temp_value = deviation_inverse[power_idx] * step_size_inverse[category];
174   int sign_idx, idx, non_zeroes, max, bits_available;
175   int current_word = 0;
176   int region_bits = 0;
177 
178   bits_available = 32;
179   for (i = 0; i < number_of_vectors[category]; i++) {
180     sign_idx = idx = non_zeroes = 0;
181     for (j = 0; j < vector_dimension[category]; j++) {
182       max = (int) ((fabs (*mlts) * temp_value) + dead_zone[category]);
183       if (max != 0) {
184         sign_idx <<= 1;
185         non_zeroes++;
186         if (*mlts > 0)
187           sign_idx++;
188         if (max > max_bin[category] || max < 0)
189           max = max_bin[category];
190 
191       }
192       mlts++;
193       idx = (idx * (max_bin[category] + 1)) + max;
194     }
195 
196     region_bits += bitcount_tables[category][idx] + non_zeroes;
197     bits_available -= bitcount_tables[category][idx] + non_zeroes;
198     if (bits_available < 0) {
199       *out++ =
200           current_word + (((code_tables[category][idx] << non_zeroes) +
201               sign_idx) >> -bits_available);
202       bits_available += 32;
203       current_word =
204           ((code_tables[category][idx] << non_zeroes) +
205           sign_idx) << bits_available;
206     } else {
207       current_word +=
208           ((code_tables[category][idx] << non_zeroes) +
209           sign_idx) << bits_available;
210     }
211 
212   }
213 
214   *out = current_word;
215   return region_bits;
216 }
217 
218 int
quantize_mlt(int number_of_regions,int rate_control_possibilities,int number_of_available_bits,float * coefs,int * absolute_region_power_index,int * power_categories,int * category_balance,int * region_mlt_bit_counts,int * region_mlt_bits)219 quantize_mlt (int number_of_regions, int rate_control_possibilities,
220     int number_of_available_bits, float *coefs,
221     int *absolute_region_power_index, int *power_categories,
222     int *category_balance, int *region_mlt_bit_counts, int *region_mlt_bits)
223 {
224   int region;
225   int mlt_bits = 0;
226   int rate_control;
227 
228   for (rate_control = 0; rate_control < ((rate_control_possibilities >> 1) - 1);
229       rate_control++)
230     power_categories[category_balance[rate_control]]++;
231 
232   for (region = 0; region < number_of_regions; region++) {
233     if (power_categories[region] > 6)
234       region_mlt_bit_counts[region] = 0;
235     else
236       region_mlt_bit_counts[region] =
237           huffman_vector (power_categories[region],
238           absolute_region_power_index[region], coefs + (region_size * region),
239           region_mlt_bits + (4 * region));
240     mlt_bits += region_mlt_bit_counts[region];
241   }
242 
243   while (mlt_bits < number_of_available_bits && rate_control > 0) {
244     rate_control--;
245     region = category_balance[rate_control];
246     power_categories[region]--;
247 
248     if (power_categories[region] < 0)
249       power_categories[region] = 0;
250 
251     mlt_bits -= region_mlt_bit_counts[region];
252 
253     if (power_categories[region] > 6)
254       region_mlt_bit_counts[region] = 0;
255     else
256       region_mlt_bit_counts[region] =
257           huffman_vector (power_categories[region],
258           absolute_region_power_index[region], coefs + (region_size * region),
259           region_mlt_bits + (4 * region));
260 
261     mlt_bits += region_mlt_bit_counts[region];
262   }
263 
264   while (mlt_bits > number_of_available_bits
265       && rate_control < rate_control_possibilities) {
266     region = category_balance[rate_control];
267     power_categories[region]++;
268     mlt_bits -= region_mlt_bit_counts[region];
269 
270     if (power_categories[region] > 6)
271       region_mlt_bit_counts[region] = 0;
272     else
273       region_mlt_bit_counts[region] =
274           huffman_vector (power_categories[region],
275           absolute_region_power_index[region], coefs + (region_size * region),
276           region_mlt_bits + (4 * region));
277 
278     mlt_bits += region_mlt_bit_counts[region];
279 
280     rate_control++;
281   }
282 
283   return rate_control;
284 }
285 
286 static int
get_dw(SirenDecoder decoder)287 get_dw (SirenDecoder decoder)
288 {
289   int ret = decoder->dw1 + decoder->dw4;
290 
291   if ((ret & 0x8000) != 0)
292     ret++;
293 
294   decoder->dw1 = decoder->dw2;
295   decoder->dw2 = decoder->dw3;
296   decoder->dw3 = decoder->dw4;
297   decoder->dw4 = ret;
298 
299   return ret;
300 }
301 
302 
303 
304 
305 int
decode_vector(SirenDecoder decoder,int number_of_regions,int number_of_available_bits,float * decoder_standard_deviation,int * power_categories,float * coefs,int scale_factor)306 decode_vector (SirenDecoder decoder, int number_of_regions,
307     int number_of_available_bits, float *decoder_standard_deviation,
308     int *power_categories, float *coefs, int scale_factor)
309 {
310   float *coefs_ptr;
311   float decoded_value;
312   float noise;
313   int *decoder_tree;
314 
315   int region;
316   int category;
317   int i, j;
318   int index;
319   int error;
320   int dw1;
321   int dw2;
322 
323   error = 0;
324   for (region = 0; region < number_of_regions; region++) {
325     category = power_categories[region];
326     coefs_ptr = coefs + (region * region_size);
327 
328     if (category < 7) {
329       decoder_tree = decoder_tables[category];
330 
331       for (i = 0; i < number_of_vectors[category]; i++) {
332         index = 0;
333         do {
334           if (number_of_available_bits <= 0) {
335             error = 1;
336             break;
337           }
338 
339           index = decoder_tree[index + next_bit ()];
340           number_of_available_bits--;
341         } while ((index & 1) == 0);
342 
343         index >>= 1;
344 
345         if (error == 0 && number_of_available_bits >= 0) {
346           for (j = 0; j < vector_dimension[category]; j++) {
347             decoded_value =
348                 mlt_quant[category][index & ((1 << index_table[category]) - 1)];
349             index >>= index_table[category];
350 
351             if (decoded_value != 0) {
352               if (next_bit () == 0)
353                 decoded_value *= -decoder_standard_deviation[region];
354               else
355                 decoded_value *= decoder_standard_deviation[region];
356               number_of_available_bits--;
357             }
358 
359             *coefs_ptr++ = decoded_value * scale_factor;
360           }
361         } else {
362           error = 1;
363           break;
364         }
365       }
366 
367       if (error == 1) {
368         for (j = region + 1; j < number_of_regions; j++)
369           power_categories[j] = 7;
370         category = 7;
371       }
372     }
373 
374 
375     coefs_ptr = coefs + (region * region_size);
376 
377     if (category == 5) {
378       i = 0;
379       for (j = 0; j < region_size; j++) {
380         if (*coefs_ptr != 0) {
381           i++;
382           if (fabs (*coefs_ptr) > 2.0 * decoder_standard_deviation[region]) {
383             i += 3;
384           }
385         }
386         coefs_ptr++;
387       }
388 
389       noise = decoder_standard_deviation[region] * noise_category5[i];
390     } else if (category == 6) {
391       i = 0;
392       for (j = 0; j < region_size; j++) {
393         if (*coefs_ptr++ != 0)
394           i++;
395       }
396 
397       noise = decoder_standard_deviation[region] * noise_category6[i];
398     } else if (category == 7) {
399       noise = decoder_standard_deviation[region] * noise_category7;
400     } else {
401       noise = 0;
402     }
403 
404     coefs_ptr = coefs + (region * region_size);
405 
406     if (category == 5 || category == 6 || category == 7) {
407       dw1 = get_dw (decoder);
408       dw2 = get_dw (decoder);
409 
410       for (j = 0; j < 10; j++) {
411         if (category == 7 || *coefs_ptr == 0) {
412           if ((dw1 & 1))
413             *coefs_ptr = noise;
414           else
415             *coefs_ptr = -noise;
416         }
417         coefs_ptr++;
418         dw1 >>= 1;
419 
420         if (category == 7 || *coefs_ptr == 0) {
421           if ((dw2 & 1))
422             *coefs_ptr = noise;
423           else
424             *coefs_ptr = -noise;
425         }
426         coefs_ptr++;
427         dw2 >>= 1;
428       }
429     }
430   }
431 
432   return error == 1 ? -1 : number_of_available_bits;
433 }
434