Initial support for a managed stack/scratchpad. Still needs some work.
[opus.git] / libcelt / quant_bands.c
1 /* (C) 2007-2008 Jean-Marc Valin, CSIRO
2 */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7    
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10    
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14    
15    - Neither the name of the Xiph.org Foundation nor the names of its
16    contributors may be used to endorse or promote products derived from
17    this software without specific prior written permission.
18    
19    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
23    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include "quant_bands.h"
37 #include "laplace.h"
38 #include <math.h>
39 #include "os_support.h"
40 #include "arch.h"
41 #include "mathops.h"
42
43 #ifdef FIXED_POINT
44 const celt_word16_t eMeans[24] = {11520, -2048, -3072, -640, 256, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
45 #else
46 const celt_word16_t eMeans[24] = {45.f, -8.f, -12.f, -2.5f, 1.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
47 #endif
48
49 /*const int frac[24] = {4, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};*/
50 const int frac[24] = {8, 6, 5, 4, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2};
51
52 #ifdef FIXED_POINT
53 static inline celt_ener_t dB2Amp(celt_ener_t dB)
54 {
55    celt_ener_t amp;
56    amp = PSHR32(celt_exp2(MULT16_16_Q14(21771,dB)),2)-QCONST16(.3f, 14);
57    if (amp < 0)
58       amp = 0;
59    return amp;
60 }
61
62 #define DBofTWO 24661
63 static inline celt_word16_t amp2dB(celt_ener_t amp)
64 {
65    /* equivalent to return 6.0207*log2(.3+amp) */
66    return ROUND(MULT16_16(24661,celt_log2(ADD32(QCONST32(.3f,14),amp))),12);
67    /* return DB_SCALING*20*log10(.3+ENER_SCALING_1*amp); */
68 }
69 #else
70 static inline celt_ener_t dB2Amp(celt_ener_t dB)
71 {
72    celt_ener_t amp;
73    amp = pow(10, .05*dB)-.3;
74    if (amp < 0)
75       amp = 0;
76    return amp;
77 }
78 static inline celt_word16_t amp2dB(celt_ener_t amp)
79 {
80    return 20*log10(.3+amp);
81 }
82 #endif
83
84 static const celt_word16_t base_resolution = QCONST16(6.f,8);
85
86 static void quant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
87 {
88    int i;
89    int bits;
90    celt_word16_t prev = 0;
91    celt_word16_t coef = m->ePredCoef;
92    VARDECL(celt_word16_t *error);
93    SAVE_STACK;
94    /* The .7 is a heuristic */
95    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
96    
97    ALLOC(error, m->nbEBands, celt_word16_t);
98    bits = ec_enc_tell(enc, 0);
99    for (i=0;i<m->nbEBands;i++)
100    {
101       int qi;
102       celt_word16_t q;   /* dB */
103       celt_word16_t x;   /* dB */
104       celt_word16_t f;   /* Q8 */
105       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
106       x = amp2dB(eBands[i]);
107       f = DIV32_16(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution);
108 #ifdef FIXED_POINT
109       /* Rounding to nearest integer here is really important! */
110       qi = (f+128)>>8;
111 #else
112       qi = (int)floor(.5+f);
113 #endif
114       /*ec_laplace_encode(enc, qi, i==0?11192:6192);*/
115       /*ec_laplace_encode(enc, qi, 8500-i*200);*/
116       /* If we don't have enough bits to encode all the energy, just assume something safe. */
117       if (ec_enc_tell(enc, 0) - bits > budget)
118          qi = -1;
119       else
120          ec_laplace_encode(enc, qi, 6000-i*200);
121       q = qi*base_resolution;
122       error[i] = f - SHL16(qi,8);
123       
124       /*printf("%d ", qi);*/
125       /*printf("%f %f ", pred+prev+q, x);*/
126       /*printf("%f ", x-pred);*/
127       
128       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
129       
130       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
131    }
132    /*bits = ec_enc_tell(enc, 0) - bits;*/
133    /*printf ("%d\n", bits);*/
134    for (i=0;i<m->nbEBands;i++)
135    {
136       int q2;
137       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
138       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
139       if (ec_enc_tell(enc, 0) - bits +EC_ILOG(frac[i])> budget)
140          break;
141 #ifdef FIXED_POINT
142       /* Has to be without rounding */
143       q2 = offset>>8;
144 #else
145       q2 = (int)floor(offset);
146 #endif
147       if (q2 > frac[i]-1)
148          q2 = frac[i]-1;
149       ec_enc_uint(enc, q2, frac[i]);
150       offset = DIV32_16(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8);
151       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
152       /*printf ("%f ", error[i] - offset);*/
153    }
154    for (i=0;i<m->nbEBands;i++)
155    {
156       eBands[i] = dB2Amp(oldEBands[i]);
157    }
158    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
159
160    /*printf ("\n");*/
161    RESTORE_STACK;
162 }
163
164 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
165 {
166    int i;
167    int bits;
168    celt_word16_t prev = 0;
169    celt_word16_t coef = m->ePredCoef;
170    /* The .7 is a heuristic */
171    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
172    bits = ec_dec_tell(dec, 0);
173    for (i=0;i<m->nbEBands;i++)
174    {
175       int qi;
176       celt_word16_t q;
177       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
178       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
179       if (ec_dec_tell(dec, 0) - bits > budget)
180          qi = -1;
181       else
182          qi = ec_laplace_decode(dec, 6000-i*200);
183       q = qi*base_resolution;
184       
185       /*printf("%d ", qi);*/
186       /*printf("%f %f ", pred+prev+q, x);*/
187       /*printf("%f ", x-pred);*/
188       
189       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
190       
191       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
192    }
193    for (i=0;i<m->nbEBands;i++)
194    {
195       int q2;
196       celt_word16_t offset;
197       if (ec_dec_tell(dec, 0) - bits +EC_ILOG(frac[i])> budget)
198          break;
199       q2 = ec_dec_uint(dec, frac[i]);
200       offset = DIV32_16(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8);
201       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
202    }
203    for (i=0;i<m->nbEBands;i++)
204    {
205       eBands[i] = dB2Amp(oldEBands[i]);
206    }
207    /*printf ("\n");*/
208 }
209
210
211
212 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
213 {
214    int C;
215    SAVE_STACK;
216    
217    C = m->nbChannels;
218
219    if (C==1)
220       quant_energy_mono(m, eBands, oldEBands, budget, enc);
221    else 
222 #if 1
223    {
224       int c;
225       VARDECL(celt_ener_t *E);
226       ALLOC(E, m->nbEBands, celt_ener_t);
227       for (c=0;c<C;c++)
228       {
229          int i;
230          for (i=0;i<m->nbEBands;i++)
231             E[i] = eBands[C*i+c];
232          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
233          for (i=0;i<m->nbEBands;i++)
234             eBands[C*i+c] = E[i];
235       }
236    }
237 #else
238       if (C==2)
239    {
240       int i;
241       int NB = m->nbEBands;
242       celt_ener_t mid[NB];
243       celt_ener_t side[NB];
244       for (i=0;i<NB;i++)
245       {
246          //left = eBands[C*i];
247          //right = eBands[C*i+1];
248          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
249          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
250          //printf ("%f %f ", mid[i], side[i]);
251       }
252       //printf ("\n");
253       quant_energy_mono(m, mid, oldEBands, enc);
254       for (i=0;i<NB;i++)
255          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
256          
257       //quant_energy_side(m, side, oldEBands+NB, enc);
258       for (i=0;i<NB;i++)
259       {
260          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
261          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
262          //printf ("%f %f ", mid[i], side[i]);
263       }
264
265    } else {
266       celt_fatal("more than 2 channels not supported");
267    }
268 #endif
269    RESTORE_STACK;
270 }
271
272
273
274 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
275 {
276    int C;   
277    SAVE_STACK;
278    C = m->nbChannels;
279
280    if (C==1)
281       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
282    else {
283       int c;
284       VARDECL(celt_ener_t *E);
285       ALLOC(E, m->nbEBands, celt_ener_t);
286       for (c=0;c<C;c++)
287       {
288          int i;
289          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
290          for (i=0;i<m->nbEBands;i++)
291             eBands[C*i+c] = E[i];
292       }
293    }
294    RESTORE_STACK;
295 }