Another C90-fying pass. Fixed some warnings in the process.
[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    celt_word16_t beta;
93    VARDECL(celt_word16_t *error);
94    SAVE_STACK;
95    /* The .7 is a heuristic */
96    beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
97    
98    ALLOC(error, m->nbEBands, celt_word16_t);
99    bits = ec_enc_tell(enc, 0);
100    for (i=0;i<m->nbEBands;i++)
101    {
102       int qi;
103       celt_word16_t q;   /* dB */
104       celt_word16_t x;   /* dB */
105       celt_word16_t f;   /* Q8 */
106       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
107       x = amp2dB(eBands[i]);
108       f = DIV32_16(SHL32(EXTEND32(x-mean-MULT16_16_Q15(coef,oldEBands[i])-prev),8),base_resolution);
109 #ifdef FIXED_POINT
110       /* Rounding to nearest integer here is really important! */
111       qi = (f+128)>>8;
112 #else
113       qi = (int)floor(.5+f);
114 #endif
115       /*ec_laplace_encode(enc, qi, i==0?11192:6192);*/
116       /*ec_laplace_encode(enc, qi, 8500-i*200);*/
117       /* If we don't have enough bits to encode all the energy, just assume something safe. */
118       if (ec_enc_tell(enc, 0) - bits > budget)
119          qi = -1;
120       else
121          ec_laplace_encode(enc, qi, 6000-i*200);
122       q = qi*base_resolution;
123       error[i] = f - SHL16(qi,8);
124       
125       /*printf("%d ", qi);*/
126       /*printf("%f %f ", pred+prev+q, x);*/
127       /*printf("%f ", x-pred);*/
128       
129       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
130       
131       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
132    }
133    /*bits = ec_enc_tell(enc, 0) - bits;*/
134    /*printf ("%d\n", bits);*/
135    for (i=0;i<m->nbEBands;i++)
136    {
137       int q2;
138       celt_word16_t offset = (error[i]+QCONST16(.5f,8))*frac[i];
139       /* FIXME: Instead of giving up without warning, we should degrade everything gracefully */
140       if (ec_enc_tell(enc, 0) - bits +EC_ILOG(frac[i])> budget)
141          break;
142 #ifdef FIXED_POINT
143       /* Has to be without rounding */
144       q2 = offset>>8;
145 #else
146       q2 = (int)floor(offset);
147 #endif
148       if (q2 > frac[i]-1)
149          q2 = frac[i]-1;
150       ec_enc_uint(enc, q2, frac[i]);
151       offset = DIV32_16(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8);
152       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
153       /*printf ("%f ", error[i] - offset);*/
154    }
155    for (i=0;i<m->nbEBands;i++)
156    {
157       eBands[i] = dB2Amp(oldEBands[i]);
158    }
159    /*printf ("%d\n", ec_enc_tell(enc, 0)-9);*/
160
161    /*printf ("\n");*/
162    RESTORE_STACK;
163 }
164
165 static void unquant_energy_mono(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
166 {
167    int i;
168    int bits;
169    celt_word16_t prev = 0;
170    celt_word16_t coef = m->ePredCoef;
171    /* The .7 is a heuristic */
172    celt_word16_t beta = MULT16_16_Q15(QCONST16(.7f,15),coef);
173    bits = ec_dec_tell(dec, 0);
174    for (i=0;i<m->nbEBands;i++)
175    {
176       int qi;
177       celt_word16_t q;
178       celt_word16_t mean = MULT16_16_Q15(Q15ONE-coef,eMeans[i]);
179       /* If we didn't have enough bits to encode all the energy, just assume something safe. */
180       if (ec_dec_tell(dec, 0) - bits > budget)
181          qi = -1;
182       else
183          qi = ec_laplace_decode(dec, 6000-i*200);
184       q = qi*base_resolution;
185       
186       /*printf("%d ", qi);*/
187       /*printf("%f %f ", pred+prev+q, x);*/
188       /*printf("%f ", x-pred);*/
189       
190       oldEBands[i] = mean+MULT16_16_Q15(coef,oldEBands[i])+prev+q;
191       
192       prev = mean+prev+MULT16_16_Q15(Q15ONE-beta,q);
193    }
194    for (i=0;i<m->nbEBands;i++)
195    {
196       int q2;
197       celt_word16_t offset;
198       if (ec_dec_tell(dec, 0) - bits +EC_ILOG(frac[i])> budget)
199          break;
200       q2 = ec_dec_uint(dec, frac[i]);
201       offset = DIV32_16(SHL16(q2,8)+QCONST16(.5,8),frac[i])-QCONST16(.5f,8);
202       oldEBands[i] += PSHR32(MULT16_16(DB_SCALING*6,offset),8);
203    }
204    for (i=0;i<m->nbEBands;i++)
205    {
206       eBands[i] = dB2Amp(oldEBands[i]);
207    }
208    /*printf ("\n");*/
209 }
210
211
212
213 void quant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_enc *enc)
214 {
215    int C;
216    SAVE_STACK;
217    
218    C = m->nbChannels;
219
220    if (C==1)
221       quant_energy_mono(m, eBands, oldEBands, budget, enc);
222    else 
223 #if 1
224    {
225       int c;
226       VARDECL(celt_ener_t *E);
227       ALLOC(E, m->nbEBands, celt_ener_t);
228       for (c=0;c<C;c++)
229       {
230          int i;
231          for (i=0;i<m->nbEBands;i++)
232             E[i] = eBands[C*i+c];
233          quant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, enc);
234          for (i=0;i<m->nbEBands;i++)
235             eBands[C*i+c] = E[i];
236       }
237    }
238 #else
239       if (C==2)
240    {
241       int i;
242       int NB = m->nbEBands;
243       celt_ener_t mid[NB];
244       celt_ener_t side[NB];
245       for (i=0;i<NB;i++)
246       {
247          //left = eBands[C*i];
248          //right = eBands[C*i+1];
249          mid[i] = ENER_SCALING_1*sqrt(eBands[C*i]*eBands[C*i] + eBands[C*i+1]*eBands[C*i+1]);
250          side[i] = 20*log10((ENER_SCALING_1*eBands[2*i]+.3)/(ENER_SCALING_1*eBands[2*i+1]+.3));
251          //printf ("%f %f ", mid[i], side[i]);
252       }
253       //printf ("\n");
254       quant_energy_mono(m, mid, oldEBands, enc);
255       for (i=0;i<NB;i++)
256          side[i] = pow(10.f,floor(.5f+side[i])/10.f);
257          
258       //quant_energy_side(m, side, oldEBands+NB, enc);
259       for (i=0;i<NB;i++)
260       {
261          eBands[C*i] = ENER_SCALING*mid[i]*sqrt(side[i]/(1.f+side[i]));
262          eBands[C*i+1] = ENER_SCALING*mid[i]*sqrt(1.f/(1.f+side[i]));
263          //printf ("%f %f ", mid[i], side[i]);
264       }
265
266    } else {
267       celt_fatal("more than 2 channels not supported");
268    }
269 #endif
270    RESTORE_STACK;
271 }
272
273
274
275 void unquant_energy(const CELTMode *m, celt_ener_t *eBands, celt_word16_t *oldEBands, int budget, ec_dec *dec)
276 {
277    int C;   
278    SAVE_STACK;
279    C = m->nbChannels;
280
281    if (C==1)
282       unquant_energy_mono(m, eBands, oldEBands, budget, dec);
283    else {
284       int c;
285       VARDECL(celt_ener_t *E);
286       ALLOC(E, m->nbEBands, celt_ener_t);
287       for (c=0;c<C;c++)
288       {
289          int i;
290          unquant_energy_mono(m, E, oldEBands+c*m->nbEBands, budget/C, dec);
291          for (i=0;i<m->nbEBands;i++)
292             eBands[C*i+c] = E[i];
293       }
294    }
295    RESTORE_STACK;
296 }