Using only a sign for the intra prediction instead of using up to 5 bits for
[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 "stack_alloc.h"
42 #include "os_support.h"
43 #include "mathops.h"
44
45 void exp_rotation(celt_norm_t *X, int len, int dir, int stride, int iter)
46 {
47    int i, k;
48    celt_word16_t c, s;
49    celt_norm_t *Xptr;
50    /* Equivalent to cos(.3) and sin(.3) */
51    c = QCONST16(0.95534,15);
52    s = dir*QCONST16(0.29552,15);
53    for (k=0;k<iter;k++)
54    {
55       /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy, 
56          but at this point, I really don't think it's necessary */
57       Xptr = X;
58       for (i=0;i<len-stride;i++)
59       {
60          celt_norm_t x1, x2;
61          x1 = Xptr[0];
62          x2 = Xptr[stride];
63          Xptr[stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
64          *Xptr++      = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
65       }
66       Xptr = &X[len-2*stride-1];
67       for (i=len-2*stride-1;i>=0;i--)
68       {
69          celt_norm_t x1, x2;
70          x1 = Xptr[0];
71          x2 = Xptr[stride];
72          Xptr[stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
73          *Xptr--      = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
74       }
75    }
76 }
77
78
79
80 const celt_word16_t sqrtC_1[2] = {QCONST16(1.f, 14), QCONST16(1.414214f, 14)};
81
82 #ifdef FIXED_POINT
83 /* Compute the amplitude (sqrt energy) in each of the bands */
84 void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank)
85 {
86    int i, c;
87    const celt_int16_t *eBands = m->eBands;
88    const int C = CHANNELS(m);
89    for (c=0;c<C;c++)
90    {
91       for (i=0;i<m->nbEBands;i++)
92       {
93          int j;
94          celt_word32_t maxval=0;
95          celt_word32_t sum = 0;
96          
97          j=eBands[i]; do {
98             maxval = MAX32(maxval, X[j*C+c]);
99             maxval = MAX32(maxval, -X[j*C+c]);
100          } while (++j<eBands[i+1]);
101          
102          if (maxval > 0)
103          {
104             int shift = celt_ilog2(maxval)-10;
105             j=eBands[i]; do {
106                sum += MULT16_16(EXTRACT16(VSHR32(X[j*C+c],shift)),
107                                 EXTRACT16(VSHR32(X[j*C+c],shift)));
108             } while (++j<eBands[i+1]);
109             /* We're adding one here to make damn sure we never end up with a pitch vector that's
110                larger than unity norm */
111             bank[i*C+c] = EPSILON+VSHR32(EXTEND32(celt_sqrt(sum)),-shift);
112          } else {
113             bank[i*C+c] = EPSILON;
114          }
115          /*printf ("%f ", bank[i*C+c]);*/
116       }
117    }
118    /*printf ("\n");*/
119 }
120
121 /* Normalise each band such that the energy is one. */
122 void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank)
123 {
124    int i, c;
125    const celt_int16_t *eBands = m->eBands;
126    const int C = CHANNELS(m);
127    for (c=0;c<C;c++)
128    {
129       i=0; do {
130          celt_word16_t g;
131          int j,shift;
132          celt_word16_t E;
133          shift = celt_zlog2(bank[i*C+c])-13;
134          E = VSHR32(bank[i*C+c], shift);
135          g = EXTRACT16(celt_rcp(SHR32(MULT16_16(E,sqrtC_1[C-1]),11)));
136          j=eBands[i]; do {
137             X[j*C+c] = MULT16_16_Q15(VSHR32(freq[j*C+c],shift-1),g);
138          } while (++j<eBands[i+1]);
139       } while (++i<m->nbEBands);
140    }
141 }
142
143 #ifndef DISABLE_STEREO
144 void renormalise_bands(const CELTMode *m, celt_norm_t * restrict X)
145 {
146    int i;
147    VARDECL(celt_ener_t, tmpE);
148    VARDECL(celt_sig_t, freq);
149    SAVE_STACK;
150    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
151    ALLOC(freq, m->nbChannels*m->eBands[m->nbEBands+1], celt_sig_t);
152    for (i=0;i<m->nbChannels*m->eBands[m->nbEBands+1];i++)
153       freq[i] = SHL32(EXTEND32(X[i]), 10);
154    compute_band_energies(m, freq, tmpE);
155    normalise_bands(m, freq, X, tmpE);
156    RESTORE_STACK;
157 }
158 #endif /* DISABLE_STEREO */
159 #else /* FIXED_POINT */
160 /* Compute the amplitude (sqrt energy) in each of the bands */
161 void compute_band_energies(const CELTMode *m, const celt_sig_t *X, celt_ener_t *bank)
162 {
163    int i, c;
164    const celt_int16_t *eBands = m->eBands;
165    const int C = CHANNELS(m);
166    for (c=0;c<C;c++)
167    {
168       for (i=0;i<m->nbEBands;i++)
169       {
170          int j;
171          celt_word32_t sum = 1e-10;
172          for (j=eBands[i];j<eBands[i+1];j++)
173             sum += X[j*C+c]*X[j*C+c];
174          bank[i*C+c] = sqrt(sum);
175          /*printf ("%f ", bank[i*C+c]);*/
176       }
177    }
178    /*printf ("\n");*/
179 }
180
181 /* Normalise each band such that the energy is one. */
182 void normalise_bands(const CELTMode *m, const celt_sig_t * restrict freq, celt_norm_t * restrict X, const celt_ener_t *bank)
183 {
184    int i, c;
185    const celt_int16_t *eBands = m->eBands;
186    const int C = CHANNELS(m);
187    for (c=0;c<C;c++)
188    {
189       for (i=0;i<m->nbEBands;i++)
190       {
191          int j;
192          celt_word16_t g = 1.f/(1e-10+bank[i*C+c]*sqrt(C));
193          for (j=eBands[i];j<eBands[i+1];j++)
194             X[j*C+c] = freq[j*C+c]*g;
195       }
196    }
197 }
198
199 #ifndef DISABLE_STEREO
200 void renormalise_bands(const CELTMode *m, celt_norm_t * restrict X)
201 {
202    VARDECL(celt_ener_t, tmpE);
203    SAVE_STACK;
204    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
205    compute_band_energies(m, X, tmpE);
206    normalise_bands(m, X, X, tmpE);
207    RESTORE_STACK;
208 }
209 #endif /* DISABLE_STEREO */
210 #endif /* FIXED_POINT */
211
212 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
213 void denormalise_bands(const CELTMode *m, const celt_norm_t * restrict X, celt_sig_t * restrict freq, const celt_ener_t *bank)
214 {
215    int i, c;
216    const celt_int16_t *eBands = m->eBands;
217    const int C = CHANNELS(m);
218    if (C>2)
219       celt_fatal("denormalise_bands() not implemented for >2 channels");
220    for (c=0;c<C;c++)
221    {
222       for (i=0;i<m->nbEBands;i++)
223       {
224          int j;
225          celt_word32_t g = MULT16_32_Q13(sqrtC_1[C-1],bank[i*C+c]);
226          j=eBands[i]; do {
227             freq[j*C+c] = MULT16_32_Q15(X[j*C+c], g);
228          } while (++j<eBands[i+1]);
229       }
230    }
231    for (i=C*eBands[m->nbEBands];i<C*eBands[m->nbEBands+1];i++)
232       freq[i] = 0;
233 }
234
235
236 /* Compute the best gain for each "pitch band" */
237 void compute_pitch_gain(const CELTMode *m, const celt_norm_t *X, const celt_norm_t *P, celt_pgain_t *gains)
238 {
239    int i;
240    const celt_int16_t *pBands = m->pBands;
241    const int C = CHANNELS(m);
242
243    for (i=0;i<m->nbPBands;i++)
244    {
245       celt_word32_t Sxy=0, Sxx=0;
246       int j;
247       /* We know we're not going to overflow because Sxx can't be more than 1 (Q28) */
248       for (j=C*pBands[i];j<C*pBands[i+1];j++)
249       {
250          Sxy = MAC16_16(Sxy, X[j], P[j]);
251          Sxx = MAC16_16(Sxx, X[j], X[j]);
252       }
253       /* No negative gain allowed */
254       if (Sxy < 0)
255          Sxy = 0;
256       /* Not sure how that would happen, just making sure */
257       if (Sxy > Sxx)
258          Sxy = Sxx;
259       /* We need to be a bit conservative (multiply gain by 0.9), otherwise the
260          residual doesn't quantise well */
261       Sxy = MULT16_32_Q15(QCONST16(.9f, 15), Sxy);
262       /* gain = Sxy/Sxx */
263       gains[i] = EXTRACT16(celt_div(Sxy,ADD32(SHR32(Sxx, PGAIN_SHIFT),EPSILON)));
264       /*printf ("%f ", 1-sqrt(1-gain*gain));*/
265    }
266    /*if(rand()%10==0)
267    {
268       for (i=0;i<m->nbPBands;i++)
269          printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
270       printf ("\n");
271    }*/
272 }
273
274 /* Apply the (quantised) gain to each "pitch band" */
275 void pitch_quant_bands(const CELTMode *m, celt_norm_t * restrict P, const celt_pgain_t * restrict gains)
276 {
277    int i;
278    const celt_int16_t *pBands = m->pBands;
279    const int C = CHANNELS(m);
280    for (i=0;i<m->nbPBands;i++)
281    {
282       int j;
283       for (j=C*pBands[i];j<C*pBands[i+1];j++)
284          P[j] = MULT16_16_Q15(gains[i], P[j]);
285       /*printf ("%f ", gain);*/
286    }
287    for (i=C*pBands[m->nbPBands];i<C*pBands[m->nbPBands+1];i++)
288       P[i] = 0;
289 }
290
291
292 /* Quantisation of the residual */
293 void quant_bands(const CELTMode *m, celt_norm_t * restrict X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
294 {
295    int i, j, bits;
296    const celt_int16_t * restrict eBands = m->eBands;
297    celt_norm_t * restrict norm;
298    VARDECL(celt_norm_t, _norm);
299    VARDECL(int, pulses);
300    VARDECL(int, offsets);
301    const int C = CHANNELS(m);
302    SAVE_STACK;
303
304    ALLOC(_norm, C*eBands[m->nbEBands+1], celt_norm_t);
305    ALLOC(pulses, m->nbEBands, int);
306    ALLOC(offsets, m->nbEBands, int);
307    norm = _norm;
308
309    for (i=0;i<m->nbEBands;i++)
310       offsets[i] = 0;
311    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
312    bits = total_bits - ec_enc_tell(enc, 0) - 1;
313    compute_allocation(m, offsets, bits, pulses);
314    
315    /*printf("bits left: %d\n", bits);
316    for (i=0;i<m->nbEBands;i++)
317       printf ("%d ", pulses[i]);
318    printf ("\n");*/
319    /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
320    for (i=0;i<m->nbEBands;i++)
321    {
322       int q;
323       celt_word16_t n;
324       q = pulses[i];
325       n = SHL16(celt_sqrt(C*(eBands[i+1]-eBands[i])),11);
326
327       /* If pitch isn't available, use intra-frame prediction */
328       if (eBands[i] >= m->pitchEnd || q<=0)
329       {
330          q -= 1;
331          if (q<0)
332             intra_fold(m, X+C*eBands[i], eBands[i+1]-eBands[i], norm, P+C*eBands[i], eBands[i], eBands[m->nbEBands+1]);
333          else
334             intra_prediction(m, X+C*eBands[i], W+C*eBands[i], eBands[i+1]-eBands[i], q, norm, P+C*eBands[i], eBands[i], eBands[m->nbEBands+1], enc);
335       }
336       
337       if (q > 0)
338       {
339          int nb_rotations = q <= 2*C ? 2*C/q : 0;
340          if (nb_rotations != 0)
341          {
342             exp_rotation(P+C*eBands[i], C*(eBands[i+1]-eBands[i]), -1, C, nb_rotations);
343             exp_rotation(X+C*eBands[i], C*(eBands[i+1]-eBands[i]), -1, C, nb_rotations);
344          }
345          alg_quant(X+C*eBands[i], W+C*eBands[i], C*(eBands[i+1]-eBands[i]), q, P+C*eBands[i], enc);
346          if (nb_rotations != 0)
347             exp_rotation(X+C*eBands[i], C*(eBands[i+1]-eBands[i]), 1, C, nb_rotations);
348       }
349       for (j=C*eBands[i];j<C*eBands[i+1];j++)
350          norm[j] = MULT16_16_Q15(n,X[j]);
351    }
352    RESTORE_STACK;
353 }
354
355 /* Decoding of the residual */
356 void unquant_bands(const CELTMode *m, celt_norm_t * restrict X, celt_norm_t *P, int total_bits, ec_dec *dec)
357 {
358    int i, j, bits;
359    const celt_int16_t * restrict eBands = m->eBands;
360    celt_norm_t * restrict norm;
361    VARDECL(celt_norm_t, _norm);
362    VARDECL(int, pulses);
363    VARDECL(int, offsets);
364    const int C = CHANNELS(m);
365    SAVE_STACK;
366
367    ALLOC(_norm, C*eBands[m->nbEBands+1], celt_norm_t);
368    ALLOC(pulses, m->nbEBands, int);
369    ALLOC(offsets, m->nbEBands, int);
370    norm = _norm;
371
372    for (i=0;i<m->nbEBands;i++)
373       offsets[i] = 0;
374    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
375    bits = total_bits - ec_dec_tell(dec, 0) - 1;
376    compute_allocation(m, offsets, bits, pulses);
377
378    for (i=0;i<m->nbEBands;i++)
379    {
380       int q;
381       celt_word16_t n;
382       q = pulses[i];
383       n = SHL16(celt_sqrt(C*(eBands[i+1]-eBands[i])),11);
384
385       /* If pitch isn't available, use intra-frame prediction */
386       if (eBands[i] >= m->pitchEnd || q<=0)
387       {
388          q -= 1;
389          if (q<0)
390             intra_fold(m, X+C*eBands[i], eBands[i+1]-eBands[i], norm, P+C*eBands[i], eBands[i], eBands[m->nbEBands+1]);
391          else
392             intra_unquant(m, X+C*eBands[i], eBands[i+1]-eBands[i], q, norm, P+C*eBands[i], eBands[i], eBands[m->nbEBands+1], dec);
393       }
394       
395       if (q > 0)
396       {
397          int nb_rotations = q <= 2*C ? 2*C/q : 0;
398          if (nb_rotations != 0)
399             exp_rotation(P+C*eBands[i], C*(eBands[i+1]-eBands[i]), -1, C, nb_rotations);
400          alg_unquant(X+C*eBands[i], C*(eBands[i+1]-eBands[i]), q, P+C*eBands[i], dec);
401          if (nb_rotations != 0)
402             exp_rotation(X+C*eBands[i], C*(eBands[i+1]-eBands[i]), 1, C, nb_rotations);
403       }
404       for (j=C*eBands[i];j<C*eBands[i+1];j++)
405          norm[j] = MULT16_16_Q15(n,X[j]);
406    }
407    RESTORE_STACK;
408 }
409
410 #ifndef DISABLE_STEREO
411 void stereo_mix(const CELTMode *m, celt_norm_t *X, const celt_ener_t *bank, int dir)
412 {
413    int i;
414    const celt_int16_t *eBands = m->eBands;
415    const int C = CHANNELS(m);
416    for (i=0;i<m->nbEBands;i++)
417    {
418       int j;
419       celt_word16_t left, right;
420       celt_word16_t a1, a2;
421       celt_word16_t norm;
422 #ifdef FIXED_POINT
423       int shift = celt_zlog2(MAX32(bank[i*C], bank[i*C+1]))-13;
424 #endif
425       left = VSHR32(bank[i*C],shift);
426       right = VSHR32(bank[i*C+1],shift);
427       norm = EPSILON + celt_sqrt(EPSILON+MULT16_16(left,left)+MULT16_16(right,right));
428       a1 = DIV32_16(SHL32(EXTEND32(left),14),norm);
429       a2 = dir*DIV32_16(SHL32(EXTEND32(right),14),norm);
430       for (j=eBands[i];j<eBands[i+1];j++)
431       {
432          celt_norm_t r, l;
433          l = X[j*C];
434          r = X[j*C+1];
435          X[j*C] = MULT16_16_Q14(a1,l) + MULT16_16_Q14(a2,r);
436          X[j*C+1] = MULT16_16_Q14(a1,r) - MULT16_16_Q14(a2,l);
437       }
438    }
439 }
440 #endif