FFmpeg
lpc.c
Go to the documentation of this file.
1 /*
2  * LPC utility code
3  * Copyright (c) 2006 Justin Ruggles <justin.ruggles@gmail.com>
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21 
22 #include "libavutil/common.h"
23 #include "libavutil/lls.h"
24 #include "libavutil/mem.h"
25 #include "libavutil/mem_internal.h"
26 
27 #define LPC_USE_DOUBLE
28 #include "lpc.h"
29 #include "lpc_functions.h"
30 #include "libavutil/avassert.h"
31 
32 /**
33  * Schur recursion.
34  * Produces reflection coefficients from autocorrelation data.
35  */
36 static inline void compute_ref_coefs(const LPC_TYPE *autoc, int max_order,
38 {
39  LPC_TYPE err;
41 
42  for (int i = 0; i < max_order; i++)
43  gen0[i] = gen1[i] = autoc[i + 1];
44 
45  err = autoc[0];
46  ref[0] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1);
47  err += gen1[0] * ref[0];
48  if (error)
49  error[0] = err;
50  for (int i = 1; i < max_order; i++) {
51  for (int j = 0; j < max_order - i; j++) {
52  gen1[j] = gen1[j + 1] + ref[i - 1] * gen0[j];
53  gen0[j] = gen1[j + 1] * ref[i - 1] + gen0[j];
54  }
55  ref[i] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1);
56  err += gen1[0] * ref[i];
57  if (error)
58  error[i] = err;
59  }
60 }
61 
62 
63 /**
64  * Apply Welch window function to audio block
65  */
66 static void lpc_apply_welch_window_c(const int32_t *data, ptrdiff_t len,
67  double *w_data)
68 {
69  int i, n2;
70  double w;
71  double c;
72 
73  if (len == 1) {
74  w_data[0] = 0.0;
75  return;
76  }
77 
78  n2 = (len >> 1);
79  c = 2.0 / (len - 1.0);
80 
81  if (len & 1) {
82  for(i=0; i<n2; i++) {
83  w = c - i - 1.0;
84  w = 1.0 - (w * w);
85  w_data[i] = data[i] * w;
86  w_data[len-1-i] = data[len-1-i] * w;
87  }
88  w_data[n2] = 0.0;
89  return;
90  }
91 
92  w_data+=n2;
93  data+=n2;
94  for(i=0; i<n2; i++) {
95  w = c - n2 + i;
96  w = 1.0 - (w * w);
97  w_data[-i-1] = data[-i-1] * w;
98  w_data[+i ] = data[+i ] * w;
99  }
100 }
101 
102 /**
103  * Calculate autocorrelation data from audio samples
104  * A Welch window function is applied before calculation.
105  */
106 static void lpc_compute_autocorr_c(const double *data, ptrdiff_t len, int lag,
107  double *autoc)
108 {
109  int i, j;
110 
111  for(j=0; j<lag; j+=2){
112  double sum0 = 1.0, sum1 = 1.0;
113  for(i=j; i<len; i++){
114  sum0 += data[i] * data[i-j];
115  sum1 += data[i] * data[i-j-1];
116  }
117  autoc[j ] = sum0;
118  autoc[j+1] = sum1;
119  }
120 
121  if(j==lag){
122  double sum = 1.0;
123  for(i=j-1; i<len; i++){
124  sum += data[i] * data[i-j];
125  }
126  autoc[j] = sum;
127  }
128 }
129 
130 /**
131  * Quantize LPC coefficients
132  */
133 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
134  int32_t *lpc_out, int *shift, int min_shift,
135  int max_shift, int zero_shift)
136 {
137  int i;
138  double cmax, error;
139  int32_t qmax;
140  int sh;
141 
142  /* define maximum levels */
143  qmax = (1 << (precision - 1)) - 1;
144 
145  /* find maximum coefficient value */
146  cmax = 0.0;
147  for(i=0; i<order; i++) {
148  cmax= FFMAX(cmax, fabs(lpc_in[i]));
149  }
150 
151  /* if maximum value quantizes to zero, return all zeros */
152  if(cmax * (1 << max_shift) < 1.0) {
153  *shift = zero_shift;
154  memset(lpc_out, 0, sizeof(int32_t) * order);
155  return;
156  }
157 
158  /* calculate level shift which scales max coeff to available bits */
159  sh = max_shift;
160  while((cmax * (1 << sh) > qmax) && (sh > min_shift)) {
161  sh--;
162  }
163 
164  /* since negative shift values are unsupported in decoder, scale down
165  coefficients instead */
166  if(sh == 0 && cmax > qmax) {
167  double scale = ((double)qmax) / cmax;
168  for(i=0; i<order; i++) {
169  lpc_in[i] *= scale;
170  }
171  }
172 
173  /* output quantized coefficients and level shift */
174  error=0;
175  for(i=0; i<order; i++) {
176  error -= lpc_in[i] * (1 << sh);
177  lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
178  error -= lpc_out[i];
179  }
180  *shift = sh;
181 }
182 
183 static int estimate_best_order(double *ref, int min_order, int max_order)
184 {
185  int i, est;
186 
187  est = min_order;
188  for(i=max_order-1; i>=min_order-1; i--) {
189  if(ref[i] > 0.10) {
190  est = i+1;
191  break;
192  }
193  }
194  return est;
195 }
196 
198  const int32_t *samples, int order, double *ref)
199 {
200  double autoc[MAX_LPC_ORDER + 1];
201 
202  s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples);
203  s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc);
204  compute_ref_coefs(autoc, order, ref, NULL);
205 
206  return order;
207 }
208 
209 double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len,
210  int order, double *ref)
211 {
212  int i;
213  double signal = 0.0f, avg_err = 0.0f;
214  double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0};
215  const double a = 0.5f, b = 1.0f - a;
216 
217  /* Apply windowing */
218  for (i = 0; i <= len / 2; i++) {
219  double weight = a - b*cos((2*M_PI*i)/(len - 1));
220  s->windowed_samples[i] = weight*samples[i];
221  s->windowed_samples[len-1-i] = weight*samples[len-1-i];
222  }
223 
224  s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc);
225  signal = autoc[0];
226  compute_ref_coefs(autoc, order, ref, error);
227  for (i = 0; i < order; i++)
228  avg_err = (avg_err + error[i])/2.0f;
229  return avg_err ? signal/avg_err : NAN;
230 }
231 
232 /**
233  * Calculate LPC coefficients for multiple orders
234  *
235  * @param lpc_type LPC method for determining coefficients,
236  * see #FFLPCType for details
237  */
239  const int32_t *samples, int blocksize, int min_order,
240  int max_order, int precision,
241  int32_t coefs[][MAX_LPC_ORDER], int *shift,
242  enum FFLPCType lpc_type, int lpc_passes,
243  int omethod, int min_shift, int max_shift, int zero_shift)
244 {
245  double autoc[MAX_LPC_ORDER+1];
246  double ref[MAX_LPC_ORDER] = { 0 };
247  double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
248  int i, j, pass = 0;
249  int opt_order;
250 
251  av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER &&
252  lpc_type > FF_LPC_TYPE_FIXED);
253  av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON);
254 
255  /* reinit LPC context if parameters have changed */
256  if (blocksize != s->blocksize || max_order != s->max_order ||
257  lpc_type != s->lpc_type) {
258  ff_lpc_end(s);
259  ff_lpc_init(s, blocksize, max_order, lpc_type);
260  }
261 
262  if(lpc_passes <= 0)
263  lpc_passes = 2;
264 
265  if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) {
266  s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples);
267 
268  s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc);
269 
270  compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);
271 
272  for(i=0; i<max_order; i++)
273  ref[i] = fabs(lpc[i][i]);
274 
275  pass++;
276  }
277 
278  if (lpc_type == FF_LPC_TYPE_CHOLESKY) {
279  LLSModel *m = s->lls_models;
280  LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]);
281  double av_uninit(weight);
282  memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var));
283 
284  /* Avoids initializing with an unused value when lpc_passes == 1 */
285  if (lpc_passes > 1)
286  for(j=0; j<max_order; j++)
287  m[0].coeff[max_order-1][j] = -lpc[max_order-1][j];
288 
289  for(; pass<lpc_passes; pass++){
290  avpriv_init_lls(&m[pass&1], max_order);
291 
292  weight=0;
293  for(i=max_order; i<blocksize; i++){
294  for(j=0; j<=max_order; j++)
295  var[j]= samples[i-j];
296 
297  if(pass){
298  double eval, inv, rinv;
299  eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
300  eval= (512>>pass) + fabs(eval - var[0]);
301  inv = 1/eval;
302  rinv = sqrt(inv);
303  for(j=0; j<=max_order; j++)
304  var[j] *= rinv;
305  weight += inv;
306  }else
307  weight++;
308 
309  m[pass&1].update_lls(&m[pass&1], var);
310  }
311  avpriv_solve_lls(&m[pass&1], 0.001, 0);
312  }
313 
314  for(i=0; i<max_order; i++){
315  for(j=0; j<max_order; j++)
316  lpc[i][j]=-m[(pass-1)&1].coeff[i][j];
317  ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
318  }
319  for(i=max_order-1; i>0; i--)
320  ref[i] = ref[i-1] - ref[i];
321  }
322 
323  opt_order = max_order;
324 
325  if(omethod == ORDER_METHOD_EST) {
326  opt_order = estimate_best_order(ref, min_order, max_order);
327  i = opt_order-1;
328  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
329  min_shift, max_shift, zero_shift);
330  } else {
331  for(i=min_order-1; i<max_order; i++) {
332  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
333  min_shift, max_shift, zero_shift);
334  }
335  }
336 
337  return opt_order;
338 }
339 
340 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order,
341  enum FFLPCType lpc_type)
342 {
343  s->blocksize = blocksize;
344  s->max_order = max_order;
345  s->lpc_type = lpc_type;
346 
347  s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) *
348  sizeof(*s->windowed_samples));
349  if (!s->windowed_buffer)
350  return AVERROR(ENOMEM);
351  s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4);
352 
353  s->lpc_apply_welch_window = lpc_apply_welch_window_c;
354  s->lpc_compute_autocorr = lpc_compute_autocorr_c;
355 
356 #if ARCH_RISCV
358 #elif ARCH_X86
360 #endif
361 
362  return 0;
363 }
364 
366 {
367  av_freep(&s->windowed_buffer);
368 }
error
static void error(const char *err)
Definition: target_bsf_fuzzer.c:32
avpriv_solve_lls
void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order)
Definition: lls.c:47
LLSModel
Linear least squares model.
Definition: lls.h:37
FFLPCType
FFLPCType
LPC analysis type.
Definition: lpc.h:42
av_clip
#define av_clip
Definition: common.h:100
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
mem_internal.h
ff_lpc_calc_coefs
int ff_lpc_calc_coefs(LPCContext *s, const int32_t *samples, int blocksize, int min_order, int max_order, int precision, int32_t coefs[][MAX_LPC_ORDER], int *shift, enum FFLPCType lpc_type, int lpc_passes, int omethod, int min_shift, int max_shift, int zero_shift)
Calculate LPC coefficients for multiple orders.
Definition: lpc.c:238
ff_lpc_init
av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, enum FFLPCType lpc_type)
Initialize LPCContext.
Definition: lpc.c:340
lpc_apply_welch_window_c
static void lpc_apply_welch_window_c(const int32_t *data, ptrdiff_t len, double *w_data)
Apply Welch window function to audio block.
Definition: lpc.c:66
w
uint8_t w
Definition: llviddspenc.c:38
FF_LPC_TYPE_CHOLESKY
@ FF_LPC_TYPE_CHOLESKY
Cholesky factorization.
Definition: lpc.h:47
b
#define b
Definition: input.c:41
data
const char data[16]
Definition: mxf.c:149
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
lpc.h
LPCContext
Definition: lpc.h:51
weight
const h264_weight_func weight
Definition: h264dsp_init.c:33
avassert.h
av_cold
#define av_cold
Definition: attributes.h:90
LOCAL_ALIGNED
#define LOCAL_ALIGNED(a, t, v,...)
Definition: mem_internal.h:122
compute_ref_coefs
static void compute_ref_coefs(const LPC_TYPE *autoc, int max_order, LPC_TYPE *ref, LPC_TYPE *error)
Schur recursion.
Definition: lpc.c:36
s
#define s(width, name)
Definition: cbs_vp9.c:198
av_assert0
#define av_assert0(cond)
assert() equivalent, that is always enabled.
Definition: avassert.h:40
lls.h
NAN
#define NAN
Definition: mathematics.h:115
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
compute_lpc_coefs
static int compute_lpc_coefs(const LPC_TYPE *autoc, int max_order, LPC_TYPE *lpc, int lpc_stride, int fail, int normalize)
Levinson-Durbin recursion.
Definition: lpc_functions.h:54
NULL
#define NULL
Definition: coverity.c:32
ff_lpc_end
av_cold void ff_lpc_end(LPCContext *s)
Uninitialize LPCContext.
Definition: lpc.c:365
double
double
Definition: af_crystalizer.c:132
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
LLSModel::update_lls
void(* update_lls)(struct LLSModel *m, const double *var)
Take the outer-product of var[] with itself, and add to the covariance matrix.
Definition: lls.h:49
shift
static int shift(int a, int b)
Definition: bonk.c:261
MAX_LPC_ORDER
#define MAX_LPC_ORDER
Definition: lpc.h:37
MIN_LPC_ORDER
#define MIN_LPC_ORDER
Definition: lpc.h:36
ff_lpc_init_x86
void ff_lpc_init_x86(LPCContext *s)
Definition: lpc_init.c:106
a
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:41
ORDER_METHOD_EST
#define ORDER_METHOD_EST
Definition: lpc.h:29
M_PI
#define M_PI
Definition: mathematics.h:67
ff_lpc_calc_ref_coefs_f
double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, int order, double *ref)
Definition: lpc.c:209
LLSModel::evaluate_lls
double(* evaluate_lls)(struct LLSModel *m, const double *var, int order)
Inner product of var[] and the LPC coefs.
Definition: lls.h:56
estimate_best_order
static int estimate_best_order(double *ref, int min_order, int max_order)
Definition: lpc.c:183
quantize_lpc_coefs
static void quantize_lpc_coefs(double *lpc_in, int order, int precision, int32_t *lpc_out, int *shift, int min_shift, int max_shift, int zero_shift)
Quantize LPC coefficients.
Definition: lpc.c:133
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:67
lrintf
#define lrintf(x)
Definition: libm_mips.h:72
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
common.h
lpc_compute_autocorr_c
static void lpc_compute_autocorr_c(const double *data, ptrdiff_t len, int lag, double *autoc)
Calculate autocorrelation data from audio samples A Welch window function is applied before calculati...
Definition: lpc.c:106
av_mallocz
void * av_mallocz(size_t size)
Allocate a memory block with alignment suitable for all memory accesses (including vectors if availab...
Definition: mem.c:256
len
int len
Definition: vorbis_enc_data.h:426
LPC_TYPE
float LPC_TYPE
Definition: lpc_functions.h:45
avpriv_init_lls
av_cold void avpriv_init_lls(LLSModel *m, int indep_count)
Definition: lls.c:109
av_uninit
#define av_uninit(x)
Definition: attributes.h:154
ff_lpc_init_riscv
void ff_lpc_init_riscv(LPCContext *s)
Definition: lpc_init.c:31
LLSModel::coeff
double coeff[32][MAX_VARS]
Definition: lls.h:39
lpc_functions.h
ref
static int ref[MAX_W *MAX_W]
Definition: jpeg2000dwt.c:117
samples
Filter the word “frame” indicates either a video frame or a group of audio samples
Definition: filter_design.txt:8
ff_lpc_calc_ref_coefs
int ff_lpc_calc_ref_coefs(LPCContext *s, const int32_t *samples, int order, double *ref)
Definition: lpc.c:197
mem.h
scale
static void scale(int *out, const int *in, const int w, const int h, const int shift)
Definition: intra.c:291
FFALIGN
#define FFALIGN(x, a)
Definition: macros.h:78
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:34
int32_t
int32_t
Definition: audioconvert.c:56
coeff
static const double coeff[2][5]
Definition: vf_owdenoise.c:80
FF_LPC_TYPE_LEVINSON
@ FF_LPC_TYPE_LEVINSON
Levinson-Durbin recursion.
Definition: lpc.h:46
FF_LPC_TYPE_FIXED
@ FF_LPC_TYPE_FIXED
fixed LPC coefficients
Definition: lpc.h:45
LPC_USE_FIXED
#define LPC_USE_FIXED
Definition: lpc_functions.h:28