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