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