1 /* (C) 2007 Jean-Marc Valin, CSIRO
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions
8 - Redistributions of source code must retain the above copyright
9 notice, this list of conditions and the following disclaimer.
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.
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.
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.
43 void exp_rotation(celt_norm_t *X, int len, float theta, int dir, int stride, int iter)
47 c = Q15ONE*cos(theta);
48 s = dir*Q15ONE*sin(theta);
51 /* We could use MULT16_16_P15 instead of MULT16_16_Q15 for more accuracy,
52 but at this point, I really don't think it's necessary */
53 for (i=0;i<len-stride;i++)
58 X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
59 X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
61 for (i=len-2*stride-1;i>=0;i--)
66 X[i] = MULT16_16_Q15(c,x1) - MULT16_16_Q15(s,x2);
67 X[i+stride] = MULT16_16_Q15(c,x2) + MULT16_16_Q15(s,x1);
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)
76 const int *eBands = m->eBands;
81 for (i=0;i<m->nbEBands;i++)
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]);*/
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)
98 const int *eBands = m->eBands;
103 for (i=0;i<m->nbEBands;i++)
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;
111 for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
116 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
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);
129 void renormalise_bands(const CELTMode *m, celt_norm_t *X)
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);
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)
142 const int *eBands = m->eBands;
147 for (i=0;i<m->nbEBands;i++)
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;
155 for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)
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)
164 const int *pBands = m->pBands;
165 B = m->nbMdctBlocks*m->nbChannels;
167 for (i=0;i<m->nbPBands;i++)
169 celt_word32_t Sxy=0, Sxx=0;
171 /* We know we're not going to overflow because Sxx can't be more than 1 (Q28) */
172 for (j=B*pBands[i];j<B*pBands[i+1];j++)
174 Sxy = MAC16_16(Sxy, X[j], P[j]);
175 Sxx = MAC16_16(Sxx, X[j], X[j]);
177 /* No negative gain allowed */
180 /* Not sure how that would happen, just making sure */
183 /* We need to be a bit conservative (multiply gain by 0.9), otherwise the
184 residual doesn't quantise well */
185 Sxy = MULT16_32_Q15(QCONST16(.9f, 15), Sxy);
187 gains[i] = DIV32_16(Sxy,ADD32(SHR32(Sxx, PGAIN_SHIFT),EPSILON));
188 /*printf ("%f ", 1-sqrt(1-gain*gain));*/
192 for (i=0;i<m->nbPBands;i++)
193 printf ("%f ", 1-sqrt(1-gains[i]*gains[i]));
196 for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
200 /* Apply the (quantised) gain to each "pitch band" */
201 void pitch_quant_bands(const CELTMode *m, celt_norm_t *P, celt_pgain_t *gains)
204 const int *pBands = m->pBands;
205 B = m->nbMdctBlocks*m->nbChannels;
206 for (i=0;i<m->nbPBands;i++)
209 for (j=B*pBands[i];j<B*pBands[i+1];j++)
210 P[j] = MULT16_16_Q15(gains[i], P[j]);
211 /*printf ("%f ", gain);*/
213 for (i=B*pBands[m->nbPBands];i<B*pBands[m->nbPBands+1];i++)
218 /* Quantisation of the residual */
219 void quant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, celt_mask_t *W, int total_bits, ec_enc *enc)
222 const int *eBands = m->eBands;
224 VARDECL(celt_norm_t *norm);
225 VARDECL(int *pulses);
226 VARDECL(int *offsets);
228 B = m->nbMdctBlocks*m->nbChannels;
230 ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
231 ALLOC(pulses, m->nbEBands, int);
232 ALLOC(offsets, m->nbEBands, int);
234 for (i=0;i<m->nbEBands;i++)
236 /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
237 bits = total_bits - ec_enc_tell(enc, 0) - 1;
238 compute_allocation(m, offsets, bits, pulses);
240 /*printf("bits left: %d\n", bits);
241 for (i=0;i<m->nbEBands;i++)
242 printf ("%d ", pulses[i]);
244 /*printf ("%d %d\n", ec_enc_tell(enc, 0), compute_allocation(m, m->nbPulses));*/
245 for (i=0;i<m->nbEBands;i++)
250 /*Scale factor of .0625f is just there to prevent overflows in fixed-point
251 (has no effect on float)*/
252 n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
253 theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
255 /* If pitch isn't available, use intra-frame prediction */
256 if (eBands[i] >= m->pitchEnd || q<=0)
261 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 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 alpha = QCONST16(.7f,15);
270 exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
271 exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
272 alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, enc);
273 exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
275 for (j=B*eBands[i];j<B*eBands[i+1];j++)
278 for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
282 /* Decoding of the residual */
283 void unquant_bands(const CELTMode *m, celt_norm_t *X, celt_norm_t *P, int total_bits, ec_dec *dec)
286 const int *eBands = m->eBands;
288 VARDECL(celt_norm_t *norm);
289 VARDECL(int *pulses);
290 VARDECL(int *offsets);
292 B = m->nbMdctBlocks*m->nbChannels;
294 ALLOC(norm, B*eBands[m->nbEBands+1], celt_norm_t);
295 ALLOC(pulses, m->nbEBands, int);
296 ALLOC(offsets, m->nbEBands, int);
298 for (i=0;i<m->nbEBands;i++)
300 /* Use a single-bit margin to guard against overrunning (make sure it's enough) */
301 bits = total_bits - ec_dec_tell(dec, 0) - 1;
302 compute_allocation(m, offsets, bits, pulses);
304 for (i=0;i<m->nbEBands;i++)
309 /*Scale factor of .0625f is just there to prevent overflows in fixed-point
310 (has no effect on float)*/
311 n = .0625f*sqrt(B*(eBands[i+1]-eBands[i]));
312 theta = .007*(B*(eBands[i+1]-eBands[i]))/(.1f+q);
314 /* If pitch isn't available, use intra-frame prediction */
315 if (eBands[i] >= m->pitchEnd || q<=0)
320 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 intra_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, norm, P+B*eBands[i], B, eBands[i], dec);
324 alpha = QCONST16(.7f,15);
329 exp_rotation(P+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, -1, B, 8);
330 alg_unquant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], alpha, dec);
331 exp_rotation(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), theta, 1, B, 8);
333 for (j=B*eBands[i];j<B*eBands[i+1];j++)
336 for (i=B*eBands[m->nbEBands];i<B*eBands[m->nbEBands+1];i++)
340 void stereo_mix(const CELTMode *m, celt_norm_t *X, celt_ener_t *bank, int dir)
343 const int *eBands = m->eBands;
346 for (i=0;i<m->nbEBands;i++)
353 a1 = left/sqrt(.01+left*left+right*right);
354 a2 = dir*right/sqrt(.01+left*left+right*right);
355 for (j=B*eBands[i];j<B*eBands[i+1];j++)
360 X[j*C] = a1*l + a2*r;
361 X[j*C+1] = a1*r - a2*l;
364 for (i=B*C*eBands[m->nbEBands];i<B*C*eBands[m->nbEBands+1];i++)