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 /*
370  * Faster than libopus's search, operates entirely in the signed domain.
371  * Slightly worse/better depending on N, K and the input vector.
372  */
373 static float ppp_pvq_search_c(float *X, int *y, int K, int N)
374 {
375  int i, y_norm = 0;
376  float res = 0.0f, xy_norm = 0.0f;
377 
378  for (i = 0; i < N; i++)
379  res += FFABS(X[i]);
380 
381  res = K/(res + FLT_EPSILON);
382 
383  for (i = 0; i < N; i++) {
384  y[i] = lrintf(res*X[i]);
385  y_norm += y[i]*y[i];
386  xy_norm += y[i]*X[i];
387  K -= FFABS(y[i]);
388  }
389 
390  while (K) {
391  int max_idx = 0, phase = FFSIGN(K);
392  float max_num = 0.0f;
393  float max_den = 1.0f;
394  y_norm += 1.0f;
395 
396  for (i = 0; i < N; i++) {
397  /* If the sum has been overshot and the best place has 0 pulses allocated
398  * to it, attempting to decrease it further will actually increase the
399  * sum. Prevent this by disregarding any 0 positions when decrementing. */
400  const int ca = 1 ^ ((y[i] == 0) & (phase < 0));
401  const int y_new = y_norm + 2*phase*FFABS(y[i]);
402  float xy_new = xy_norm + 1*phase*FFABS(X[i]);
403  xy_new = xy_new * xy_new;
404  if (ca && (max_den*xy_new) > (y_new*max_num)) {
405  max_den = y_new;
406  max_num = xy_new;
407  max_idx = i;
408  }
409  }
410 
411  K -= phase;
412 
413  phase *= FFSIGN(X[max_idx]);
414  xy_norm += 1*phase*X[max_idx];
415  y_norm += 2*phase*y[max_idx];
416  y[max_idx] += phase;
417  }
418 
419  return (float)y_norm;
420 }
421 
422 static uint32_t celt_alg_quant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
423  enum CeltSpread spread, uint32_t blocks, float gain,
424  CeltPVQ *pvq)
425 {
426  int *y = pvq->qcoeff;
427 
428  celt_exp_rotation(X, N, blocks, K, spread, 1);
429  gain /= sqrtf(pvq->pvq_search(X, y, K, N));
430  celt_encode_pulses(rc, y, N, K);
431  celt_normalize_residual(y, X, N, gain);
432  celt_exp_rotation(X, N, blocks, K, spread, 0);
433  return celt_extract_collapse_mask(y, N, blocks);
434 }
435 
436 /** Decode pulse vector and combine the result with the pitch vector to produce
437  the final normalised signal in the current band. */
438 static uint32_t celt_alg_unquant(OpusRangeCoder *rc, float *X, uint32_t N, uint32_t K,
439  enum CeltSpread spread, uint32_t blocks, float gain,
440  CeltPVQ *pvq)
441 {
442  int *y = pvq->qcoeff;
443 
444  gain /= sqrtf(celt_decode_pulses(rc, y, N, K));
445  celt_normalize_residual(y, X, N, gain);
446  celt_exp_rotation(X, N, blocks, K, spread, 0);
447  return celt_extract_collapse_mask(y, N, blocks);
448 }
449 
450 static int celt_calc_theta(const float *X, const float *Y, int coupling, int N)
451 {
452  int i;
453  float e[2] = { 0.0f, 0.0f };
454  if (coupling) { /* Coupling case */
455  for (i = 0; i < N; i++) {
456  e[0] += (X[i] + Y[i])*(X[i] + Y[i]);
457  e[1] += (X[i] - Y[i])*(X[i] - Y[i]);
458  }
459  } else {
460  for (i = 0; i < N; i++) {
461  e[0] += X[i]*X[i];
462  e[1] += Y[i]*Y[i];
463  }
464  }
465  return lrintf(32768.0f*atan2f(sqrtf(e[1]), sqrtf(e[0]))/M_PI);
466 }
467 
468 static void celt_stereo_is_decouple(float *X, float *Y, float e_l, float e_r, int N)
469 {
470  int i;
471  const float energy_n = 1.0f/(sqrtf(e_l*e_l + e_r*e_r) + FLT_EPSILON);
472  e_l *= energy_n;
473  e_r *= energy_n;
474  for (i = 0; i < N; i++)
475  X[i] = e_l*X[i] + e_r*Y[i];
476 }
477 
478 static void celt_stereo_ms_decouple(float *X, float *Y, int N)
479 {
480  int i;
481  for (i = 0; i < N; i++) {
482  const float Xret = X[i];
483  X[i] = (X[i] + Y[i])*M_SQRT1_2;
484  Y[i] = (Y[i] - Xret)*M_SQRT1_2;
485  }
486 }
487 
489  OpusRangeCoder *rc,
490  const int band, float *X,
491  float *Y, int N, int b,
492  uint32_t blocks, float *lowband,
493  int duration, float *lowband_out,
494  int level, float gain,
495  float *lowband_scratch,
496  int fill, int quant)
497 {
498  int i;
499  const uint8_t *cache;
500  int stereo = !!Y, split = stereo;
501  int imid = 0, iside = 0;
502  uint32_t N0 = N;
503  int N_B = N / blocks;
504  int N_B0 = N_B;
505  int B0 = blocks;
506  int time_divide = 0;
507  int recombine = 0;
508  int inv = 0;
509  float mid = 0, side = 0;
510  int longblocks = (B0 == 1);
511  uint32_t cm = 0;
512 
513  if (N == 1) {
514  float *x = X;
515  for (i = 0; i <= stereo; i++) {
516  int sign = 0;
517  if (f->remaining2 >= 1 << 3) {
518  if (quant) {
519  sign = x[0] < 0;
520  ff_opus_rc_put_raw(rc, sign, 1);
521  } else {
522  sign = ff_opus_rc_get_raw(rc, 1);
523  }
524  f->remaining2 -= 1 << 3;
525  }
526  x[0] = 1.0f - 2.0f*sign;
527  x = Y;
528  }
529  if (lowband_out)
530  lowband_out[0] = X[0];
531  return 1;
532  }
533 
534  if (!stereo && level == 0) {
535  int tf_change = f->tf_change[band];
536  int k;
537  if (tf_change > 0)
538  recombine = tf_change;
539  /* Band recombining to increase frequency resolution */
540 
541  if (lowband &&
542  (recombine || ((N_B & 1) == 0 && tf_change < 0) || B0 > 1)) {
543  for (i = 0; i < N; i++)
544  lowband_scratch[i] = lowband[i];
545  lowband = lowband_scratch;
546  }
547 
548  for (k = 0; k < recombine; k++) {
549  if (quant || lowband)
550  celt_haar1(quant ? X : lowband, N >> k, 1 << k);
551  fill = ff_celt_bit_interleave[fill & 0xF] | ff_celt_bit_interleave[fill >> 4] << 2;
552  }
553  blocks >>= recombine;
554  N_B <<= recombine;
555 
556  /* Increasing the time resolution */
557  while ((N_B & 1) == 0 && tf_change < 0) {
558  if (quant || lowband)
559  celt_haar1(quant ? X : lowband, N_B, blocks);
560  fill |= fill << blocks;
561  blocks <<= 1;
562  N_B >>= 1;
563  time_divide++;
564  tf_change++;
565  }
566  B0 = blocks;
567  N_B0 = N_B;
568 
569  /* Reorganize the samples in time order instead of frequency order */
570  if (B0 > 1 && (quant || lowband))
572  N_B >> recombine, B0 << recombine,
573  longblocks);
574  }
575 
576  /* If we need 1.5 more bit than we can produce, split the band in two. */
577  cache = ff_celt_cache_bits +
579  if (!stereo && duration >= 0 && b > cache[cache[0]] + 12 && N > 2) {
580  N >>= 1;
581  Y = X + N;
582  split = 1;
583  duration -= 1;
584  if (blocks == 1)
585  fill = (fill & 1) | (fill << 1);
586  blocks = (blocks + 1) >> 1;
587  }
588 
589  if (split) {
590  int qn;
591  int itheta = quant ? celt_calc_theta(X, Y, stereo, N) : 0;
592  int mbits, sbits, delta;
593  int qalloc;
594  int pulse_cap;
595  int offset;
596  int orig_fill;
597  int tell;
598 
599  /* Decide on the resolution to give to the split parameter theta */
600  pulse_cap = ff_celt_log_freq_range[band] + duration * 8;
601  offset = (pulse_cap >> 1) - (stereo && N == 2 ? CELT_QTHETA_OFFSET_TWOPHASE :
603  qn = (stereo && band >= f->intensity_stereo) ? 1 :
604  celt_compute_qn(N, b, offset, pulse_cap, stereo);
605  tell = opus_rc_tell_frac(rc);
606  if (qn != 1) {
607  if (quant)
608  itheta = (itheta*qn + 8192) >> 14;
609  /* Entropy coding of the angle. We use a uniform pdf for the
610  * time split, a step for stereo, and a triangular one for the rest. */
611  if (quant) {
612  if (stereo && N > 2)
613  ff_opus_rc_enc_uint_step(rc, itheta, qn / 2);
614  else if (stereo || B0 > 1)
615  ff_opus_rc_enc_uint(rc, itheta, qn + 1);
616  else
617  ff_opus_rc_enc_uint_tri(rc, itheta, qn);
618  itheta = itheta * 16384 / qn;
619  if (stereo) {
620  if (itheta == 0)
621  celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
622  f->block[1].lin_energy[band], N);
623  else
625  }
626  } else {
627  if (stereo && N > 2)
628  itheta = ff_opus_rc_dec_uint_step(rc, qn / 2);
629  else if (stereo || B0 > 1)
630  itheta = ff_opus_rc_dec_uint(rc, qn+1);
631  else
632  itheta = ff_opus_rc_dec_uint_tri(rc, qn);
633  itheta = itheta * 16384 / qn;
634  }
635  } else if (stereo) {
636  if (quant) {
637  inv = f->apply_phase_inv ? itheta > 8192 : 0;
638  if (inv) {
639  for (i = 0; i < N; i++)
640  Y[i] *= -1;
641  }
642  celt_stereo_is_decouple(X, Y, f->block[0].lin_energy[band],
643  f->block[1].lin_energy[band], N);
644 
645  if (b > 2 << 3 && f->remaining2 > 2 << 3) {
646  ff_opus_rc_enc_log(rc, inv, 2);
647  } else {
648  inv = 0;
649  }
650  } else {
651  inv = (b > 2 << 3 && f->remaining2 > 2 << 3) ? ff_opus_rc_dec_log(rc, 2) : 0;
652  inv = f->apply_phase_inv ? inv : 0;
653  }
654  itheta = 0;
655  }
656  qalloc = opus_rc_tell_frac(rc) - tell;
657  b -= qalloc;
658 
659  orig_fill = fill;
660  if (itheta == 0) {
661  imid = 32767;
662  iside = 0;
663  fill = av_mod_uintp2(fill, blocks);
664  delta = -16384;
665  } else if (itheta == 16384) {
666  imid = 0;
667  iside = 32767;
668  fill &= ((1 << blocks) - 1) << blocks;
669  delta = 16384;
670  } else {
671  imid = celt_cos(itheta);
672  iside = celt_cos(16384-itheta);
673  /* This is the mid vs side allocation that minimizes squared error
674  in that band. */
675  delta = ROUND_MUL16((N - 1) << 7, celt_log2tan(iside, imid));
676  }
677 
678  mid = imid / 32768.0f;
679  side = iside / 32768.0f;
680 
681  /* This is a special case for N=2 that only works for stereo and takes
682  advantage of the fact that mid and side are orthogonal to encode
683  the side with just one bit. */
684  if (N == 2 && stereo) {
685  int c;
686  int sign = 0;
687  float tmp;
688  float *x2, *y2;
689  mbits = b;
690  /* Only need one bit for the side */
691  sbits = (itheta != 0 && itheta != 16384) ? 1 << 3 : 0;
692  mbits -= sbits;
693  c = (itheta > 8192);
694  f->remaining2 -= qalloc+sbits;
695 
696  x2 = c ? Y : X;
697  y2 = c ? X : Y;
698  if (sbits) {
699  if (quant) {
700  sign = x2[0]*y2[1] - x2[1]*y2[0] < 0;
701  ff_opus_rc_put_raw(rc, sign, 1);
702  } else {
703  sign = ff_opus_rc_get_raw(rc, 1);
704  }
705  }
706  sign = 1 - 2 * sign;
707  /* We use orig_fill here because we want to fold the side, but if
708  itheta==16384, we'll have cleared the low bits of fill. */
709  cm = pvq->quant_band(pvq, f, rc, band, x2, NULL, N, mbits, blocks, lowband, duration,
710  lowband_out, level, gain, lowband_scratch, orig_fill);
711  /* We don't split N=2 bands, so cm is either 1 or 0 (for a fold-collapse),
712  and there's no need to worry about mixing with the other channel. */
713  y2[0] = -sign * x2[1];
714  y2[1] = sign * x2[0];
715  X[0] *= mid;
716  X[1] *= mid;
717  Y[0] *= side;
718  Y[1] *= side;
719  tmp = X[0];
720  X[0] = tmp - Y[0];
721  Y[0] = tmp + Y[0];
722  tmp = X[1];
723  X[1] = tmp - Y[1];
724  Y[1] = tmp + Y[1];
725  } else {
726  /* "Normal" split code */
727  float *next_lowband2 = NULL;
728  float *next_lowband_out1 = NULL;
729  int next_level = 0;
730  int rebalance;
731  uint32_t cmt;
732 
733  /* Give more bits to low-energy MDCTs than they would
734  * otherwise deserve */
735  if (B0 > 1 && !stereo && (itheta & 0x3fff)) {
736  if (itheta > 8192)
737  /* Rough approximation for pre-echo masking */
738  delta -= delta >> (4 - duration);
739  else
740  /* Corresponds to a forward-masking slope of
741  * 1.5 dB per 10 ms */
742  delta = FFMIN(0, delta + (N << 3 >> (5 - duration)));
743  }
744  mbits = av_clip((b - delta) / 2, 0, b);
745  sbits = b - mbits;
746  f->remaining2 -= qalloc;
747 
748  if (lowband && !stereo)
749  next_lowband2 = lowband + N; /* >32-bit split case */
750 
751  /* Only stereo needs to pass on lowband_out.
752  * Otherwise, it's handled at the end */
753  if (stereo)
754  next_lowband_out1 = lowband_out;
755  else
756  next_level = level + 1;
757 
758  rebalance = f->remaining2;
759  if (mbits >= sbits) {
760  /* In stereo mode, we do not apply a scaling to the mid
761  * because we need the normalized mid for folding later */
762  cm = pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
763  lowband, duration, next_lowband_out1, next_level,
764  stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
765  rebalance = mbits - (rebalance - f->remaining2);
766  if (rebalance > 3 << 3 && itheta != 0)
767  sbits += rebalance - (3 << 3);
768 
769  /* For a stereo split, the high bits of fill are always zero,
770  * so no folding will be done to the side. */
771  cmt = 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 |= cmt << ((B0 >> 1) & (stereo - 1));
775  } else {
776  /* For a stereo split, the high bits of fill are always zero,
777  * so no folding will be done to the side. */
778  cm = pvq->quant_band(pvq, f, rc, band, Y, NULL, N, sbits, blocks,
779  next_lowband2, duration, NULL, next_level,
780  gain * side, NULL, fill >> blocks);
781  cm <<= ((B0 >> 1) & (stereo - 1));
782  rebalance = sbits - (rebalance - f->remaining2);
783  if (rebalance > 3 << 3 && itheta != 16384)
784  mbits += rebalance - (3 << 3);
785 
786  /* In stereo mode, we do not apply a scaling to the mid because
787  * we need the normalized mid for folding later */
788  cm |= pvq->quant_band(pvq, f, rc, band, X, NULL, N, mbits, blocks,
789  lowband, duration, next_lowband_out1, next_level,
790  stereo ? 1.0f : (gain * mid), lowband_scratch, fill);
791  }
792  }
793  } else {
794  /* This is the basic no-split case */
795  uint32_t q = celt_bits2pulses(cache, b);
796  uint32_t curr_bits = celt_pulses2bits(cache, q);
797  f->remaining2 -= curr_bits;
798 
799  /* Ensures we can never bust the budget */
800  while (f->remaining2 < 0 && q > 0) {
801  f->remaining2 += curr_bits;
802  curr_bits = celt_pulses2bits(cache, --q);
803  f->remaining2 -= curr_bits;
804  }
805 
806  if (q != 0) {
807  /* Finally do the actual (de)quantization */
808  if (quant) {
809  cm = celt_alg_quant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
810  f->spread, blocks, gain, pvq);
811  } else {
812  cm = celt_alg_unquant(rc, X, N, (q < 8) ? q : (8 + (q & 7)) << ((q >> 3) - 1),
813  f->spread, blocks, gain, pvq);
814  }
815  } else {
816  /* If there's no pulse, fill the band anyway */
817  uint32_t cm_mask = (1 << blocks) - 1;
818  fill &= cm_mask;
819  if (fill) {
820  if (!lowband) {
821  /* Noise */
822  for (i = 0; i < N; i++)
823  X[i] = (((int32_t)celt_rng(f)) >> 20);
824  cm = cm_mask;
825  } else {
826  /* Folded spectrum */
827  for (i = 0; i < N; i++) {
828  /* About 48 dB below the "normal" folding level */
829  X[i] = lowband[i] + (((celt_rng(f)) & 0x8000) ? 1.0f / 256 : -1.0f / 256);
830  }
831  cm = fill;
832  }
833  celt_renormalize_vector(X, N, gain);
834  } else {
835  memset(X, 0, N*sizeof(float));
836  }
837  }
838  }
839 
840  /* This code is used by the decoder and by the resynthesis-enabled encoder */
841  if (stereo) {
842  if (N > 2)
843  celt_stereo_merge(X, Y, mid, N);
844  if (inv) {
845  for (i = 0; i < N; i++)
846  Y[i] *= -1;
847  }
848  } else if (level == 0) {
849  int k;
850 
851  /* Undo the sample reorganization going from time order to frequency order */
852  if (B0 > 1)
853  celt_interleave_hadamard(pvq->hadamard_tmp, X, N_B >> recombine,
854  B0 << recombine, longblocks);
855 
856  /* Undo time-freq changes that we did earlier */
857  N_B = N_B0;
858  blocks = B0;
859  for (k = 0; k < time_divide; k++) {
860  blocks >>= 1;
861  N_B <<= 1;
862  cm |= cm >> blocks;
863  celt_haar1(X, N_B, blocks);
864  }
865 
866  for (k = 0; k < recombine; k++) {
868  celt_haar1(X, N0>>k, 1<<k);
869  }
870  blocks <<= recombine;
871 
872  /* Scale output for later folding */
873  if (lowband_out) {
874  float n = sqrtf(N0);
875  for (i = 0; i < N0; i++)
876  lowband_out[i] = n * X[i];
877  }
878  cm = av_mod_uintp2(cm, blocks);
879  }
880 
881  return cm;
882 }
883 
884 static QUANT_FN(pvq_decode_band)
885 {
886 #if CONFIG_OPUS_DECODER
887  return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
888  lowband_out, level, gain, lowband_scratch, fill, 0);
889 #else
890  return 0;
891 #endif
892 }
893 
894 static QUANT_FN(pvq_encode_band)
895 {
896 #if CONFIG_OPUS_ENCODER
897  return quant_band_template(pvq, f, rc, band, X, Y, N, b, blocks, lowband, duration,
898  lowband_out, level, gain, lowband_scratch, fill, 1);
899 #else
900  return 0;
901 #endif
902 }
903 
905 {
906  CeltPVQ *s = av_malloc(sizeof(CeltPVQ));
907  if (!s)
908  return AVERROR(ENOMEM);
909 
910  s->pvq_search = ppp_pvq_search_c;
911  s->quant_band = encode ? pvq_encode_band : pvq_decode_band;
912 
913 #if CONFIG_OPUS_ENCODER && ARCH_X86
915 #endif
916 
917  *pvq = s;
918 
919  return 0;
920 }
921 
923 {
924  av_freep(pvq);
925 }
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:438
ppp_pvq_search_c
static float ppp_pvq_search_c(float *X, int *y, int K, int N)
Definition: opus_pvq.c:373
opus_ilog
#define opus_ilog(i)
Definition: opus_rc.h:31
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:488
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:468
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:450
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:884
ff_celt_pvq_init
int av_cold ff_celt_pvq_init(CeltPVQ **pvq, int encode)
Definition: opus_pvq.c:904
ff_celt_pvq_u_row
const uint32_t *const ff_celt_pvq_u_row[15]
Definition: opustab.c:1157
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
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:422
g
const char * g
Definition: vf_curves.c:117
bits
uint8_t bits
Definition: vp3data.h:141
CELT_QTHETA_OFFSET_TWOPHASE
#define CELT_QTHETA_OFFSET_TWOPHASE
Definition: opus_celt.h:50
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:922
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
X
@ X
Definition: vf_addroi.c:26
OpusRangeCoder
Definition: opus_rc.h:40
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:66
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
B
#define B
Definition: huffyuvdsp.h:32
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:478
CeltFrame
Definition: opus_celt.h:97
CELT_PVQ_V
#define CELT_PVQ_V(n, k)
Definition: opus_pvq.c:37