Switch iteration over channels to the do{}while(); construct in order to inform the...
[opus.git] / libcelt / quant_bands.c
1 /* Copyright (c) 2007-2008 CSIRO
2    Copyright (c) 2007-2009 Xiph.Org Foundation
3    Written by Jean-Marc Valin */
4 /*
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8    
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11    
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15    
16    - Neither the name of the Xiph.org Foundation nor the names of its
17    contributors may be used to endorse or promote products derived from
18    this software without specific prior written permission.
19    
20    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 #ifdef HAVE_CONFIG_H
34 #include "config.h"
35 #endif
36
37 #include "quant_bands.h"
38 #include "laplace.h"
39 #include <math.h>
40 #include "os_support.h"
41 #include "arch.h"
42 #include "mathops.h"
43 #include "stack_alloc.h"
44
45 #ifdef FIXED_POINT
46 /* Mean energy in each band quantized in Q6 */
47 static const signed char eMeans[25] = {
48       103,100, 92, 85, 81,
49        77, 72, 70, 78, 75,
50        73, 71, 78, 74, 69,
51        72, 70, 74, 76, 71,
52        60, 60, 60, 60, 60
53 };
54 #else
55 /* Mean energy in each band quantized in Q6 and converted back to float */
56 static const celt_word16 eMeans[25] = {
57       6.437500f, 6.250000f, 5.750000f, 5.312500f, 5.062500f,
58       4.812500f, 4.500000f, 4.375000f, 4.875000f, 4.687500f,
59       4.562500f, 4.437500f, 4.875000f, 4.625000f, 4.312500f,
60       4.500000f, 4.375000f, 4.625000f, 4.750000f, 4.437500f,
61       3.750000f, 3.750000f, 3.750000f, 3.750000f, 3.750000f
62 };
63 #endif
64 /* prediction coefficients: 0.9, 0.8, 0.65, 0.5 */
65 #ifdef FIXED_POINT
66 static const celt_word16 pred_coef[4] = {29440, 26112, 21248, 16384};
67 static const celt_word16 beta_coef[4] = {30147, 22282, 12124, 6554};
68 #else
69 static const celt_word16 pred_coef[4] = {29440/32768., 26112/32768., 21248/32768., 16384/32768.};
70 static const celt_word16 beta_coef[4] = {30147/32768., 22282/32768., 12124/32768., 6554/32768.};
71 #endif
72
73 static int intra_decision(const celt_word16 *eBands, celt_word16 *oldEBands, int start, int end, int len, int C)
74 {
75    int c, i;
76    celt_word32 dist = 0;
77    c=0; do {
78       for (i=start;i<end;i++)
79       {
80          celt_word16 d = SHR16(SUB16(eBands[i+c*len], oldEBands[i+c*len]),2);
81          dist = MAC16_16(dist, d,d);
82       }
83    } while (++c<C);
84    return SHR32(dist,2*DB_SHIFT-4) > 2*C*(end-start);
85 }
86
87 #ifndef STATIC_MODES
88
89 celt_int16 *quant_prob_alloc(const CELTMode *m)
90 {
91    int i;
92    celt_int16 *prob;
93    prob = celt_alloc(4*m->nbEBands*sizeof(celt_int16));
94    if (prob==NULL)
95      return NULL;
96    for (i=0;i<m->nbEBands;i++)
97    {
98       prob[2*i] = 7000-i*200;
99       prob[2*i+1] = ec_laplace_get_start_freq(prob[2*i]);
100    }
101    for (i=0;i<m->nbEBands;i++)
102    {
103       prob[2*m->nbEBands+2*i] = 9000-i*220;
104       prob[2*m->nbEBands+2*i+1] = ec_laplace_get_start_freq(prob[2*m->nbEBands+2*i]);
105    }
106    return prob;
107 }
108
109 void quant_prob_free(const celt_int16 *freq)
110 {
111    celt_free((celt_int16*)freq);
112 }
113 #endif
114
115 static void quant_coarse_energy_impl(const CELTMode *m, int start, int end,
116       const celt_word16 *eBands, celt_word16 *oldEBands, int budget,
117       const celt_int16 *prob, celt_word16 *error, ec_enc *enc, int _C, int LM,
118       int intra, celt_word16 max_decay)
119 {
120    const int C = CHANNELS(_C);
121    int i, c;
122    celt_word32 prev[2] = {0,0};
123    celt_word16 coef;
124    celt_word16 beta;
125
126    ec_enc_bit_prob(enc, intra, 8192);
127    if (intra)
128    {
129       coef = 0;
130       prob += 2*m->nbEBands;
131       beta = QCONST16(.15f,15);
132    } else {
133       beta = beta_coef[LM];
134       coef = pred_coef[LM];
135    }
136
137    /* Encode at a fixed coarse resolution */
138    for (i=start;i<end;i++)
139    {
140       c=0;
141       do {
142          int bits_left;
143          int qi;
144          celt_word16 q;
145          celt_word16 x;
146          celt_word32 f;
147          x = eBands[i+c*m->nbEBands];
148 #ifdef FIXED_POINT
149          f = SHL32(EXTEND32(x),15) -MULT16_16(coef,oldEBands[i+c*m->nbEBands])-prev[c];
150          /* Rounding to nearest integer here is really important! */
151          qi = (f+QCONST32(.5,DB_SHIFT+15))>>(DB_SHIFT+15);
152 #else
153          f = x-coef*oldEBands[i+c*m->nbEBands]-prev[c];
154          /* Rounding to nearest integer here is really important! */
155          qi = (int)floor(.5f+f);
156 #endif
157          /* Prevent the energy from going down too quickly (e.g. for bands
158             that have just one bin) */
159          if (qi < 0 && x < oldEBands[i+c*m->nbEBands]-max_decay)
160          {
161             qi += (int)SHR16(oldEBands[i+c*m->nbEBands]-max_decay-x, DB_SHIFT);
162             if (qi > 0)
163                qi = 0;
164          }
165          /* If we don't have enough bits to encode all the energy, just assume something safe.
166             We allow slightly busting the budget here */
167          bits_left = budget-(int)ec_enc_tell(enc, 0)-2*C*(end-i);
168          if (bits_left < 24)
169          {
170             if (qi > 1)
171                qi = 1;
172             if (qi < -1)
173                qi = -1;
174             if (bits_left<8)
175                qi = 0;
176          }
177          ec_laplace_encode_start(enc, &qi, prob[2*i], prob[2*i+1]);
178          error[i+c*m->nbEBands] = PSHR32(f,15) - SHL16(qi,DB_SHIFT);
179          q = SHL16(qi,DB_SHIFT);
180          
181          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
182          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
183       } while (++c < C);
184    }
185 }
186
187 void quant_coarse_energy(const CELTMode *m, int start, int end, int effEnd,
188       const celt_word16 *eBands, celt_word16 *oldEBands, int budget,
189       const celt_int16 *prob, celt_word16 *error, ec_enc *enc, int _C, int LM,
190       int nbAvailableBytes, int force_intra, int *delayedIntra, int two_pass)
191 {
192    const int C = CHANNELS(_C);
193    int intra;
194    celt_word16 max_decay;
195    VARDECL(celt_word16, oldEBands_intra);
196    VARDECL(celt_word16, error_intra);
197    ec_enc enc_start_state;
198    ec_byte_buffer buf_start_state;
199    SAVE_STACK;
200
201    intra = force_intra || (*delayedIntra && nbAvailableBytes > end);
202    if (/*shortBlocks || */intra_decision(eBands, oldEBands, start, effEnd, m->nbEBands, C))
203       *delayedIntra = 1;
204    else
205       *delayedIntra = 0;
206
207    /* Encode the global flags using a simple probability model
208       (first symbols in the stream) */
209
210 #ifdef FIXED_POINT
211       max_decay = MIN32(QCONST16(16,DB_SHIFT), SHL32(EXTEND32(nbAvailableBytes),DB_SHIFT-3));
212 #else
213    max_decay = MIN32(16.f, .125f*nbAvailableBytes);
214 #endif
215
216    enc_start_state = *enc;
217    buf_start_state = *(enc->buf);
218
219    ALLOC(oldEBands_intra, C*m->nbEBands, celt_word16);
220    ALLOC(error_intra, C*m->nbEBands, celt_word16);
221    CELT_COPY(oldEBands_intra, oldEBands, C*end);
222
223    if (two_pass || intra)
224    {
225       quant_coarse_energy_impl(m, start, end, eBands, oldEBands_intra, budget,
226             prob, error_intra, enc, C, LM, 1, max_decay);
227    }
228
229    if (!intra)
230    {
231       ec_enc enc_intra_state;
232       ec_byte_buffer buf_intra_state;
233       int tell_intra;
234       VARDECL(unsigned char, intra_bits);
235
236       tell_intra = ec_enc_tell(enc, 3);
237
238       enc_intra_state = *enc;
239       buf_intra_state = *(enc->buf);
240
241       ALLOC(intra_bits, buf_intra_state.ptr-buf_start_state.ptr, unsigned char);
242       /* Copy bits from intra bit-stream */
243       CELT_COPY(intra_bits, buf_start_state.ptr, buf_intra_state.ptr-buf_start_state.ptr);
244
245       *enc = enc_start_state;
246       *(enc->buf) = buf_start_state;
247
248       quant_coarse_energy_impl(m, start, end, eBands, oldEBands, budget,
249             prob, error, enc, C, LM, 0, max_decay);
250
251       if (two_pass && ec_enc_tell(enc, 3) > tell_intra)
252       {
253          *enc = enc_intra_state;
254          *(enc->buf) = buf_intra_state;
255          /* Copy bits from to bit-stream */
256          CELT_COPY(buf_start_state.ptr, intra_bits, buf_intra_state.ptr-buf_start_state.ptr);
257          CELT_COPY(oldEBands, oldEBands_intra, C*end);
258          CELT_COPY(error, error_intra, C*end);
259       }
260    } else {
261       CELT_COPY(oldEBands, oldEBands_intra, C*end);
262       CELT_COPY(error, error_intra, C*end);
263    }
264    RESTORE_STACK;
265 }
266
267 void quant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, ec_enc *enc, int _C)
268 {
269    int i, c;
270    const int C = CHANNELS(_C);
271
272    /* Encode finer resolution */
273    for (i=start;i<end;i++)
274    {
275       celt_int16 frac = 1<<fine_quant[i];
276       if (fine_quant[i] <= 0)
277          continue;
278       c=0;
279       do {
280          int q2;
281          celt_word16 offset;
282 #ifdef FIXED_POINT
283          /* Has to be without rounding */
284          q2 = (error[i+c*m->nbEBands]+QCONST16(.5f,DB_SHIFT))>>(DB_SHIFT-fine_quant[i]);
285 #else
286          q2 = (int)floor((error[i+c*m->nbEBands]+.5f)*frac);
287 #endif
288          if (q2 > frac-1)
289             q2 = frac-1;
290          if (q2<0)
291             q2 = 0;
292          ec_enc_bits(enc, q2, fine_quant[i]);
293 #ifdef FIXED_POINT
294          offset = SUB16(SHR32(SHL32(EXTEND32(q2),DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
295 #else
296          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
297 #endif
298          oldEBands[i+c*m->nbEBands] += offset;
299          error[i+c*m->nbEBands] -= offset;
300          /*printf ("%f ", error[i] - offset);*/
301       } while (++c < C);
302    }
303 }
304
305 void quant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, celt_word16 *error, int *fine_quant, int *fine_priority, int bits_left, ec_enc *enc, int _C)
306 {
307    int i, prio, c;
308    const int C = CHANNELS(_C);
309
310    /* Use up the remaining bits */
311    for (prio=0;prio<2;prio++)
312    {
313       for (i=start;i<end && bits_left>=C ;i++)
314       {
315          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
316             continue;
317          c=0;
318          do {
319             int q2;
320             celt_word16 offset;
321             q2 = error[i+c*m->nbEBands]<0 ? 0 : 1;
322             ec_enc_bits(enc, q2, 1);
323 #ifdef FIXED_POINT
324             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
325 #else
326             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
327 #endif
328             oldEBands[i+c*m->nbEBands] += offset;
329             bits_left--;
330          } while (++c < C);
331       }
332    }
333 }
334
335 void unquant_coarse_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int intra, const celt_int16 *prob, ec_dec *dec, int _C, int LM)
336 {
337    int i, c;
338    celt_word32 prev[2] = {0, 0};
339    celt_word16 coef;
340    celt_word16 beta;
341    const int C = CHANNELS(_C);
342
343
344    if (intra)
345    {
346       coef = 0;
347       beta = QCONST16(.15f,15);
348       prob += 2*m->nbEBands;
349    } else {
350       beta = beta_coef[LM];
351       coef = pred_coef[LM];
352    }
353
354    /* Decode at a fixed coarse resolution */
355    for (i=start;i<end;i++)
356    {
357       c=0;
358       do {
359          int qi;
360          celt_word16 q;
361          qi = ec_laplace_decode_start(dec, prob[2*i], prob[2*i+1]);
362          q = SHL16(qi,DB_SHIFT);
363
364          oldEBands[i+c*m->nbEBands] = PSHR32(MULT16_16(coef,oldEBands[i+c*m->nbEBands]) + prev[c] + SHL32(EXTEND32(q),15), 15);
365          prev[c] = prev[c] + SHL32(EXTEND32(q),15) - MULT16_16(beta,q);
366       } while (++c < C);
367    }
368 }
369
370 void unquant_fine_energy(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant, ec_dec *dec, int _C)
371 {
372    int i, c;
373    const int C = CHANNELS(_C);
374    /* Decode finer resolution */
375    for (i=start;i<end;i++)
376    {
377       if (fine_quant[i] <= 0)
378          continue;
379       c=0; 
380       do {
381          int q2;
382          celt_word16 offset;
383          q2 = ec_dec_bits(dec, fine_quant[i]);
384 #ifdef FIXED_POINT
385          offset = SUB16(SHR32(SHL32(EXTEND32(q2),DB_SHIFT)+QCONST16(.5,DB_SHIFT),fine_quant[i]),QCONST16(.5f,DB_SHIFT));
386 #else
387          offset = (q2+.5f)*(1<<(14-fine_quant[i]))*(1.f/16384) - .5f;
388 #endif
389          oldEBands[i+c*m->nbEBands] += offset;
390       } while (++c < C);
391    }
392 }
393
394 void unquant_energy_finalise(const CELTMode *m, int start, int end, celt_ener *eBands, celt_word16 *oldEBands, int *fine_quant,  int *fine_priority, int bits_left, ec_dec *dec, int _C)
395 {
396    int i, prio, c;
397    const int C = CHANNELS(_C);
398
399    /* Use up the remaining bits */
400    for (prio=0;prio<2;prio++)
401    {
402       for (i=start;i<end && bits_left>=C ;i++)
403       {
404          if (fine_quant[i] >= 7 || fine_priority[i]!=prio)
405             continue;
406          c=0;
407          do {
408             int q2;
409             celt_word16 offset;
410             q2 = ec_dec_bits(dec, 1);
411 #ifdef FIXED_POINT
412             offset = SHR16(SHL16(q2,DB_SHIFT)-QCONST16(.5,DB_SHIFT),fine_quant[i]+1);
413 #else
414             offset = (q2-.5f)*(1<<(14-fine_quant[i]-1))*(1.f/16384);
415 #endif
416             oldEBands[i+c*m->nbEBands] += offset;
417             bits_left--;
418          } while (++c < C);
419       }
420    }
421 }
422
423 void log2Amp(const CELTMode *m, int start, int end,
424       celt_ener *eBands, celt_word16 *oldEBands, int _C)
425 {
426    int c, i;
427    const int C = CHANNELS(_C);
428    c=0;
429    do {
430       for (i=start;i<m->nbEBands;i++)
431       {
432          celt_word16 lg = oldEBands[i+c*m->nbEBands]
433                         + SHL16((celt_word16)eMeans[i],6);
434          eBands[i+c*m->nbEBands] = PSHR32(celt_exp2(SHL16(lg,11-DB_SHIFT)),4);
435          if (oldEBands[i+c*m->nbEBands] < -QCONST16(14.f,DB_SHIFT))
436             oldEBands[i+c*m->nbEBands] = -QCONST16(14.f,DB_SHIFT);
437       }
438    } while (++c < C);
439 }
440
441 void amp2Log2(const CELTMode *m, int effEnd, int end,
442       celt_ener *bandE, celt_word16 *bandLogE, int _C)
443 {
444    int c, i;
445    const int C = CHANNELS(_C);
446    c=0;
447    do {
448       for (i=0;i<effEnd;i++)
449          bandLogE[i+c*m->nbEBands] =
450                celt_log2(MAX32(QCONST32(.001f,14),SHL32(bandE[i+c*m->nbEBands],2)))
451                - SHL16((celt_word16)eMeans[i],6);
452       for (i=effEnd;i<end;i++)
453          bandLogE[c*m->nbEBands+i] = -QCONST16(14.f,DB_SHIFT);
454    } while (++c < C);
455 }