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