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 
165  s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples);
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 
314  s->lpc_apply_welch_window = lpc_apply_welch_window_c;
315  s->lpc_compute_autocorr = lpc_compute_autocorr_c;
316 
317  if (ARCH_X86)
319 
320  return 0;
321 }
322 
324 {
325  av_freep(&s->windowed_buffer);
326 }
error
static void error(const char *err)
Definition: target_bsf_fuzzer.c:30
LLSModel
Linear least squares model.
Definition: lls.h:38
FFLPCType
FFLPCType
LPC analysis type.
Definition: lpc.h:43
av_clip
#define av_clip
Definition: common.h:122
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
w
uint8_t w
Definition: llviddspenc.c:39
FF_LPC_TYPE_CHOLESKY
@ FF_LPC_TYPE_CHOLESKY
Cholesky factorization.
Definition: lpc.h:48
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:166
b
#define b
Definition: input.c:41
data
const char data[16]
Definition: mxf.c:142
estimate_best_order
static int estimate_best_order(double *ref, int min_order, int max_order)
Definition: lpc.c:146
lpc.h
lpc_apply_welch_window_c
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
LPCContext
Definition: lpc.h:52
avassert.h
av_cold
#define av_cold
Definition: attributes.h:90
LOCAL_ALIGNED
#define LOCAL_ALIGNED(a, t, v,...)
Definition: mem_internal.h:113
s
#define s(width, name)
Definition: cbs_vp9.c:257
av_assert0
#define av_assert0(cond)
assert() equivalent, that is always enabled.
Definition: avassert.h:37
lls.h
pass
#define pass
Definition: fft_template.c:603
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
NULL
#define NULL
Definition: coverity.c:32
avpriv_init_lls
av_cold void avpriv_init_lls(LLSModel *m, int indep_count)
Definition: lls.c:115
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:201
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
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:172
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:50
MAX_LPC_ORDER
#define MAX_LPC_ORDER
Definition: lpc.h:38
FFMAX
#define FFMAX(a, b)
Definition: common.h:103
MIN_LPC_ORDER
#define MIN_LPC_ORDER
Definition: lpc.h:37
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:30
M_PI
#define M_PI
Definition: mathematics.h:52
ff_lpc_calc_ref_coefs
int ff_lpc_calc_ref_coefs(LPCContext *s, const int32_t *samples, int order, double *ref)
Definition: lpc.c:160
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:57
ff_lpc_end
av_cold void ff_lpc_end(LPCContext *s)
Uninitialize LPCContext.
Definition: lpc.c:323
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
i
int i
Definition: input.c:407
lrintf
#define lrintf(x)
Definition: libm_mips.h:72
common.h
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:243
len
int len
Definition: vorbis_enc_data.h:452
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:135
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:96
av_uninit
#define av_uninit(x)
Definition: attributes.h:154
lpc_compute_autocorr_c
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
ff_lpc_init_x86
av_cold void ff_lpc_init_x86(LPCContext *c)
Definition: lpc.c:152
LLSModel::coeff
double coeff[32][MAX_VARS]
Definition: lls.h:40
ref
static int ref[MAX_W *MAX_W]
Definition: jpeg2000dwt.c:107
samples
Filter the word “frame” indicates either a video frame or a group of audio samples
Definition: filter_design.txt:8
shift
static int shift(int a, int b)
Definition: sonic.c:82
FFALIGN
#define FFALIGN(x, a)
Definition: macros.h:48
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:35
int32_t
int32_t
Definition: audioconvert.c:56
coeff
static const double coeff[2][5]
Definition: vf_owdenoise.c:73
FF_LPC_TYPE_LEVINSON
@ FF_LPC_TYPE_LEVINSON
Levinson-Durbin recursion.
Definition: lpc.h:47
avpriv_solve_lls
void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order)
Definition: lls.c:47
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:301
FF_LPC_TYPE_FIXED
@ FF_LPC_TYPE_FIXED
fixed LPC coefficients
Definition: lpc.h:46