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