FFmpeg
opus_pvq.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2007-2008 CSIRO
3  * Copyright (c) 2007-2009 Xiph.Org Foundation
4  * Copyright (c) 2008-2009 Gregory Maxwell
5  * Copyright (c) 2012 Andrew D'Addesio
6  * Copyright (c) 2013-2014 Mozilla Corporation
7  * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker@gmail.com>
8  *
9  * This file is part of FFmpeg.
10  *
11  * FFmpeg is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU Lesser General Public
13  * License as published by the Free Software Foundation; either
14  * version 2.1 of the License, or (at your option) any later version.
15  *
16  * FFmpeg is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  * Lesser General Public License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with FFmpeg; if not, write to the Free Software
23  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24  */
25 
26 #include "opustab.h"
27 #include "opus_pvq.h"
28 
29 #define CELT_PVQ_U(n, k) (ff_celt_pvq_u_row[FFMIN(n, k)][FFMAX(n, k)])
30 #define CELT_PVQ_V(n, k) (CELT_PVQ_U(n, k) + CELT_PVQ_U(n, (k) + 1))
31 
32 static inline int16_t celt_cos(int16_t x)
33 {
34  x = (MUL16(x, x) + 4096) >> 13;
35  x = (32767-x) + ROUND_MUL16(x, (-7651 + ROUND_MUL16(x, (8277 + ROUND_MUL16(-626, x)))));
36  return x + 1;
37 }
38 
39 static inline int celt_log2tan(int isin, int icos)
40 {
41  int lc, ls;
42  lc = opus_ilog(icos);
43  ls = opus_ilog(isin);
44  icos <<= 15 - lc;
45  isin <<= 15 - ls;
46  return (ls << 11) - (lc << 11) +
47  ROUND_MUL16(isin, ROUND_MUL16(isin, -2597) + 7932) -
48  ROUND_MUL16(icos, ROUND_MUL16(icos, -2597) + 7932);
49 }
50 
51 static inline int celt_bits2pulses(const uint8_t *cache, int bits)
52 {
53  // TODO: Find the size of cache and make it into an array in the parameters list
54  int i, low = 0, high;
55 
56  high = cache[0];
57  bits--;
58 
59  for (i = 0; i < 6; i++) {
60  int center = (low + high + 1) >> 1;
61  if (cache[center] >= bits)
62  high = center;
63  else
64  low = center;
65  }
66 
67  return (bits - (low == 0 ? -1 : cache[low]) <= cache[high] - bits) ? low : high;
68 }
69 
70 static inline int celt_pulses2bits(const uint8_t *cache, int pulses)
71 {
72  // TODO: Find the size of cache and make it into an array in the parameters list
73  return (pulses == 0) ? 0 : cache[pulses] + 1;
74 }
75 
76 static inline void celt_normalize_residual(const int * av_restrict iy, float * av_restrict X,
77  int N, float g)
78 {
79  int i;
80  for (i = 0; i < N; i++)
81  X[i] = g * iy[i];
82 }
83 
84 static void celt_exp_rotation_impl(float *X, uint32_t len, uint32_t stride,
85  float c, float s)
86 {
87  float *Xptr;
88  int i;
89 
90  Xptr = X;
91  for (i = 0; i < len - stride; i++) {
92  float x1 = Xptr[0];
93  float x2 = Xptr[stride];
94  Xptr[stride] = c * x2 + s * x1;
95  *Xptr++ = c * x1 - s * x2;
96  }
97 
98  Xptr = &X[len - 2 * stride - 1];
99  for (i = len - 2 * stride - 1; i >= 0; i--) {
100  float x1 = Xptr[0];
101  float x2 = Xptr[stride];
102  Xptr[stride] = c * x2 + s * x1;
103  *Xptr-- = c * x1 - s * x2;
104  }
105 }
106 
107 static inline void celt_exp_rotation(float *X, uint32_t len,
108  uint32_t stride, uint32_t K,
109  enum CeltSpread spread, const int encode)
110 {
111  uint32_t stride2 = 0;
112  float c, s;
113  float gain, theta;
114  int i;
115 
116  if (2*K >= len || spread == CELT_SPREAD_NONE)
117  return;
118 
119  gain = (float)len / (len + (20 - 5*spread) * K);
120  theta = M_PI * gain * gain / 4;
121 
122  c = cosf(theta);
123  s = sinf(theta);
124 
125  if (len >= stride << 3) {
126  stride2 = 1;
127  /* This is just a simple (equivalent) way of computing sqrt(len/stride) with rounding.
128  It's basically incrementing long as (stride2+0.5)^2 < len/stride. */
129  while ((stride2 * stride2 + stride2) * stride + (stride >> 2) < len)
130  stride2++;
131  }
132 
133  len /= stride;
134  for (i = 0; i < stride; i++) {
135  if (encode) {
136  celt_exp_rotation_impl(X + i * len, len, 1, c, -s);
137  if (stride2)
138  celt_exp_rotation_impl(X + i * len, len, stride2, s, -c);
139  } else {
140  if (stride2)
141  celt_exp_rotation_impl(X + i * len, len, stride2, s, c);
142  celt_exp_rotation_impl(X + i * len, len, 1, c, s);
143  }
144  }
145 }
146 
147 static inline uint32_t celt_extract_collapse_mask(const int *iy, uint32_t N, uint32_t B)
148 {
149  int i, j, N0 = N / B;
150  uint32_t collapse_mask = 0;
151 
152  if (B <= 1)
153  return 1;
154 
155  for (i = 0; i < B; i++)
156  for (j = 0; j < N0; j++)
157  collapse_mask |= (!!iy[i*N0+j]) << i;
158  return collapse_mask;
159 }
160 
161 static inline void celt_stereo_merge(float *X, float *Y, float mid, int N)
162 {
163  int i;
164  float xp = 0, side = 0;
165  float E[2];
166  float mid2;
167  float gain[2];
168 
169  /* Compute the norm of X+Y and X-Y as |X|^2 + |Y|^2 +/- sum(xy) */
170  for (i = 0; i < N; i++) {
171  xp += X[i] * Y[i];
172  side += Y[i] * Y[i];
173  }
174 
175  /* Compensating for the mid normalization */
176  xp *= mid;
177  mid2 = mid;
178  E[0] = mid2 * mid2 + side - 2 * xp;
179  E[1] = mid2 * mid2 + side + 2 * xp;
180  if (E[0] < 6e-4f || E[1] < 6e-4f) {
181  for (i = 0; i < N; i++)
182  Y[i] = X[i];
183  return;
184  }
185 
186  gain[0] = 1.0f / sqrtf(E[0]);
187  gain[1] = 1.0f / sqrtf(E[1]);
188 
189  for (i = 0; i < N; i++) {
190  float value[2];
191  /* Apply mid scaling (side is already scaled) */
192  value[0] = mid * X[i];
193  value[1] = Y[i];
194  X[i] = gain[0] * (value[0] - value[1]);
195  Y[i] = gain[1] * (value[0] + value[1]);
196  }
197 }
198 
199 static void celt_interleave_hadamard(float *tmp, float *X, int N0,
200  int stride, int hadamard)
201 {
202  int i, j, N = N0*stride;
203  const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
204 
205  for (i = 0; i < stride; i++)
206  for (j = 0; j < N0; j++)
207  tmp[j*stride+i] = X[order[i]*N0+j];
208 
209  memcpy(X, tmp, N*sizeof(float));
210 }
211 
212 static void celt_deinterleave_hadamard(float *tmp, float *X, int N0,
213  int stride, int hadamard)
214 {
215  int i, j, N = N0*stride;
216  const uint8_t *order = &ff_celt_hadamard_order[hadamard ? stride - 2 : 30];
217 
218  for (i = 0; i < stride; i++)
219  for (j = 0; j < N0; j++)
220  tmp[order[i]*N0+j] = X[j*stride+i];
221 
222  memcpy(X, tmp, N*sizeof(float));
223 }
224 
225 static void celt_haar1(float *X, int N0, int stride)
226 {
227  int i, j;
228  N0 >>= 1;
229  for (i = 0; i < stride; i++) {
230  for (j = 0; j < N0; j++) {
231  float x0 = X[stride * (2 * j + 0) + i];
232  float x1 = X[stride * (2 * j + 1) + i];
233  X[stride * (2 * j + 0) + i] = (x0 + x1) * M_SQRT1_2;
234  X[stride * (2 * j + 1) + i] = (x0 - x1) * M_SQRT1_2;
235  }
236  }
237 }
238 
239 static inline int celt_compute_qn(int N, int b, int offset, int pulse_cap,
240  int stereo)
241 {
242  int qn, qb;
243  int N2 = 2 * N - 1;
244  if (stereo && N == 2)
245  N2--;
246 
247  /* The upper limit ensures that in a stereo split with itheta==16384, we'll
248  * always have enough bits left over to code at least one pulse in the
249  * side; otherwise it would collapse, since it doesn't get folded. */
250  qb = FFMIN3(b - pulse_cap - (4 << 3), (b + N2 * offset) / N2, 8 << 3);
251  qn = (qb < (1 << 3 >> 1)) ? 1 : ((ff_celt_qn_exp2[qb & 0x7] >> (14 - (qb >> 3))) + 1) >> 1 << 1;
252  return qn;
253 }
254 
255 /* Convert the quantized vector to an index */
256 static inline uint32_t celt_icwrsi(uint32_t N, uint32_t K, const int *y)
257 {
258  int i, idx = 0, sum = 0;
259  for (i = N - 1; i >= 0; i--) {
260  const uint32_t i_s = CELT_PVQ_U(N - i, sum + FFABS(y[i]) + 1);
261  idx += CELT_PVQ_U(N - i, sum) + (y[i] < 0)*i_s;
262  sum += FFABS(y[i]);
263  }
264  return idx;
265 }
266 
267 // this code was adapted from libopus
268 static inline uint64_t celt_cwrsi(uint32_t N, uint32_t K, uint32_t i, int *y)
269 {
270  uint64_t norm = 0;
271  uint32_t q, p;
272  int s, val;
273  int k0;
274 
275  while (N > 2) {
276  /*Lots of pulses case:*/
277  if (K >= N) {
278  const uint32_t *row = ff_celt_pvq_u_row[N];
279 
280  /* Are the pulses in this dimension negative? */
281  p = row[K + 1];
282  s = -(i >= p);
283  i -= p & s;
284 
285  /*Count how many pulses were placed in this dimension.*/
286  k0 = K;
287  q = row[N];
288  if (q > i) {
289  K = N;
290  do {
291  p = ff_celt_pvq_u_row[--K][N];
292  } while (p > i);
293  } else
294  for (p = row[K]; p > i; p = row[K])
295  K--;
296 
297  i -= p;
298  val = (k0 - K + s) ^ s;
299  norm += val * val;
300  *y++ = val;
301  } else { /*Lots of dimensions case:*/
302  /*Are there any pulses in this dimension at all?*/
303  p = ff_celt_pvq_u_row[K ][N];
304  q = ff_celt_pvq_u_row[K + 1][N];
305 
306  if (p <= i && i < q) {
307  i -= p;
308  *y++ = 0;
309  } else {
310  /*Are the pulses in this dimension negative?*/
311  s = -(i >= q);
312  i -= q & s;
313 
314  /*Count how many pulses were placed in this dimension.*/
315  k0 = K;
316  do p = ff_celt_pvq_u_row[--K][N];
317  while (p > i);
318 
319  i -= p;
320  val = (k0 - K + s) ^ s;
321  norm += val * val;
322  *y++ = val;
323  }
324  }
325  N--;
326  }
327 
328  /* N == 2 */
329  p = 2 * K + 1;
330  s = -(i >= p);
331  i -= p & s;
332  k0 = K;
333  K = (i + 1) / 2;
334 
335  if (K)
336  i -= 2 * K - 1;
337 
338  val = (k0 - K + s) ^ s;
339  norm += val * val;
340  *y++ = val;
341 
342  /* N==1 */
343  s = -i;
344  val = (K + s) ^ s;
345  norm += val * val;
346  *y = val;
347 
348  return norm;
349 }
350 
351 static inline void celt_encode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
352 {
353  ff_opus_rc_enc_uint(rc, celt_icwrsi(N, K, y), CELT_PVQ_V(N, K));
354 }
355 
356 static inline float celt_decode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
357 {
358  const uint32_t idx = ff_opus_rc_dec_uint(rc, CELT_PVQ_V(N, K));
359  return celt_cwrsi(N, K, idx, y);
360 }
361 
362 /*
363  * Faster than libopus's search, operates entirely in the signed domain.
364  * Slightly worse/better depending on N, K and the input vector.
365  */
366 static float ppp_pvq_search_c(float *X, int *y, int K, int N)
367 {
368  int i, y_norm = 0;
369  float res = 0.0f, xy_norm = 0.0f;
370 
371  for (i = 0; i < N; i++)
372  res += FFABS(X[i]);
373 
374  res = K/(res + FLT_EPSILON);
375 
376  for (i = 0; i < N; i++) {
377  y[i] = lrintf(res*X[i]);
378  y_norm += y[i]*y[i];
379  xy_norm += y[i]*X[i];
380  K -= FFABS(y[i]);
381  }
382 
383  while (K) {
384  int max_idx = 0, phase = FFSIGN(K);
385  float max_num = 0.0f;
386  float max_den = 1.0f;
387  y_norm += 1.0f;
388 
389  for (i = 0; i < N; i++) {
390  /* If the sum has been overshot and the best place has 0 pulses allocated
391  * to it, attempting to decrease it further will actually increase the
392  * sum. Prevent this by disregarding any 0 positions when decrementing. */
393  const int ca = 1 ^ ((y[i] == 0) & (phase < 0));
394  const int y_new = y_norm + 2*phase*FFABS(y[i]);
395  float xy_new = xy_norm + 1*phase*FFABS(X[i]);
396  xy_new = xy_new * xy_new;
397  if (ca && (max_den*xy_new) > (y_new*max_num)) {
398  max_den = y_new;
399  max_num = xy_new;
400  max_idx = i;
401  }
402  }
403 
404  K -= phase;
405 
406  phase *= FFSIGN(X[max_idx]);
407  xy_norm += 1*phase*X[max_idx];
408  y_norm += 2*phase*y[max_idx];
409  y[max_idx] += phase;
410  }
411 
412  return (float)y_norm;
413 }
414 
415 static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
416  enum CeltSpread spread, uint32_t blocks, float gain,
417  CeltPVQ *pvq)
418 {
419  int *y = pvq->qcoeff;
420 
421  celt_exp_rotation(X, N, blocks, K, spread, 1);
422  gain /= sqrtf(pvq->pvq_search(X, y, K, N));
423  celt_encode_pulses(rc, y, N, K);
424  celt_normalize_residual(y, X, N, gain);
425  celt_exp_rotation(X, N, blocks, K, spread, 0);
426  return celt_extract_collapse_mask(y, N, blocks);
427 }
428 
429 /** Decode pulse vector and combine the result with the pitch vector to produce
430  the final normalised signal in the current band. */
431 static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
432  enum CeltSpread spread, uint32_t blocks, float gain,
433  CeltPVQ *pvq)
434 {
435  int *y = pvq->qcoeff;
436 
437  gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
438  celt_normalize_residual(y, X, N, gain);
439  celt_exp_rotation(X, N, blocks, K, spread, 0);
440  return celt_extract_collapse_mask(y, N, blocks);
441 }
442 
443 static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
444 {
445  int i;
446  float e[2] = { 0.0f, 0.0f };
447  if (coupling) { /* Coupling case */
448  for (i = 0; i < N; i++) {
449  e[0] += (X[i] + Y[i])*(X[i] + Y[i]);
450  e[1] += (X[i] - Y[i])*(X[i] - Y[i]);
451  }
452  } else {
453  for (i = 0; i < N; i++) {
454  e[0] += X[i]*X[i];
455  e[1] += Y[i]*Y[i];
456  }
457  }
458  return lrintf(32768.0f*atan2f(sqrtf(e[1]), sqrtf(e[0]))/M_PI);
459 }
460 
461 static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
462 {
463  int i;
464  const float energy_n = 1.0f/(sqrtf(e_l*e_l + e_r*e_r) + FLT_EPSILON);
465  e_l *= energy_n;
466  e_r *= energy_n;
467  for (i = 0; i < N; i++)
468  X[i] = e_l*X[i] + e_r*Y[i];
469 }
470 
471 static void celt_stereo_ms_decouple(float *X, float *Y, int N)
472 {
473  int i;
474  for (i = 0; i < N; i++) {
475  const float Xret = X[i];
476  X[i] = (X[i] + Y[i])*M_SQRT1_2;
477  Y[i] = (Y[i] - Xret)*M_SQRT1_2;
478  }
479 }
480 
482  OpusRangeCoder *rc,
483  const int band, float *X,
484  float *Y, int N, int b,
485  uint32_t blocks, float *lowband,
486  int duration, float *lowband_out,
487  int level, float gain,
488  float *lowband_scratch,
489  int fill, int quant)
490 {
491  int i;
492  const uint8_t *cache;
493  int stereo = !!Y, split = stereo;
494  int imid = 0, iside = 0;
495  uint32_t N0 = N;
496  int N_B = N / blocks;
497  int N_B0 = N_B;
498  int B0 = blocks;
499  int time_divide = 0;
500  int recombine = 0;
501  int inv = 0;
502  float mid = 0, side = 0;
503  int longblocks = (B0 == 1);
504  uint32_t cm = 0;
505 
506  if (N == 1) {
507  float *x = X;
508  for (i = 0; i <= stereo; i++) {
509  int sign = 0;
510  if (f->remaining2 >= 1 << 3) {
511  if (quant) {
512  sign = x[0] < 0;
513  ff_opus_rc_put_raw(rc, sign, 1);
514  } else {
515  sign = ff_opus_rc_get_raw(rc, 1);
516  }
517  f->remaining2 -= 1 << 3;
518  }
519  x[0] = 1.0f - 2.0f*sign;
520  x = Y;
521  }
522  if (lowband_out)
523  lowband_out[0] = X[0];
524  return 1;
525  }
526 
527  if (!stereo && level == 0) {
528  int tf_change = f->tf_change[band];
529  int k;
530  if (tf_change > 0)
531  recombine = tf_change;
532  /* Band recombining to increase frequency resolution */
533 
534  if (lowband &&
535  (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
536  for (i = 0; i < N; i++)
537  lowband_scratch[i] = lowband[i];
538  lowband = lowband_scratch;
539  }
540 
541  for (k = 0; k < recombine; k++) {
542  if (quant || lowband)
543  celt_haar1(quant ? X : lowband, N >> k, 1 << k);
544  fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
545  }
546  blocks >>= recombine;
547  N_B <<= recombine;
548 
549  /* Increasing the time resolution */
550  while ((N_B & 1) == 0 && tf_change < 0) {
551  if (quant || lowband)
552  celt_haar1(quant ? X : lowband, N_B, blocks);
553  fill |= fill << blocks;
554  blocks <<= 1;
555  N_B >>= 1;
556  time_divide++;
557  tf_change++;
558  }
559  B0 = blocks;
560  N_B0 = N_B;
561 
562  /* Reorganize the samples in time order instead of frequency order */
563  if (B0 > 1 && (quant || lowband))
564  celt_deinterleave_hadamard(pvq->hadamard_tmp, quant ? X : lowband,
565  N_B >> recombine, B0 << recombine,
566  longblocks);
567  }
568 
569  /* If we need 1.5 more bit than we can produce, split the band in two. */
570  cache = ff_celt_cache_bits +
571  ff_celt_cache_index[(duration + 1) * CELT_MAX_BANDS + band];
572  if (!stereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
573  N >>= 1;
574  Y = X + N;
575  split = 1;
576  duration -= 1;
577  if (blocks == 1)
578  fill = (fill & 1) | (fill << 1);
579  blocks = (blocks + 1) >> 1;
580  }
581 
582  if (split) {
583  int qn;
584  int itheta = quant ? celt_calc_theta(X, Y, stereo, N) : 0;
585  int mbits, sbits, delta;
586  int qalloc;
587  int pulse_cap;
588  int offset;
589  int orig_fill;
590  int tell;
591 
592  /* Decide on the resolution to give to the split parameter theta */
593  pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
594  offset = (pulse_cap >> 1) - (stereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
596  qn = (stereo && band >= f->intensity_stereo) ? 1 :
597  celt_compute_qn(N, b, offset, pulse_cap, stereo);
598  tell = opus_rc_tell_frac(rc);
599  if (qn != 1) {
600  if (quant)
601  itheta = (itheta*qn + 8192) >> 14;
602  /* Entropy coding of the angle. We use a uniform pdf for the
603  * time split, a step for stereo, and a triangular one for the rest. */
604  if (quant) {
605  if (stereo && N > 2)
606  ff_opus_rc_enc_uint_step(rc, itheta, qn / 2);
607  else if (stereo || B0 > 1)
608  ff_opus_rc_enc_uint(rc, itheta, qn + 1);
609  else
610  ff_opus_rc_enc_uint_tri(rc, itheta, qn);
611  itheta = itheta * 16384 / qn;
612  if (stereo) {
613  if (itheta == 0)
614  celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
615  f->block[1].lin_energy[band], N);
616  else
617  celt_stereo_ms_decouple(X, Y, N);
618  }
619  } else {
620  if (stereo && N > 2)
621  itheta = ff_opus_rc_dec_uint_step(rc, qn / 2);
622  else if (stereo || B0 > 1)
623  itheta = ff_opus_rc_dec_uint(rc, qn+1);
624  else
625  itheta = ff_opus_rc_dec_uint_tri(rc, qn);
626  itheta = itheta * 16384 / qn;
627  }
628  } else if (stereo) {
629  if (quant) {
630  inv = itheta > 8192;
631  if (inv) {
632  for (i = 0; i < N; i++)
633  Y[i] *= -1;
634  }
635  celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
636  f->block[1].lin_energy[band], N);
637 
638  if (b > 2 << 3 && f->remaining2 > 2 << 3) {
639  ff_opus_rc_enc_log(rc, inv, 2);
640  } else {
641  inv = 0;
642  }
643  } else {
644  inv = (b > 2 << 3 && f->remaining2 > 2 << 3) ? ff_opus_rc_dec_log(rc, 2) : 0;
645  inv = f->apply_phase_inv ? inv : 0;
646  }
647  itheta = 0;
648  }
649  qalloc = opus_rc_tell_frac(rc) - tell;
650  b -= qalloc;
651 
652  orig_fill = fill;
653  if (itheta == 0) {
654  imid = 32767;
655  iside = 0;
656  fill = av_mod_uintp2(fill, blocks);
657  delta = -16384;
658  } else if (itheta == 16384) {
659  imid = 0;
660  iside = 32767;
661  fill &= ((1 << blocks) - 1) << blocks;
662  delta = 16384;
663  } else {
664  imid = celt_cos(itheta);
665  iside = celt_cos(16384-itheta);
666  /* This is the mid vs side allocation that minimizes squared error
667  in that band. */
668  delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
669  }
670 
671  mid = imid / 32768.0f;
672  side = iside / 32768.0f;
673 
674  /* This is a special case for N=2 that only works for stereo and takes
675  advantage of the fact that mid and side are orthogonal to encode
676  the side with just one bit. */
677  if (N == 2 && stereo) {
678  int c;
679  int sign = 0;
680  float tmp;
681  float *x2, *y2;
682  mbits = b;
683  /* Only need one bit for the side */
684  sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
685  mbits -= sbits;
686  c = (itheta > 8192);
687  f->remaining2 -= qalloc+sbits;
688 
689  x2 = c ? Y : X;
690  y2 = c ? X : Y;
691  if (sbits) {
692  if (quant) {
693  sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
694  ff_opus_rc_put_raw(rc, sign, 1);
695  } else {
696  sign = ff_opus_rc_get_raw(rc, 1);
697  }
698  }
699  sign = 1 - 2 * sign;
700  /* We use orig_fill here because we want to fold the side, but if
701  itheta==16384, we'll have cleared the low bits of fill. */
702  cm = pvq->quant_band(pvq, f, rc, band, x2, NULL, N, mbits, blocks, lowband, duration,
703  lowband_out, level, gain, lowband_scratch, orig_fill);
704  /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
705  and there's no need to worry about mixing with the other channel. */
706  y2[0] = -sign * x2[1];
707  y2[1] = sign * x2[0];
708  X[0] *= mid;
709  X[1] *= mid;
710  Y[0] *= side;
711  Y[1] *= side;
712  tmp = X[0];
713  X[0] = tmp - Y[0];
714  Y[0] = tmp + Y[0];
715  tmp = X[1];
716  X[1] = tmp - Y[1];
717  Y[1] = tmp + Y[1];
718  } else {
719  /* "Normal" split code */
720  float *next_lowband2 = NULL;
721  float *next_lowband_out1 = NULL;
722  int next_level = 0;
723  int rebalance;
724  uint32_t cmt;
725 
726  /* Give more bits to low-energy MDCTs than they would
727  * otherwise deserve */
728  if (B0 > 1 && !stereo && (itheta & 0x3fff)) {
729  if (itheta > 8192)
730  /* Rough approximation for pre-echo masking */
731  delta -= delta >> (4 - duration);
732  else
733  /* Corresponds to a forward-masking slope of
734  * 1.5 dB per 10 ms */
735  delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
736  }
737  mbits = av_clip((b - delta) / 2, 0, b);
738  sbits = b - mbits;
739  f->remaining2 -= qalloc;
740 
741  if (lowband && !stereo)
742  next_lowband2 = lowband + N; /* >32-bit split case */
743 
744  /* Only stereo needs to pass on lowband_out.
745  * Otherwise, it's handled at the end */
746  if (stereo)
747  next_lowband_out1 = lowband_out;
748  else
749  next_level = level + 1;
750 
751  rebalance = f->remaining2;
752  if (mbits >= sbits) {
753  /* In stereo mode, we do not apply a scaling to the mid
754  * because we need the normalized mid for folding later */
755  cm = pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
756  lowband, duration, next_lowband_out1, next_level,
757  stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
758  rebalance = mbits - (rebalance - f->remaining2);
759  if (rebalance > 3 << 3 && itheta != 0)
760  sbits += rebalance - (3 << 3);
761 
762  /* For a stereo split, the high bits of fill are always zero,
763  * so no folding will be done to the side. */
764  cmt = pvq->quant_band(pvq, f, rc, band, Y, NULL, N, sbits, blocks,
765  next_lowband2, duration, NULL, next_level,
766  gain * side, NULL, fill >> blocks);
767  cm |= cmt << ((B0 >> 1) & (stereo - 1));
768  } else {
769  /* For a stereo split, the high bits of fill are always zero,
770  * so no folding will be done to the side. */
771  cm = pvq->quant_band(pvq, f, rc, band, Y, NULL, N, sbits, blocks,
772  next_lowband2, duration, NULL, next_level,
773  gain * side, NULL, fill >> blocks);
774  cm <<= ((B0 >> 1) & (stereo - 1));
775  rebalance = sbits - (rebalance - f->remaining2);
776  if (rebalance > 3 << 3 && itheta != 16384)
777  mbits += rebalance - (3 << 3);
778 
779  /* In stereo mode, we do not apply a scaling to the mid because
780  * we need the normalized mid for folding later */
781  cm |= pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
782  lowband, duration, next_lowband_out1, next_level,
783  stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
784  }
785  }
786  } else {
787  /* This is the basic no-split case */
788  uint32_t q = celt_bits2pulses(cache, b);
789  uint32_t curr_bits = celt_pulses2bits(cache, q);
790  f->remaining2 -= curr_bits;
791 
792  /* Ensures we can never bust the budget */
793  while (f->remaining2 < 0 && q > 0) {
794  f->remaining2 += curr_bits;
795  curr_bits = celt_pulses2bits(cache, --q);
796  f->remaining2 -= curr_bits;
797  }
798 
799  if (q != 0) {
800  /* Finally do the actual (de)quantization */
801  if (quant) {
802  cm = celt_alg_quant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
803  f->spread, blocks, gain, pvq);
804  } else {
805  cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
806  f->spread, blocks, gain, pvq);
807  }
808  } else {
809  /* If there's no pulse, fill the band anyway */
810  uint32_t cm_mask = (1 << blocks) - 1;
811  fill &= cm_mask;
812  if (fill) {
813  if (!lowband) {
814  /* Noise */
815  for (i = 0; i < N; i++)
816  X[i] = (((int32_t)celt_rng(f)) >> 20);
817  cm = cm_mask;
818  } else {
819  /* Folded spectrum */
820  for (i = 0; i < N; i++) {
821  /* About 48 dB below the "normal" folding level */
822  X[i] = lowband[i] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
823  }
824  cm = fill;
825  }
826  celt_renormalize_vector(X, N, gain);
827  } else {
828  memset(X, 0, N*sizeof(float));
829  }
830  }
831  }
832 
833  /* This code is used by the decoder and by the resynthesis-enabled encoder */
834  if (stereo) {
835  if (N > 2)
836  celt_stereo_merge(X, Y, mid, N);
837  if (inv) {
838  for (i = 0; i < N; i++)
839  Y[i] *= -1;
840  }
841  } else if (level == 0) {
842  int k;
843 
844  /* Undo the sample reorganization going from time order to frequency order */
845  if (B0 > 1)
846  celt_interleave_hadamard(pvq->hadamard_tmp, X, N_B >> recombine,
847  B0 << recombine, longblocks);
848 
849  /* Undo time-freq changes that we did earlier */
850  N_B = N_B0;
851  blocks = B0;
852  for (k = 0; k < time_divide; k++) {
853  blocks >>= 1;
854  N_B <<= 1;
855  cm |= cm >> blocks;
856  celt_haar1(X, N_B, blocks);
857  }
858 
859  for (k = 0; k < recombine; k++) {
861  celt_haar1(X, N0>>k, 1<<k);
862  }
863  blocks <<= recombine;
864 
865  /* Scale output for later folding */
866  if (lowband_out) {
867  float n = sqrtf(N0);
868  for (i = 0; i < N0; i++)
869  lowband_out[i] = n * X[i];
870  }
871  cm = av_mod_uintp2(cm, blocks);
872  }
873 
874  return cm;
875 }
876 
877 static QUANT_FN(pvq_decode_band)
878 {
879 #if CONFIG_OPUS_DECODER
880  return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
881  lowband_out, level, gain, lowband_scratch, fill, 0);
882 #else
883  return 0;
884 #endif
885 }
886 
887 static QUANT_FN(pvq_encode_band)
888 {
889 #if CONFIG_OPUS_ENCODER
890  return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
891  lowband_out, level, gain, lowband_scratch, fill, 1);
892 #else
893  return 0;
894 #endif
895 }
896 
898 {
899  CeltPVQ *s = av_malloc(sizeof(CeltPVQ));
900  if (!s)
901  return AVERROR(ENOMEM);
902 
904  s->quant_band = encode ? pvq_encode_band : pvq_decode_band;
905 
906  if (CONFIG_OPUS_ENCODER && ARCH_X86)
908 
909  *pvq = s;
910 
911  return 0;
912 }
913 
915 {
916  av_freep(pvq);
917 }
static void celt_stereo_merge(float *X, float *Y, float mid, int N)
Definition: opus_pvq.c:161
#define NULL
Definition: coverity.c:32
const char const char void * val
Definition: avisynth_c.h:863
const uint8_t ff_celt_cache_bits[392]
Definition: opustab.c:885
const uint32_t *const ff_celt_pvq_u_row[15]
Definition: opustab.c:1152
static int celt_compute_qn(int N, int b, int offset, int pulse_cap, int stereo)
Definition: opus_pvq.c:239
const uint8_t ff_celt_log_freq_range[]
Definition: opustab.c:771
int remaining2
Definition: opus_celt.h:134
void ff_celt_pvq_init_x86(struct CeltPVQ *s)
Definition: celt_pvq_init.c:31
const char * g
Definition: vf_curves.c:115
static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
Definition: opus_pvq.c:443
#define M_SQRT1_2
Definition: mathematics.h:58
static float ppp_pvq_search_c(float *X, int *y, int K, int N)
Definition: opus_pvq.c:366
static const int8_t pulses[4]
Number of non-zero pulses in the MP-MLQ excitation.
Definition: g723_1.h:723
Definition: vf_addroi.c:26
void ff_opus_rc_enc_log(OpusRangeCoder *rc, int val, uint32_t bits)
Definition: opus_rc.c:131
uint32_t ff_opus_rc_dec_log(OpusRangeCoder *rc, uint32_t bits)
Definition: opus_rc.c:114
int av_cold ff_celt_pvq_init(CeltPVQ **pvq, int encode)
Definition: opus_pvq.c:897
const uint8_t ff_celt_bit_deinterleave[]
Definition: opustab.c:933
static uint32_t celt_icwrsi(uint32_t N, uint32_t K, const int *y)
Definition: opus_pvq.c:256
#define N
Definition: af_mcompand.c:54
const uint16_t ff_celt_qn_exp2[]
Definition: opustab.c:946
void ff_opus_rc_enc_uint(OpusRangeCoder *rc, uint32_t val, uint32_t size)
CELT: write a uniformly distributed integer.
Definition: opus_rc.c:204
static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K, enum CeltSpread spread, uint32_t blocks, float gain, CeltPVQ *pvq)
Decode pulse vector and combine the result with the pitch vector to produce the final normalised sign...
Definition: opus_pvq.c:431
static void celt_interleave_hadamard(float *tmp, float *X, int N0, int stride, int hadamard)
Definition: opus_pvq.c:199
const uint8_t ff_celt_hadamard_order[]
Definition: opustab.c:938
uint32_t ff_opus_rc_dec_uint_tri(OpusRangeCoder *rc, int qn)
Definition: opus_rc.c:234
#define opus_ilog(i)
Definition: opus_rc.h:31
CeltBlock block[2]
Definition: opus_celt.h:97
uint8_t
static void celt_stereo_ms_decouple(float *X, float *Y, int N)
Definition: opus_pvq.c:471
#define av_cold
Definition: attributes.h:82
#define av_malloc(s)
float delta
it s the only field you need to keep assuming you have a context There is some magic you don t need to care about around this just let it vf offset
#define f(width, name)
Definition: cbs_vp9.c:255
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
GLsizei GLboolean const GLfloat * value
Definition: opengl_enc.c:108
int qcoeff[256]
Definition: opus_pvq.h:36
#define cosf(x)
Definition: libm.h:78
int64_t duration
Definition: movenc.c:63
static void celt_exp_rotation_impl(float *X, uint32_t len, uint32_t stride, float c, float s)
Definition: opus_pvq.c:84
#define FFMIN3(a, b, c)
Definition: common.h:97
#define atan2f(y, x)
Definition: libm.h:45
static void celt_deinterleave_hadamard(float *tmp, float *X, int N0, int stride, int hadamard)
Definition: opus_pvq.c:212
static float celt_decode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
Definition: opus_pvq.c:356
#define lrintf(x)
Definition: libm_mips.h:70
#define N2
Definition: vf_pp7.c:69
static void celt_haar1(float *X, int N0, int stride)
Definition: opus_pvq.c:225
float lin_energy[CELT_MAX_BANDS]
Definition: opus_celt.h:67
#define cm
Definition: dvbsubdec.c:37
void ff_opus_rc_enc_uint_tri(OpusRangeCoder *rc, uint32_t k, int qn)
Definition: opus_rc.c:258
const int16_t ff_celt_cache_index[105]
Definition: opustab.c:915
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:259
static int celt_log2tan(int isin, int icos)
Definition: opus_pvq.c:39
static void celt_encode_pulses(OpusRangeCoder *rc, int *y, uint32_t N, uint32_t K)
Definition: opus_pvq.c:351
#define B
Definition: huffyuvdsp.h:32
int tf_change[CELT_MAX_BANDS]
Definition: opus_celt.h:139
uint8_t bits
Definition: vp3data.h:202
int apply_phase_inv
Definition: opus_celt.h:102
#define MUL16(ra, rb)
Definition: mathops.h:88
static av_always_inline uint32_t quant_band_template(CeltPVQ *pvq, CeltFrame *f, OpusRangeCoder *rc, const int band, float *X, float *Y, int N, int b, uint32_t blocks, float *lowband, int duration, float *lowband_out, int level, float gain, float *lowband_scratch, int fill, int quant)
Definition: opus_pvq.c:481
float hadamard_tmp[256]
Definition: opus_pvq.h:37
static char * split(char *message, char delim)
Definition: af_channelmap.c:81
#define b
Definition: input.c:41
#define Y
Definition: boxblur.h:38
#define E
Definition: avdct.c:32
QUANT_FN * quant_band
Definition: opus_pvq.h:40
#define FFMIN(a, b)
Definition: common.h:96
#define FFSIGN(a)
Definition: common.h:73
static int16_t celt_cos(int16_t x)
Definition: opus_pvq.c:32
#define N0
Definition: vf_pp7.c:67
int32_t
static uint64_t celt_cwrsi(uint32_t N, uint32_t K, uint32_t i, int *y)
Definition: opus_pvq.c:268
#define FFABS(a)
Absolute value, Note, INT_MIN / INT64_MIN result in undefined behavior as they are not representable ...
Definition: common.h:72
#define s(width, name)
Definition: cbs_vp9.c:257
static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K, enum CeltSpread spread, uint32_t blocks, float gain, CeltPVQ *pvq)
Definition: opus_pvq.c:415
int n
Definition: avisynth_c.h:760
#define sinf(x)
Definition: libm.h:419
#define CELT_MAX_BANDS
Definition: opus.h:45
static void encode(AVCodecContext *ctx, AVFrame *frame, AVPacket *pkt, FILE *output)
Definition: encode_audio.c:95
#define CELT_PVQ_U(n, k)
Definition: opus_pvq.c:29
uint32_t ff_opus_rc_get_raw(OpusRangeCoder *rc, uint32_t count)
CELT: read 1-25 raw bits at the end of the frame, backwards byte-wise.
Definition: opus_rc.c:140
uint32_t ff_opus_rc_dec_uint_step(OpusRangeCoder *rc, int k0)
Definition: opus_rc.c:211
CeltSpread
Definition: opus_celt.h:49
static av_always_inline void celt_renormalize_vector(float *X, int N, float gain)
Definition: opus_celt.h:149
static uint32_t celt_extract_collapse_mask(const int *iy, uint32_t N, uint32_t B)
Definition: opus_pvq.c:147
static int celt_pulses2bits(const uint8_t *cache, int pulses)
Definition: opus_pvq.c:70
const uint8_t * quant
#define CELT_PVQ_V(n, k)
Definition: opus_pvq.c:30
uint8_t level
Definition: svq3.c:207
#define CELT_QTHETA_OFFSET
Definition: opus_celt.h:42
static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
Definition: opus_pvq.c:461
GLint GLenum GLboolean GLsizei stride
Definition: opengl_enc.c:104
#define ROUND_MUL16(a, b)
Definition: opus.h:51
void ff_opus_rc_put_raw(OpusRangeCoder *rc, uint32_t val, uint32_t count)
CELT: write 0 - 31 bits to the rawbits buffer.
Definition: opus_rc.c:161
static int celt_bits2pulses(const uint8_t *cache, int bits)
Definition: opus_pvq.c:51
static av_always_inline uint32_t celt_rng(CeltFrame *f)
Definition: opus_celt.h:143
enum CeltSpread spread
Definition: opus_celt.h:123
const uint8_t ff_celt_bit_interleave[]
Definition: opustab.c:929
int len
static void celt_exp_rotation(float *X, uint32_t len, uint32_t stride, uint32_t K, enum CeltSpread spread, const int encode)
Definition: opus_pvq.c:107
static void celt_normalize_residual(const int *av_restrict iy, float *av_restrict X, int N, float g)
Definition: opus_pvq.c:76
uint32_t ff_opus_rc_dec_uint(OpusRangeCoder *rc, uint32_t size)
CELT: read a uniform distribution.
Definition: opus_rc.c:182
#define av_freep(p)
#define av_always_inline
Definition: attributes.h:39
#define M_PI
Definition: mathematics.h:52
void ff_opus_rc_enc_uint_step(OpusRangeCoder *rc, uint32_t val, int k0)
Definition: opus_rc.c:226
#define stride
static QUANT_FN(pvq_decode_band)
Definition: opus_pvq.c:877
Filter the word “frame” indicates either a video frame or a group of audio as stored in an AVFrame structure Format for each input and each output the list of supported formats For video that means pixel format For audio that means channel sample they are references to shared objects When the negotiation mechanism computes the intersection of the formats supported at each end of a all references to both lists are replaced with a reference to the intersection And when a single format is eventually chosen for a link amongst the remaining all references to the list are updated That means that if a filter requires that its input and output have the same format amongst a supported all it has to do is use a reference to the same list of formats query_formats can leave some formats unset and return AVERROR(EAGAIN) to cause the negotiation mechanism toagain later.That can be used by filters with complex requirements to use the format negotiated on one link to set the formats supported on another.Frame references ownership and permissions
#define B0
Definition: faandct.c:40
void av_cold ff_celt_pvq_uninit(CeltPVQ **pvq)
Definition: opus_pvq.c:914
int intensity_stereo
Definition: opus_celt.h:119
#define CELT_QTHETA_OFFSET_TWOPHASE
Definition: opus_celt.h:43
float(* pvq_search)(float *X, int *y, int K, int N)
Definition: opus_pvq.h:39
static av_always_inline uint32_t opus_rc_tell_frac(const OpusRangeCoder *rc)
Definition: opus_rc.h:66
static uint8_t tmp[11]
Definition: aes_ctr.c:26