fixed-point: got stereo to work again by fixing renormalise_bands()
[opus.git] / libcelt / bands.c
1 /* (C) 2007 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
42 /** Applies a series of rotations so that pulses are spread like a two-sided
43 exponential. The effect of this is to reduce the tonal noise created by the
44 sparse spectrum resulting from the pulse codebook */
45 static void exp_rotation(celt_norm_t *X, int len, float theta, int dir, int stride, int iter)
46 {
47    int i, k;
48    float c, s;
49    c = cos(theta);
50    s = dir*sin(theta);
51    for (k=0;k<iter;k++)
52    {
53       for (i=0;i<len-stride;i++)
54       {
55          float x1, x2;
56          x1 = X[i];
57          x2 = X[i+stride];
58          X[i] = c*x1 - s*x2;
59          X[i+stride] = c*x2 + s*x1;
60       }
61       for (i=len-2*stride-1;i>=0;i--)
62       {
63          float x1, x2;
64          x1 = X[i];
65          x2 = X[i+stride];
66          X[i] = c*x1 - s*x2;
67          X[i+stride] = c*x2 + s*x1;
68       }
69    }
70 }
71
72 /* Compute the amplitude (sqrt energy) in each of the bands */
73 void compute_band_energies(const CELTMode *m, celt_sig_t *X, celt_ener_t *bank)
74 {
75    int i, c, B, C;
76    const int *eBands = m->eBands;
77    B = m->nbMdctBlocks;
78    C = m->nbChannels;
79    for (c=0;c<C;c++)
80    {
81       for (i=0;i<m->nbEBands;i++)
82       {
83          int j;
84          float sum = 1e-10;
85          for (j=B*eBands[i];j<B*eBands[i+1];j++)
86             sum += SIG_SCALING_1*SIG_SCALING_1*X[j*C+c]*X[j*C+c];
87          bank[i*C+c] = ENER_SCALING*sqrt(sum);
88          /*printf ("%f ", bank[i*C+c]);*/
89       }
90    }
91    /*printf ("\n");*/
92 }
93
94 /* Normalise each band such that the energy is one. */
95 void normalise_bands(const CELTMode *m, celt_sig_t *freq, celt_norm_t *X, celt_ener_t *bank)
96 {
97    int i, c, B, C;
98    const int *eBands = m->eBands;
99    B = m->nbMdctBlocks;
100    C = m->nbChannels;
101    for (c=0;c<C;c++)
102    {
103       for (i=0;i<m->nbEBands;i++)
104       {
105          int j;
106          float g = 1.f/(1e-10+ENER_SCALING_1*bank[i*C+c]*sqrt(C));
107          for (j=B*eBands[i];j<B*eBands[i+1];j++)
108             X[j*C+c] = NORM_SCALING*SIG_SCALING_1*freq[j*C+c]*g;
109       }
110    }
111    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
112       X[i] = 0;
113 }
114
115 #ifdef FIXED_POINT
116 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
117 {
118    int i;
119    VARDECL(celt_ener_t *tmpE);
120    VARDECL(celt_sig_t *freq);
121    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
122    ALLOC(freq, m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1], celt_sig_t);
123    for (i=0;i<m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1];i++)
124       freq[i] = SHL32(EXTEND32(X[i]), 10);
125    compute_band_energies(m, freq, tmpE);
126    normalise_bands(m, freq, X, tmpE);
127 }
128 #else
129 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
130 {
131    VARDECL(celt_ener_t *tmpE);
132    ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
133    compute_band_energies(m, X, tmpE);
134    normalise_bands(m, X, X, tmpE);
135 }
136 #endif
137
138 /* De-normalise the energy to produce the synthesis from the unit-energy bands */
139 void denormalise_bands(const CELTMode *m, celt_norm_t *X, celt_sig_t *freq, celt_ener_t *bank)
140 {
141    int i, c, B, C;
142    const int *eBands = m->eBands;
143    B = m->nbMdctBlocks;
144    C = m->nbChannels;
145    for (c=0;c<C;c++)
146    {
147       for (i=0;i<m->nbEBands;i++)
148       {
149          int j;
150          float g = ENER_SCALING_1*sqrt(C)*bank[i*C+c];
151          for (j=B*eBands[i];j<B*eBands[i+1];j++)
152             freq[j*C+c] = NORM_SCALING_1*SIG_SCALING*X[j*C+c] * g;
153       }
154    }
155    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
156       freq[i] = 0;
157 }
158
159
160 /* Compute the best gain for each "pitch band" */
161 void compute_pitch_gain(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_pgain_t *gains, celt_ener_t *bank)
162 {
163    int i, B;
164    const int *eBands = m->eBands;
165    const int *pBands = m->pBands;
166    VARDECL(float *w);
167    B = m->nbMdctBlocks*m->nbChannels;
168    ALLOC(w, B*eBands[m->nbEBands], float);
169    for (i=0;i<m->nbEBands;i++)
170    {
171       int j;
172       for (j=B*eBands[i];j<B*eBands[i+1];j++)
173          w[j] = bank[i]*ENER_SCALING_1;
174    }
175
176    
177    for (i=0;i<m->nbPBands;i++)
178    {
179       float Sxy=0;
180       float Sxx = 0;
181       int j;
182       float gain;
183       for (j=B*pBands[i];j<B*pBands[i+1];j++)
184       {
185          Sxy += 1.f*X[j]*P[j]*w[j];
186          Sxx += 1.f*X[j]*X[j]*w[j];
187       }
188       gain = Sxy/(1e-10*NORM_SCALING*NORM_SCALING+Sxx);
189       if (gain > 1.f)
190          gain = 1.f;
191       if (gain < 0.0f)
192          gain = 0.0f;
193       /* We need to be a bit conservative, otherwise residual doesn't quantise well */
194       gain *= .9f;
195       gains[i] = PGAIN_SCALING*gain;
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    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
205       P[i] = 0;
206 }
207
208 /* Apply the (quantised) gain to each "pitch band" */
209 void pitch_quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_pgain_t *gains)
210 {
211    int i, B;
212    const int *pBands = m->pBands;
213    B = m->nbMdctBlocks*m->nbChannels;
214    for (i=0;i<m->nbPBands;i++)
215    {
216       int j;
217       for (j=B*pBands[i];j<B*pBands[i+1];j++)
218          P[j] *= PGAIN_SCALING_1*gains[i];
219       /*printf ("%f ", gain);*/
220    }
221    for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
222       P[i] = 0;
223 }
224
225
226 /* Quantisation of the residual */
227 void quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
228 {
229    int i, j, B, bits;
230    const int *eBands = m->eBands;
231    float alpha = .7;
232    VARDECL(celt_norm_t *norm);
233    VARDECL(int *pulses);
234    VARDECL(int *offsets);
235    
236    B = m->nbMdctBlocks*m->nbChannels;
237    
238    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
239    ALLOC(pulses, m->nbEBands, int);
240    ALLOC(offsets, m->nbEBands, int);
241
242    for (i=0;i<m->nbEBands;i++)
243       offsets[i] = 0;
244    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
245    bits = total_bits - ec_enc_tell(enc, 0) - 1;
246    compute_allocation(m, offsets, bits, pulses);
247    
248    /*printf("bits left: %d\n", bits);
249    for (i=0;i<m->nbEBands;i++)
250       printf ("%d ", pulses[i]);
251    printf ("\n");*/
252    /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
253    for (i=0;i<m->nbEBands;i++)
254    {
255       int q;
256       float theta, 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 = .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 = .7;
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 }
289
290 /* Decoding of the residual */
291 void unquant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, int total_bits, ec_dec *dec)
292 {
293    int i, j, B, bits;
294    const int *eBands = m->eBands;
295    float alpha = .7;
296    VARDECL(celt_norm_t *norm);
297    VARDECL(int *pulses);
298    VARDECL(int *offsets);
299    
300    B = m->nbMdctBlocks*m->nbChannels;
301    
302    ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
303    ALLOC(pulses, m->nbEBands, int);
304    ALLOC(offsets, m->nbEBands, int);
305
306    for (i=0;i<m->nbEBands;i++)
307       offsets[i] = 0;
308    /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
309    bits = total_bits - ec_dec_tell(dec, 0) - 1;
310    compute_allocation(m, offsets, bits, pulses);
311
312    for (i=0;i<m->nbEBands;i++)
313    {
314       int q;
315       float theta, n;
316       q = pulses[i];
317       /*Scale factor of .0625f is just there to prevent overflows in fixed-point
318       (has no effect on float)*/
319       n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
320       theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
321
322       /* If pitch isn't available, use intra-frame prediction */
323       if (eBands[i] >= m->pitchEnd || q<=0)
324       {
325          q -= 1;
326          alpha = 0;
327          if (q<0)
328             intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
329          else
330             intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
331       } else {
332          alpha = .7;
333       }
334       
335       if (q > 0)
336       {
337          exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
338          alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, dec);
339          exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
340       }
341       for (j=B*eBands[i];j<B*eBands[i+1];j++)
342          norm[j] = X[j] * n;
343    }
344    for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
345       X[i] = 0;
346 }
347
348 void stereo_mix(const CELTMode *m, celt_norm_t *X, celt_ener_t *bank, int dir)
349 {
350    int i, B, C;
351    const int *eBands = m->eBands;
352    B = m->nbMdctBlocks;
353    C = m->nbChannels;
354    for (i=0;i<m->nbEBands;i++)
355    {
356       int j;
357       float left, right;
358       float a1, a2;
359       left = bank[i*C];
360       right = bank[i*C+1];
361       a1 = left/sqrt(.01+left*left+right*right);
362       a2 = dir*right/sqrt(.01+left*left+right*right);
363       for (j=B*eBands[i];j<B*eBands[i+1];j++)
364       {
365          float r, l;
366          l = X[j*C];
367          r = X[j*C+1];         
368          X[j*C] = a1*l + a2*r;
369          X[j*C+1] = a1*r - a2*l;
370       }
371    }
372    for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
373       X[i] = 0;
374
375 }