Removed deprecated mode interface and added missing include
[opus.git] / libcelt / 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 <math.h>
37 #include "bands.h"
38 #include "modes.h"
39 #include "vq.h"
40 #include "cwrs.h"
41 #include "os_support.h"
42
43 void exp_rotation(celt_norm_t *X, int len, celt_word16_t theta, int dir, int stride, int iter)
44 {
45    int i, k;
46    celt_word16_t c, s;
47    /* c = cos(theta); s = dir*sin(theta); but we're approximating here for small theta */
48    c = Q15ONE-MULT16_16_Q15(QCONST16(.5f,15),MULT16_16_Q15(theta,theta));
49    s = dir*theta;
50    for (k=0;k<iter;k++)
51    {
52       /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy, 
53          but at this point, I really don't think it's necessary */
54       for (i=0;i<len-stride;i++)
55       {
56          celt_norm_t x1, x2;
57          x1 = X[i];
58          x2 = X[i+stride];
59          X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
60          X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
61       }
62       for (i=len-2*stride-1;i>=0;i--)
63       {
64          celt_norm_t x1, x2;
65          x1 = X[i];
66          x2 = X[i+stride];
67          X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
68          X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
69       }
70    }
71 }
72
73 /* Compute the amplitude (sqrt energy) in each of the bands */
74 void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank)
75 {
76    int i, c, B, C;
77    const int *eBands = m->eBands;
78    B = m->nbMdctBlocks;
79    C = m->nbChannels;
80    for (c=0;c<C;c++)
81    {
82       for (i=0;i<m->nbEBands;i++)
83       {
84          int j;
85          float sum = 1e-10;
86          for (j=B*eBands[i];j<B*eBands[i+1];j++)
87             sum += SIG_SCALING_1*SIG_SCALING_1*X[j*C+c]*X[j*C+c];
88          bank[i*C+c] = ENER_SCALING*sqrt(sum);
89          /*printf ("%f ", bank[i*C+c]);*/
90       }
91    }
92    /*printf ("\n");*/
93 }
94
95 /* Normalise each band such that the energy is one. */
96 void normalise_bands(const CELTMode *m, const celt_sig_t *freq, celt_norm_t *X, const celt_ener_t *bank)
97 {
98    int i, c, B, C;
99    const int *eBands = m->eBands;
100    B = m->nbMdctBlocks;
101    C = m->nbChannels;
102    for (c=0;c<C;c++)
103    {
104       for (i=0;i<m->nbEBands;i++)
105       {
106          int j;
107          float g = 1.f/(1e-10+ENER_SCALING_1*bank[i*C+c]*sqrt(C));
108          for (j=B*eBands[i];j<B*eBands[i+1];j++)
109             X[j*C+c] = NORM_SCALING*SIG_SCALING_1*freq[j*C+c]*g;
110       }
111    }
112    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
113       X[i] = 0;
114 }
115
116 #ifdef FIXED_POINT
117 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
118 {
119    int i;
120    VARDECL(celt_ener_t *tmpE);
121    VARDECL(celt_sig_t *freq);
122    SAVE_STACK;
123    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
124    ALLOC(freq, m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1], celt_sig_t);
125    for (i=0;i<m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1];i++)
126       freq[i] = SHL32(EXTEND32(X[i]), 10);
127    compute_band_energies(m, freq, tmpE);
128    normalise_bands(m, freq, X, tmpE);
129    RESTORE_STACK;
130 }
131 #else
132 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
133 {
134    VARDECL(celt_ener_t *tmpE);
135    SAVE_STACK;
136    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
137    compute_band_energies(m, X, tmpE);
138    normalise_bands(m, X, X, tmpE);
139    RESTORE_STACK;
140 }
141 #endif
142
143 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
144 void denormalise_bands(const CELTMode *m, const celt_norm_t *X, celt_sig_t *freq, const celt_ener_t *bank)
145 {
146    int i, c, B, C;
147    const celt_word16_t sqrtC_1[2] = {QCONST16(1.f, 14), QCONST16(1.414214f, 14)};
148    const int *eBands = m->eBands;
149    B = m->nbMdctBlocks;
150    C = m->nbChannels;
151    if (C>2)
152       celt_fatal("denormalise_bands() not implemented for >2 channels");
153    for (c=0;c<C;c++)
154    {
155       for (i=0;i<m->nbEBands;i++)
156       {
157          int j;
158          celt_word32_t g = MULT16_32_Q14(sqrtC_1[C-1],bank[i*C+c]);
159          for (j=B*eBands[i];j<B*eBands[i+1];j++)
160             freq[j*C+c] = MULT16_32_Q14(X[j*C+c], g);
161       }
162    }
163    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
164       freq[i] = 0;
165 }
166
167
168 /* Compute the best gain for each "pitch band" */
169 void compute_pitch_gain(const CELTMode *m, const celt_norm_t *X, const celt_norm_t *P, celt_pgain_t *gains)
170 {
171    int i, B;
172    const int *pBands = m->pBands;
173    B = m->nbMdctBlocks*m->nbChannels;
174    
175    for (i=0;i<m->nbPBands;i++)
176    {
177       celt_word32_t Sxy=0, Sxx=0;
178       int j;
179       /* We know we're not going to overflow because Sxx can't be more than 1 (Q28) */
180       for (j=B*pBands[i];j<B*pBands[i+1];j++)
181       {
182          Sxy = MAC16_16(Sxy, X[j], P[j]);
183          Sxx = MAC16_16(Sxx, X[j], X[j]);
184       }
185       /* No negative gain allowed */
186       if (Sxy < 0)
187          Sxy = 0;
188       /* Not sure how that would happen, just making sure */
189       if (Sxy > Sxx)
190          Sxy = Sxx;
191       /* We need to be a bit conservative (multiply gain by 0.9), otherwise the
192          residual doesn't quantise well */
193       Sxy = MULT16_32_Q15(QCONST16(.9f, 15), Sxy);
194       /* gain = Sxy/Sxx */
195       gains[i] = DIV32_16(Sxy,ADD32(SHR32(Sxx, PGAIN_SHIFT),EPSILON));
196       /*printf ("%f ", 1-sqrt(1-gain*gain));*/
197    }
198    /*if(rand()%10==0)
199    {
200       for (i=0;i<m->nbPBands;i++)
201          printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
202       printf ("\n");
203    }*/
204 }
205
206 /* Apply the (quantised) gain to each "pitch band" */
207 void pitch_quant_bands(const CELTMode *m, celt_norm_t *P, const celt_pgain_t *gains)
208 {
209    int i, B;
210    const int *pBands = m->pBands;
211    B = m->nbMdctBlocks*m->nbChannels;
212    for (i=0;i<m->nbPBands;i++)
213    {
214       int j;
215       for (j=B*pBands[i];j<B*pBands[i+1];j++)
216          P[j] = MULT16_16_Q15(gains[i], P[j]);
217       /*printf ("%f ", gain);*/
218    }
219    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
220       P[i] = 0;
221 }
222
223
224 /* Quantisation of the residual */
225 void quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
226 {
227    int i, j, B, bits;
228    const int *eBands = m->eBands;
229    celt_word16_t alpha;
230    VARDECL(celt_norm_t *norm);
231    VARDECL(int *pulses);
232    VARDECL(int *offsets);
233    SAVE_STACK;
234
235    B = m->nbMdctBlocks*m->nbChannels;
236    
237    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
238    ALLOC(pulses, m->nbEBands, int);
239    ALLOC(offsets, m->nbEBands, int);
240
241    for (i=0;i<m->nbEBands;i++)
242       offsets[i] = 0;
243    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
244    bits = total_bits - ec_enc_tell(enc, 0) - 1;
245    compute_allocation(m, offsets, bits, pulses);
246    
247    /*printf("bits left: %d\n", bits);
248    for (i=0;i<m->nbEBands;i++)
249       printf ("%d ", pulses[i]);
250    printf ("\n");*/
251    /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
252    for (i=0;i<m->nbEBands;i++)
253    {
254       int q;
255       celt_word16_t theta;
256       float n;
257       q = pulses[i];
258       /*Scale factor of .0625f is just there to prevent overflows in fixed-point
259        (has no effect on float)*/
260       n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
261       theta = Q15ONE*.007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
262
263       /* If pitch isn't available, use intra-frame prediction */
264       if (eBands[i] >= m->pitchEnd || q<=0)
265       {
266          q -= 1;
267          alpha = 0;
268          if (q<0)
269             intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
270          else
271             intra_prediction(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], enc);
272       } else {
273          alpha = QCONST16(.7f,15);
274       }
275       
276       if (q > 0)
277       {
278          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
279          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
280          alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, enc);
281          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
282       }
283       for (j=B*eBands[i];j<B*eBands[i+1];j++)
284          norm[j] = X[j] * n;
285    }
286    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
287       X[i] = 0;
288    RESTORE_STACK;
289 }
290
291 /* Decoding of the residual */
292 void unquant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, int total_bits, ec_dec *dec)
293 {
294    int i, j, B, bits;
295    const int *eBands = m->eBands;
296    celt_word16_t alpha;
297    VARDECL(celt_norm_t *norm);
298    VARDECL(int *pulses);
299    VARDECL(int *offsets);
300    SAVE_STACK;
301
302    B = m->nbMdctBlocks*m->nbChannels;
303    
304    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
305    ALLOC(pulses, m->nbEBands, int);
306    ALLOC(offsets, m->nbEBands, int);
307
308    for (i=0;i<m->nbEBands;i++)
309       offsets[i] = 0;
310    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
311    bits = total_bits - ec_dec_tell(dec, 0) - 1;
312    compute_allocation(m, offsets, bits, pulses);
313
314    for (i=0;i<m->nbEBands;i++)
315    {
316       int q;
317       celt_word16_t theta;
318       float n;
319       q = pulses[i];
320       /*Scale factor of .0625f is just there to prevent overflows in fixed-point
321       (has no effect on float)*/
322       n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
323       theta = Q15ONE*.007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
324
325       /* If pitch isn't available, use intra-frame prediction */
326       if (eBands[i] >= m->pitchEnd || q<=0)
327       {
328          q -= 1;
329          alpha = 0;
330          if (q<0)
331             intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
332          else
333             intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
334       } else {
335          alpha = QCONST16(.7f,15);
336       }
337       
338       if (q > 0)
339       {
340          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
341          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, dec);
342          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
343       }
344       for (j=B*eBands[i];j<B*eBands[i+1];j++)
345          norm[j] = X[j] * n;
346    }
347    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
348       X[i] = 0;
349    RESTORE_STACK;
350 }
351
352 void stereo_mix(const CELTMode *m, celt_norm_t *X, const celt_ener_t *bank, int dir)
353 {
354    int i, B, C;
355    const int *eBands = m->eBands;
356    B = m->nbMdctBlocks;
357    C = m->nbChannels;
358    for (i=0;i<m->nbEBands;i++)
359    {
360       int j;
361       float left, right;
362       float a1, a2;
363       left = bank[i*C];
364       right = bank[i*C+1];
365       a1 = left/sqrt(.01+left*left+right*right);
366       a2 = dir*right/sqrt(.01+left*left+right*right);
367       for (j=B*eBands[i];j<B*eBands[i+1];j++)
368       {
369          float r, l;
370          l = X[j*C];
371          r = X[j*C+1];         
372          X[j*C] = a1*l + a2*r;
373          X[j*C+1] = a1*r - a2*l;
374       }
375    }
376    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
377       X[i] = 0;
378
379 }