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_internal.h"
25 
26 #define LPC_USE_DOUBLE
27 #include "lpc.h"
28 #include "libavutil/avassert.h"
29 
30 
31 /**
32  * Apply Welch window function to audio block
33  */
34 static void lpc_apply_welch_window_c(const int32_t *data, ptrdiff_t len,
35  double *w_data)
36 {
37  int i, n2;
38  double w;
39  double c;
40 
41  if (len == 1) {
42  w_data[0] = 0.0;
43  return;
44  }
45 
46  n2 = (len >> 1);
47  c = 2.0 / (len - 1.0);
48 
49  if (len & 1) {
50  for(i=0; i<n2; i++) {
51  w = c - i - 1.0;
52  w = 1.0 - (w * w);
53  w_data[i] = data[i] * w;
54  w_data[len-1-i] = data[len-1-i] * w;
55  }
56  w_data[n2] = 0.0;
57  return;
58  }
59 
60  w_data+=n2;
61  data+=n2;
62  for(i=0; i<n2; i++) {
63  w = c - n2 + i;
64  w = 1.0 - (w * w);
65  w_data[-i-1] = data[-i-1] * w;
66  w_data[+i ] = data[+i ] * w;
67  }
68 }
69 
70 /**
71  * Calculate autocorrelation data from audio samples
72  * A Welch window function is applied before calculation.
73  */
74 static void lpc_compute_autocorr_c(const double *data, ptrdiff_t len, int lag,
75  double *autoc)
76 {
77  int i, j;
78 
79  for(j=0; j<lag; j+=2){
80  double sum0 = 1.0, sum1 = 1.0;
81  for(i=j; i<len; i++){
82  sum0 += data[i] * data[i-j];
83  sum1 += data[i] * data[i-j-1];
84  }
85  autoc[j ] = sum0;
86  autoc[j+1] = sum1;
87  }
88 
89  if(j==lag){
90  double sum = 1.0;
91  for(i=j-1; i<len; i+=2){
92  sum += data[i ] * data[i-j ]
93  + data[i+1] * data[i-j+1];
94  }
95  autoc[j] = sum;
96  }
97 }
98 
99 /**
100  * Quantize LPC coefficients
101  */
102 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
103  int32_t *lpc_out, int *shift, int min_shift,
104  int max_shift, int zero_shift)
105 {
106  int i;
107  double cmax, error;
108  int32_t qmax;
109  int sh;
110 
111  /* define maximum levels */
112  qmax = (1 << (precision - 1)) - 1;
113 
114  /* find maximum coefficient value */
115  cmax = 0.0;
116  for(i=0; i<order; i++) {
117  cmax= FFMAX(cmax, fabs(lpc_in[i]));
118  }
119 
120  /* if maximum value quantizes to zero, return all zeros */
121  if(cmax * (1 << max_shift) < 1.0) {
122  *shift = zero_shift;
123  memset(lpc_out, 0, sizeof(int32_t) * order);
124  return;
125  }
126 
127  /* calculate level shift which scales max coeff to available bits */
128  sh = max_shift;
129  while((cmax * (1 << sh) > qmax) && (sh > min_shift)) {
130  sh--;
131  }
132 
133  /* since negative shift values are unsupported in decoder, scale down
134  coefficients instead */
135  if(sh == 0 && cmax > qmax) {
136  double scale = ((double)qmax) / cmax;
137  for(i=0; i<order; i++) {
138  lpc_in[i] *= scale;
139  }
140  }
141 
142  /* output quantized coefficients and level shift */
143  error=0;
144  for(i=0; i<order; i++) {
145  error -= lpc_in[i] * (1 << sh);
146  lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
147  error -= lpc_out[i];
148  }
149  *shift = sh;
150 }
151 
152 static int estimate_best_order(double *ref, int min_order, int max_order)
153 {
154  int i, est;
155 
156  est = min_order;
157  for(i=max_order-1; i>=min_order-1; i--) {
158  if(ref[i] > 0.10) {
159  est = i+1;
160  break;
161  }
162  }
163  return est;
164 }
165 
167  const int32_t *samples, int order, double *ref)
168 {
169  double autoc[MAX_LPC_ORDER + 1];
170 
171  s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples);
172  s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc);
173  compute_ref_coefs(autoc, order, ref, NULL);
174 
175  return order;
176 }
177 
178 double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len,
179  int order, double *ref)
180 {
181  int i;
182  double signal = 0.0f, avg_err = 0.0f;
183  double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0};
184  const double a = 0.5f, b = 1.0f - a;
185 
186  /* Apply windowing */
187  for (i = 0; i <= len / 2; i++) {
188  double weight = a - b*cos((2*M_PI*i)/(len - 1));
189  s->windowed_samples[i] = weight*samples[i];
190  s->windowed_samples[len-1-i] = weight*samples[len-1-i];
191  }
192 
193  s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc);
194  signal = autoc[0];
195  compute_ref_coefs(autoc, order, ref, error);
196  for (i = 0; i < order; i++)
197  avg_err = (avg_err + error[i])/2.0f;
198  return avg_err ? signal/avg_err : NAN;
199 }
200 
201 /**
202  * Calculate LPC coefficients for multiple orders
203  *
204  * @param lpc_type LPC method for determining coefficients,
205  * see #FFLPCType for details
206  */
208  const int32_t *samples, int blocksize, int min_order,
209  int max_order, int precision,
210  int32_t coefs[][MAX_LPC_ORDER], int *shift,
211  enum FFLPCType lpc_type, int lpc_passes,
212  int omethod, int min_shift, int max_shift, int zero_shift)
213 {
214  double autoc[MAX_LPC_ORDER+1];
215  double ref[MAX_LPC_ORDER] = { 0 };
216  double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
217  int i, j, pass = 0;
218  int opt_order;
219 
220  av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER &&
221  lpc_type > FF_LPC_TYPE_FIXED);
222  av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON);
223 
224  /* reinit LPC context if parameters have changed */
225  if (blocksize != s->blocksize || max_order != s->max_order ||
226  lpc_type != s->lpc_type) {
227  ff_lpc_end(s);
228  ff_lpc_init(s, blocksize, max_order, lpc_type);
229  }
230 
231  if(lpc_passes <= 0)
232  lpc_passes = 2;
233 
234  if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) {
235  s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples);
236 
237  s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc);
238 
239  compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);
240 
241  for(i=0; i<max_order; i++)
242  ref[i] = fabs(lpc[i][i]);
243 
244  pass++;
245  }
246 
247  if (lpc_type == FF_LPC_TYPE_CHOLESKY) {
248  LLSModel *m = s->lls_models;
249  LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]);
250  double av_uninit(weight);
251  memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var));
252 
253  for(j=0; j<max_order; j++)
254  m[0].coeff[max_order-1][j] = -lpc[max_order-1][j];
255 
256  for(; pass<lpc_passes; pass++){
257  avpriv_init_lls(&m[pass&1], max_order);
258 
259  weight=0;
260  for(i=max_order; i<blocksize; i++){
261  for(j=0; j<=max_order; j++)
262  var[j]= samples[i-j];
263 
264  if(pass){
265  double eval, inv, rinv;
266  eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
267  eval= (512>>pass) + fabs(eval - var[0]);
268  inv = 1/eval;
269  rinv = sqrt(inv);
270  for(j=0; j<=max_order; j++)
271  var[j] *= rinv;
272  weight += inv;
273  }else
274  weight++;
275 
276  m[pass&1].update_lls(&m[pass&1], var);
277  }
278  avpriv_solve_lls(&m[pass&1], 0.001, 0);
279  }
280 
281  for(i=0; i<max_order; i++){
282  for(j=0; j<max_order; j++)
283  lpc[i][j]=-m[(pass-1)&1].coeff[i][j];
284  ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
285  }
286  for(i=max_order-1; i>0; i--)
287  ref[i] = ref[i-1] - ref[i];
288  }
289 
290  opt_order = max_order;
291 
292  if(omethod == ORDER_METHOD_EST) {
293  opt_order = estimate_best_order(ref, min_order, max_order);
294  i = opt_order-1;
295  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
296  min_shift, max_shift, zero_shift);
297  } else {
298  for(i=min_order-1; i<max_order; i++) {
299  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
300  min_shift, max_shift, zero_shift);
301  }
302  }
303 
304  return opt_order;
305 }
306 
307 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order,
308  enum FFLPCType lpc_type)
309 {
310  s->blocksize = blocksize;
311  s->max_order = max_order;
312  s->lpc_type = lpc_type;
313 
314  s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) *
315  sizeof(*s->windowed_samples));
316  if (!s->windowed_buffer)
317  return AVERROR(ENOMEM);
318  s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4);
319 
320  s->lpc_apply_welch_window = lpc_apply_welch_window_c;
321  s->lpc_compute_autocorr = lpc_compute_autocorr_c;
322 
323 #if ARCH_X86
325 #endif
326 
327  return 0;
328 }
329 
331 {
332  av_freep(&s->windowed_buffer);
333 }
error
static void error(const char *err)
Definition: target_bsf_fuzzer.c:31
LLSModel
Linear least squares model.
Definition: lls.h:37
FFLPCType
FFLPCType
LPC analysis type.
Definition: lpc.h:44
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
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:207
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:307
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:34
w
uint8_t w
Definition: llviddspenc.c:38
FF_LPC_TYPE_CHOLESKY
@ FF_LPC_TYPE_CHOLESKY
Cholesky factorization.
Definition: lpc.h:49
compute_lpc_coefs
static int AAC_RENAME() 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.h:167
b
#define b
Definition: input.c:41
data
const char data[16]
Definition: mxf.c:146
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
lpc.h
LPCContext
Definition: lpc.h:53
scale
static av_always_inline float scale(float x, float s)
Definition: vf_v360.c:1389
avassert.h
av_cold
#define av_cold
Definition: attributes.h:90
LOCAL_ALIGNED
#define LOCAL_ALIGNED(a, t, v,...)
Definition: mem_internal.h:114
s
#define s(width, name)
Definition: cbs_vp9.c:256
av_assert0
#define av_assert0(cond)
assert() equivalent, that is always enabled.
Definition: avassert.h:37
lls.h
NAN
#define NAN
Definition: mathematics.h:64
pass
#define pass
Definition: fft_template.c:608
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
NULL
#define NULL
Definition: coverity.c:32
ff_lpc_end
av_cold void ff_lpc_end(LPCContext *s)
Uninitialize LPCContext.
Definition: lpc.c:330
double
double
Definition: af_crystalizer.c:132
avpriv_init_lls
av_cold void avpriv_init_lls(LLSModel *m, int indep_count)
Definition: lls.c:114
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
weight
static int weight(int i, int blen, int offset)
Definition: diracdec.c:1561
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:253
MAX_LPC_ORDER
#define MAX_LPC_ORDER
Definition: lpc.h:39
MIN_LPC_ORDER
#define MIN_LPC_ORDER
Definition: lpc.h:38
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:31
M_PI
#define M_PI
Definition: mathematics.h:52
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:178
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:152
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:102
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
lrintf
#define lrintf(x)
Definition: libm_mips.h:72
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
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:74
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:264
len
int len
Definition: vorbis_enc_data.h:426
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.h:136
av_uninit
#define av_uninit(x)
Definition: attributes.h:154
LLSModel::coeff
double coeff[32][MAX_VARS]
Definition: lls.h:39
ref
static int ref[MAX_W *MAX_W]
Definition: jpeg2000dwt.c:112
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:166
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:78
FF_LPC_TYPE_LEVINSON
@ FF_LPC_TYPE_LEVINSON
Levinson-Durbin recursion.
Definition: lpc.h:48
avpriv_solve_lls
void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order)
Definition: lls.c:46
FF_LPC_TYPE_FIXED
@ FF_LPC_TYPE_FIXED
fixed LPC coefficients
Definition: lpc.h:47