FFmpeg  4.4.6
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 avg_err ? signal/avg_err : NAN;
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  /* Avoids initializing with an unused value when lpc_passes == 1 */
248  if (lpc_passes > 1)
249  for(j=0; j<max_order; j++)
250  m[0].coeff[max_order-1][j] = -lpc[max_order-1][j];
251 
252  for(; pass<lpc_passes; pass++){
253  avpriv_init_lls(&m[pass&1], max_order);
254 
255  weight=0;
256  for(i=max_order; i<blocksize; i++){
257  for(j=0; j<=max_order; j++)
258  var[j]= samples[i-j];
259 
260  if(pass){
261  double eval, inv, rinv;
262  eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
263  eval= (512>>pass) + fabs(eval - var[0]);
264  inv = 1/eval;
265  rinv = sqrt(inv);
266  for(j=0; j<=max_order; j++)
267  var[j] *= rinv;
268  weight += inv;
269  }else
270  weight++;
271 
272  m[pass&1].update_lls(&m[pass&1], var);
273  }
274  avpriv_solve_lls(&m[pass&1], 0.001, 0);
275  }
276 
277  for(i=0; i<max_order; i++){
278  for(j=0; j<max_order; j++)
279  lpc[i][j]=-m[(pass-1)&1].coeff[i][j];
280  ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
281  }
282  for(i=max_order-1; i>0; i--)
283  ref[i] = ref[i-1] - ref[i];
284  }
285 
286  opt_order = max_order;
287 
288  if(omethod == ORDER_METHOD_EST) {
289  opt_order = estimate_best_order(ref, min_order, max_order);
290  i = opt_order-1;
291  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
292  min_shift, max_shift, zero_shift);
293  } else {
294  for(i=min_order-1; i<max_order; i++) {
295  quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],
296  min_shift, max_shift, zero_shift);
297  }
298  }
299 
300  return opt_order;
301 }
302 
303 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order,
304  enum FFLPCType lpc_type)
305 {
306  s->blocksize = blocksize;
307  s->max_order = max_order;
308  s->lpc_type = lpc_type;
309 
310  s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) *
311  sizeof(*s->windowed_samples));
312  if (!s->windowed_buffer)
313  return AVERROR(ENOMEM);
314  s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4);
315 
316  s->lpc_apply_welch_window = lpc_apply_welch_window_c;
317  s->lpc_compute_autocorr = lpc_compute_autocorr_c;
318 
319  if (ARCH_X86)
321 
322  return 0;
323 }
324 
326 {
327  av_freep(&s->windowed_buffer);
328 }
#define av_uninit(x)
Definition: attributes.h:154
#define av_cold
Definition: attributes.h:88
int32_t
simple assert() macros that are a bit more flexible than ISO C assert().
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
#define av_assert0(cond)
assert() equivalent, that is always enabled.
Definition: avassert.h:37
#define s(width, name)
Definition: cbs_vp9.c:257
common internal and external API header
#define av_clip
Definition: common.h:122
#define FFMAX(a, b)
Definition: common.h:103
#define ARCH_X86
Definition: config.h:39
#define NULL
Definition: coverity.c:32
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
#define AVERROR(e)
Definition: error.h:43
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
int i
Definition: input.c:407
static int weight(int i, int blen, int offset)
Definition: diracdec.c:1561
#define lrintf(x)
Definition: libm_mips.h:70
av_cold void avpriv_init_lls(LLSModel *m, int indep_count)
Definition: lls.c:115
void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order)
Definition: lls.c:47
uint8_t w
Definition: llviddspenc.c:39
static int estimate_best_order(double *ref, int min_order, int max_order)
Definition: lpc.c:146
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_cold void ff_lpc_end(LPCContext *s)
Uninitialize LPCContext.
Definition: lpc.c:325
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
double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, int order, double *ref)
Definition: lpc.c:172
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
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
int ff_lpc_calc_ref_coefs(LPCContext *s, const int32_t *samples, int order, double *ref)
Definition: lpc.c:160
av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, enum FFLPCType lpc_type)
Initialize LPCContext.
Definition: lpc.c:303
#define ORDER_METHOD_EST
Definition: lpc.h:30
#define MAX_LPC_ORDER
Definition: lpc.h:38
FFLPCType
LPC analysis type.
Definition: lpc.h:43
@ FF_LPC_TYPE_CHOLESKY
Cholesky factorization.
Definition: lpc.h:48
@ FF_LPC_TYPE_FIXED
fixed LPC coefficients
Definition: lpc.h:46
@ FF_LPC_TYPE_LEVINSON
Levinson-Durbin recursion.
Definition: lpc.h:47
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 compute_ref_coefs(const LPC_TYPE *autoc, int max_order, LPC_TYPE *ref, LPC_TYPE *error)
Schur recursion.
Definition: lpc.h:135
#define MIN_LPC_ORDER
Definition: lpc.h:37
#define FFALIGN(x, a)
Definition: macros.h:48
#define NAN
Definition: mathematics.h:64
#define M_PI
Definition: mathematics.h:52
#define LOCAL_ALIGNED(a, t, v,...)
Definition: mem_internal.h:113
const char data[16]
Definition: mxf.c:142
static int shift(int a, int b)
Definition: sonic.c:82
Linear least squares model.
Definition: lls.h:38
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
double coeff[32][32]
Definition: lls.h:40
double(* evaluate_lls)(struct LLSModel *m, const double *var, int order)
Inner product of var[] and the LPC coefs.
Definition: lls.h:57
Definition: lpc.h:52
#define av_freep(p)
static void error(const char *err)
static int ref[MAX_W *MAX_W]
Definition: jpeg2000dwt.c:107
#define pass
Definition: tx_template.c:347
const char * b
Definition: vf_curves.c:118
static const double coeff[2][5]
Definition: vf_owdenoise.c:73
int len
static double c[64]
av_cold void ff_lpc_init_x86(LPCContext *c)
Definition: lpc.c:152