bands.c 11.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
/* (C) 2007 Jean-Marc Valin, CSIRO
*/
/*
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
   
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
   
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
   
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
   
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

32
33
34
35
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

36
37
#include <math.h>
#include "bands.h"
38
#include "modes.h"
39
#include "vq.h"
Jean-Marc Valin's avatar
Jean-Marc Valin committed
40
#include "cwrs.h"
41

42
43

void exp_rotation(celt_norm_t *X, int len, float theta, int dir, int stride, int iter)
44
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
45
   int i, k;
46
47
48
   celt_word16_t c, s;
   c = Q15ONE*cos(theta);
   s = dir*Q15ONE*sin(theta);
49
   for (k=0;k<iter;k++)
50
   {
51
52
      /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy, 
         but at this point, I really don't think it's necessary */
53
      for (i=0;i<len-stride;i++)
54
      {
55
         celt_norm_t x1, x2;
56
57
         x1 = X[i];
         x2 = X[i+stride];
58
59
         X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
         X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
60
      }
61
      for (i=len-2*stride-1;i>=0;i--)
62
      {
63
         celt_norm_t x1, x2;
64
65
         x1 = X[i];
         x2 = X[i+stride];
66
67
         X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
         X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(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
78
79
   B = m->nbMdctBlocks;
   C = m->nbChannels;
   for (c=0;c<C;c++)
80
   {
81
82
83
84
85
      for (i=0;i<m->nbEBands;i++)
      {
         int j;
         float sum = 1e-10;
         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
100
101
   B = m->nbMdctBlocks;
   C = m->nbChannels;
   for (c=0;c<C;c++)
102
   {
103
104
105
      for (i=0;i<m->nbEBands;i++)
      {
         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
113
114
      X[i] = 0;
}

115
116
117
118
119
120
121
122
123
124
125
126
127
128
#ifdef FIXED_POINT
void renormalise_bands(const CELTMode *m, celt_norm_t *X)
{
   int i;
   VARDECL(celt_ener_t *tmpE);
   VARDECL(celt_sig_t *freq);
   ALLOC(tmpE, m->nbEBands*m->nbChannels, celt_ener_t);
   ALLOC(freq, m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1], celt_sig_t);
   for (i=0;i<m->nbMdctBlocks*m->nbChannels*m->eBands[m->nbEBands+1];i++)
      freq[i] = SHL32(EXTEND32(X[i]), 10);
   compute_band_energies(m, freq, tmpE);
   normalise_bands(m, freq, X, tmpE);
}
#else
129
void renormalise_bands(const CELTMode *m, celt_norm_t *X)
130
{
131
132
   VARDECL(celt_ener_t *tmpE);
   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;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
142
   const int *eBands = m->eBands;
143
144
145
   B = m->nbMdctBlocks;
   C = m->nbChannels;
   for (c=0;c<C;c++)
146
   {
147
148
149
      for (i=0;i<m->nbEBands;i++)
      {
         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
}


Jean-Marc Valin's avatar
Jean-Marc Valin committed
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)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
162
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
163
164
165
   int i, B;
   const int *eBands = m->eBands;
   const int *pBands = m->pBands;
166
   B = m->nbMdctBlocks*m->nbChannels;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
167
   
Jean-Marc Valin's avatar
Jean-Marc Valin committed
168
   for (i=0;i<m->nbPBands;i++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
169
   {
170
      celt_word32_t Sxy=0, Sxx=0;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
171
      int j;
172
      /* We know we're not going to overflow because Sxx can't be more than 1 (Q28) */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
173
      for (j=B*pBands[i];j<B*pBands[i+1];j++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
174
      {
175
176
         Sxy = MAC16_16(Sxy, X[j], P[j]);
         Sxx = MAC16_16(Sxx, X[j], X[j]);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
177
      }
178
179
180
181
182
183
184
185
186
187
188
      /* No negative gain allowed */
      if (Sxy < 0)
         Sxy = 0;
      /* Not sure how that would happen, just making sure */
      if (Sxy > Sxx)
         Sxy = Sxx;
      /* We need to be a bit conservative (multiply gain by 0.9), otherwise the
         residual doesn't quantise well */
      Sxy = MULT16_32_Q15(QCONST16(.9f, 15), Sxy);
      /* gain = Sxy/Sxx */
      gains[i] = DIV32_16(Sxy,ADD32(SHR32(Sxx, PGAIN_SHIFT),EPSILON));
189
      /*printf ("%f ", 1-sqrt(1-gain*gain));*/
Jean-Marc Valin's avatar
Jean-Marc Valin committed
190
   }
191
192
193
194
195
196
   /*if(rand()%10==0)
   {
      for (i=0;i<m->nbPBands;i++)
         printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
      printf ("\n");
   }*/
Jean-Marc Valin's avatar
Jean-Marc Valin committed
197
   for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
198
199
200
201
      P[i] = 0;
}

/* Apply the (quantised) gain to each "pitch band" */
202
void pitch_quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_pgain_t *gains)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
203
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
204
205
   int i, B;
   const int *pBands = m->pBands;
206
   B = m->nbMdctBlocks*m->nbChannels;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
207
   for (i=0;i<m->nbPBands;i++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
208
209
   {
      int j;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
210
      for (j=B*pBands[i];j<B*pBands[i+1];j++)
211
         P[j] = MULT16_16_Q15(gains[i], P[j]);
212
      /*printf ("%f ", gain);*/
Jean-Marc Valin's avatar
Jean-Marc Valin committed
213
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
214
   for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
215
216
217
      P[i] = 0;
}

218

219
/* Quantisation of the residual */
220
void quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
221
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
222
   int i, j, B, bits;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
223
   const int *eBands = m->eBands;
224
   float alpha = .7;
225
   VARDECL(celt_norm_t *norm);
226
227
228
229
230
   VARDECL(int *pulses);
   VARDECL(int *offsets);
   
   B = m->nbMdctBlocks*m->nbChannels;
   
231
   ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
232
233
   ALLOC(pulses, m->nbEBands, int);
   ALLOC(offsets, m->nbEBands, int);
234

Jean-Marc Valin's avatar
Jean-Marc Valin committed
235
236
   for (i=0;i<m->nbEBands;i++)
      offsets[i] = 0;
237
   /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
238
   bits = total_bits - ec_enc_tell(enc, 0) - 1;
239
   compute_allocation(m, offsets, bits, pulses);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
240
241
242
243
244
   
   /*printf("bits left: %d\n", bits);
   for (i=0;i<m->nbEBands;i++)
      printf ("%d ", pulses[i]);
   printf ("\n");*/
245
   /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
Jean-Marc Valin's avatar
Jean-Marc Valin committed
246
   for (i=0;i<m->nbEBands;i++)
247
   {
248
      int q;
249
      float theta, n;
250
      q = pulses[i];
251
252
253
      /*Scale factor of .0625f is just there to prevent overflows in fixed-point
       (has no effect on float)*/
      n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
Jean-Marc Valin's avatar
Jean-Marc Valin committed
254
      theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
255
256

      /* If pitch isn't available, use intra-frame prediction */
257
258
259
      if (eBands[i] >= m->pitchEnd || q<=0)
      {
         q -= 1;
260
         alpha = 0;
261
         if (q<0)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
262
            intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
263
264
         else
            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);
265
266
      } else {
         alpha = .7;
267
268
269
      }
      
      if (q > 0)
270
      {
271
272
         exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
         exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
273
         alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, enc);
274
         exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
275
      }
276
277
      for (j=B*eBands[i];j<B*eBands[i+1];j++)
         norm[j] = X[j] * n;
278
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
279
   for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
280
281
282
      X[i] = 0;
}

283
/* Decoding of the residual */
284
void unquant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, int total_bits, ec_dec *dec)
285
{
Jean-Marc Valin's avatar
Jean-Marc Valin committed
286
   int i, j, B, bits;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
287
   const int *eBands = m->eBands;
288
   float alpha = .7;
289
   VARDECL(celt_norm_t *norm);
290
291
292
293
294
   VARDECL(int *pulses);
   VARDECL(int *offsets);
   
   B = m->nbMdctBlocks*m->nbChannels;
   
295
   ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
296
297
   ALLOC(pulses, m->nbEBands, int);
   ALLOC(offsets, m->nbEBands, int);
298

Jean-Marc Valin's avatar
Jean-Marc Valin committed
299
300
   for (i=0;i<m->nbEBands;i++)
      offsets[i] = 0;
301
   /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
302
   bits = total_bits - ec_dec_tell(dec, 0) - 1;
303
   compute_allocation(m, offsets, bits, pulses);
Jean-Marc Valin's avatar
Jean-Marc Valin committed
304

Jean-Marc Valin's avatar
Jean-Marc Valin committed
305
   for (i=0;i<m->nbEBands;i++)
306
   {
307
      int q;
308
      float theta, n;
309
      q = pulses[i];
310
311
312
      /*Scale factor of .0625f is just there to prevent overflows in fixed-point
      (has no effect on float)*/
      n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
Jean-Marc Valin's avatar
Jean-Marc Valin committed
313
      theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
314

315
      /* If pitch isn't available, use intra-frame prediction */
316
317
318
      if (eBands[i] >= m->pitchEnd || q<=0)
      {
         q -= 1;
319
         alpha = 0;
320
         if (q<0)
Jean-Marc Valin's avatar
Jean-Marc Valin committed
321
            intra_fold(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), norm, P+B*eBands[i], B, eBands[i], eBands[m->nbEBands+1]);
322
323
         else
            intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
324
325
      } else {
         alpha = .7;
326
327
328
      }
      
      if (q > 0)
329
      {
330
         exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
331
         alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, dec);
332
         exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
333
      }
334
335
      for (j=B*eBands[i];j<B*eBands[i+1];j++)
         norm[j] = X[j] * n;
336
   }
Jean-Marc Valin's avatar
Jean-Marc Valin committed
337
   for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
338
339
      X[i] = 0;
}
340

341
void stereo_mix(const CELTMode *m, celt_norm_t *X, celt_ener_t *bank, int dir)
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
{
   int i, B, C;
   const int *eBands = m->eBands;
   B = m->nbMdctBlocks;
   C = m->nbChannels;
   for (i=0;i<m->nbEBands;i++)
   {
      int j;
      float left, right;
      float a1, a2;
      left = bank[i*C];
      right = bank[i*C+1];
      a1 = left/sqrt(.01+left*left+right*right);
      a2 = dir*right/sqrt(.01+left*left+right*right);
      for (j=B*eBands[i];j<B*eBands[i+1];j++)
      {
         float r, l;
         l = X[j*C];
         r = X[j*C+1];         
         X[j*C] = a1*l + a2*r;
         X[j*C+1] = a1*r - a2*l;
      }
   }
   for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
      X[i] = 0;

368
}