Line | Branch | Exec | Source |
---|---|---|---|
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.h" | ||
25 | #include "libavutil/mem_internal.h" | ||
26 | |||
27 | #define LPC_USE_DOUBLE | ||
28 | #include "lpc.h" | ||
29 | #include "lpc_functions.h" | ||
30 | #include "libavutil/avassert.h" | ||
31 | |||
32 | /** | ||
33 | * Schur recursion. | ||
34 | * Produces reflection coefficients from autocorrelation data. | ||
35 | */ | ||
36 | 2147 | static inline void compute_ref_coefs(const LPC_TYPE *autoc, int max_order, | |
37 | LPC_TYPE *ref, LPC_TYPE *error) | ||
38 | { | ||
39 | LPC_TYPE err; | ||
40 | LPC_TYPE gen0[MAX_LPC_ORDER], gen1[MAX_LPC_ORDER]; | ||
41 | |||
42 |
2/2✓ Branch 0 taken 24284 times.
✓ Branch 1 taken 2147 times.
|
26431 | for (int i = 0; i < max_order; i++) |
43 | 24284 | gen0[i] = gen1[i] = autoc[i + 1]; | |
44 | |||
45 | 2147 | err = autoc[0]; | |
46 |
1/2✓ Branch 0 taken 2147 times.
✗ Branch 1 not taken.
|
2147 | ref[0] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1); |
47 | 2147 | err += gen1[0] * ref[0]; | |
48 |
1/2✓ Branch 0 taken 2147 times.
✗ Branch 1 not taken.
|
2147 | if (error) |
49 | 2147 | error[0] = err; | |
50 |
2/2✓ Branch 0 taken 22137 times.
✓ Branch 1 taken 2147 times.
|
24284 | for (int i = 1; i < max_order; i++) { |
51 |
2/2✓ Branch 0 taken 128382 times.
✓ Branch 1 taken 22137 times.
|
150519 | for (int j = 0; j < max_order - i; j++) { |
52 | 128382 | gen1[j] = gen1[j + 1] + ref[i - 1] * gen0[j]; | |
53 | 128382 | gen0[j] = gen1[j + 1] * ref[i - 1] + gen0[j]; | |
54 | } | ||
55 |
2/2✓ Branch 0 taken 21939 times.
✓ Branch 1 taken 198 times.
|
22137 | ref[i] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1); |
56 | 22137 | err += gen1[0] * ref[i]; | |
57 |
1/2✓ Branch 0 taken 22137 times.
✗ Branch 1 not taken.
|
22137 | if (error) |
58 | 22137 | error[i] = err; | |
59 | } | ||
60 | 2147 | } | |
61 | |||
62 | |||
63 | /** | ||
64 | * Apply Welch window function to audio block | ||
65 | */ | ||
66 | 4429 | static void lpc_apply_welch_window_c(const int32_t *data, ptrdiff_t len, | |
67 | double *w_data) | ||
68 | { | ||
69 | int i, n2; | ||
70 | double w; | ||
71 | double c; | ||
72 | |||
73 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4429 times.
|
4429 | if (len == 1) { |
74 | ✗ | w_data[0] = 0.0; | |
75 | ✗ | return; | |
76 | } | ||
77 | |||
78 | 4429 | n2 = (len >> 1); | |
79 | 4429 | c = 2.0 / (len - 1.0); | |
80 | |||
81 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4422 times.
|
4429 | if (len & 1) { |
82 |
2/2✓ Branch 0 taken 9400 times.
✓ Branch 1 taken 7 times.
|
9407 | for(i=0; i<n2; i++) { |
83 | 9400 | w = c - i - 1.0; | |
84 | 9400 | w = 1.0 - (w * w); | |
85 | 9400 | w_data[i] = data[i] * w; | |
86 | 9400 | w_data[len-1-i] = data[len-1-i] * w; | |
87 | } | ||
88 | 7 | w_data[n2] = 0.0; | |
89 | 7 | return; | |
90 | } | ||
91 | |||
92 | 4422 | w_data+=n2; | |
93 | 4422 | data+=n2; | |
94 |
2/2✓ Branch 0 taken 10253628 times.
✓ Branch 1 taken 4422 times.
|
10258050 | for(i=0; i<n2; i++) { |
95 | 10253628 | w = c - n2 + i; | |
96 | 10253628 | w = 1.0 - (w * w); | |
97 | 10253628 | w_data[-i-1] = data[-i-1] * w; | |
98 | 10253628 | w_data[+i ] = data[+i ] * w; | |
99 | } | ||
100 | } | ||
101 | |||
102 | /** | ||
103 | * Calculate autocorrelation data from audio samples | ||
104 | * A Welch window function is applied before calculation. | ||
105 | */ | ||
106 | 6576 | static void lpc_compute_autocorr_c(const double *data, ptrdiff_t len, int lag, | |
107 | double *autoc) | ||
108 | { | ||
109 | int i, j; | ||
110 | |||
111 |
2/2✓ Branch 0 taken 42225 times.
✓ Branch 1 taken 6576 times.
|
48801 | for(j=0; j<lag; j+=2){ |
112 | 42225 | double sum0 = 1.0, sum1 = 1.0; | |
113 |
2/2✓ Branch 0 taken 143414412 times.
✓ Branch 1 taken 42225 times.
|
143456637 | for(i=j; i<len; i++){ |
114 | 143414412 | sum0 += data[i] * data[i-j]; | |
115 | 143414412 | sum1 += data[i] * data[i-j-1]; | |
116 | } | ||
117 | 42225 | autoc[j ] = sum0; | |
118 | 42225 | autoc[j+1] = sum1; | |
119 | } | ||
120 | |||
121 |
2/2✓ Branch 0 taken 6280 times.
✓ Branch 1 taken 296 times.
|
6576 | if(j==lag){ |
122 | 6280 | double sum = 1.0; | |
123 |
2/2✓ Branch 0 taken 21615802 times.
✓ Branch 1 taken 6280 times.
|
21622082 | for(i=j-1; i<len; i++){ |
124 | 21615802 | sum += data[i] * data[i-j]; | |
125 | } | ||
126 | 6280 | autoc[j] = sum; | |
127 | } | ||
128 | 6576 | } | |
129 | |||
130 | /** | ||
131 | * Quantize LPC coefficients | ||
132 | */ | ||
133 | 12219 | static void quantize_lpc_coefs(double *lpc_in, int order, int precision, | |
134 | int32_t *lpc_out, int *shift, int min_shift, | ||
135 | int max_shift, int zero_shift) | ||
136 | { | ||
137 | int i; | ||
138 | double cmax, error; | ||
139 | int32_t qmax; | ||
140 | int sh; | ||
141 | |||
142 | /* define maximum levels */ | ||
143 | 12219 | qmax = (1 << (precision - 1)) - 1; | |
144 | |||
145 | /* find maximum coefficient value */ | ||
146 | 12219 | cmax = 0.0; | |
147 |
2/2✓ Branch 0 taken 61485 times.
✓ Branch 1 taken 12219 times.
|
73704 | for(i=0; i<order; i++) { |
148 |
2/2✓ Branch 0 taken 43393 times.
✓ Branch 1 taken 18092 times.
|
61485 | cmax= FFMAX(cmax, fabs(lpc_in[i])); |
149 | } | ||
150 | |||
151 | /* if maximum value quantizes to zero, return all zeros */ | ||
152 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12219 times.
|
12219 | if(cmax * (1 << max_shift) < 1.0) { |
153 | ✗ | *shift = zero_shift; | |
154 | ✗ | memset(lpc_out, 0, sizeof(int32_t) * order); | |
155 | ✗ | return; | |
156 | } | ||
157 | |||
158 | /* calculate level shift which scales max coeff to available bits */ | ||
159 | 12219 | sh = max_shift; | |
160 |
3/4✓ Branch 0 taken 21649 times.
✓ Branch 1 taken 12219 times.
✓ Branch 2 taken 21649 times.
✗ Branch 3 not taken.
|
33868 | while((cmax * (1 << sh) > qmax) && (sh > min_shift)) { |
161 | 21649 | sh--; | |
162 | } | ||
163 | |||
164 | /* since negative shift values are unsupported in decoder, scale down | ||
165 | coefficients instead */ | ||
166 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 12219 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12219 | if(sh == 0 && cmax > qmax) { |
167 | ✗ | double scale = ((double)qmax) / cmax; | |
168 | ✗ | for(i=0; i<order; i++) { | |
169 | ✗ | lpc_in[i] *= scale; | |
170 | } | ||
171 | } | ||
172 | |||
173 | /* output quantized coefficients and level shift */ | ||
174 | 12219 | error=0; | |
175 |
2/2✓ Branch 0 taken 61485 times.
✓ Branch 1 taken 12219 times.
|
73704 | for(i=0; i<order; i++) { |
176 | 61485 | error -= lpc_in[i] * (1 << sh); | |
177 | 61485 | lpc_out[i] = av_clip(lrintf(error), -qmax, qmax); | |
178 | 61485 | error -= lpc_out[i]; | |
179 | } | ||
180 | 12219 | *shift = sh; | |
181 | } | ||
182 | |||
183 | 9387 | static int estimate_best_order(double *ref, int min_order, int max_order) | |
184 | { | ||
185 | int i, est; | ||
186 | |||
187 | 9387 | est = min_order; | |
188 |
2/2✓ Branch 0 taken 64888 times.
✓ Branch 1 taken 174 times.
|
65062 | for(i=max_order-1; i>=min_order-1; i--) { |
189 |
2/2✓ Branch 0 taken 9213 times.
✓ Branch 1 taken 55675 times.
|
64888 | if(ref[i] > 0.10) { |
190 | 9213 | est = i+1; | |
191 | 9213 | break; | |
192 | } | ||
193 | } | ||
194 | 9387 | return est; | |
195 | } | ||
196 | |||
197 | ✗ | int ff_lpc_calc_ref_coefs(LPCContext *s, | |
198 | const int32_t *samples, int order, double *ref) | ||
199 | { | ||
200 | double autoc[MAX_LPC_ORDER + 1]; | ||
201 | |||
202 | ✗ | s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples); | |
203 | ✗ | s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc); | |
204 | ✗ | compute_ref_coefs(autoc, order, ref, NULL); | |
205 | |||
206 | ✗ | return order; | |
207 | } | ||
208 | |||
209 | 2147 | double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, | |
210 | int order, double *ref) | ||
211 | { | ||
212 | int i; | ||
213 | 2147 | double signal = 0.0f, avg_err = 0.0f; | |
214 | 2147 | double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0}; | |
215 | 2147 | const double a = 0.5f, b = 1.0f - a; | |
216 | |||
217 | /* Apply windowing */ | ||
218 |
2/2✓ Branch 0 taken 600615 times.
✓ Branch 1 taken 2147 times.
|
602762 | for (i = 0; i <= len / 2; i++) { |
219 | 600615 | double weight = a - b*cos((2*M_PI*i)/(len - 1)); | |
220 | 600615 | s->windowed_samples[i] = weight*samples[i]; | |
221 | 600615 | s->windowed_samples[len-1-i] = weight*samples[len-1-i]; | |
222 | } | ||
223 | |||
224 | 2147 | s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc); | |
225 | 2147 | signal = autoc[0]; | |
226 | 2147 | compute_ref_coefs(autoc, order, ref, error); | |
227 |
2/2✓ Branch 0 taken 24284 times.
✓ Branch 1 taken 2147 times.
|
26431 | for (i = 0; i < order; i++) |
228 | 24284 | avg_err = (avg_err + error[i])/2.0f; | |
229 |
2/2✓ Branch 0 taken 2114 times.
✓ Branch 1 taken 33 times.
|
2147 | return avg_err ? signal/avg_err : NAN; |
230 | } | ||
231 | |||
232 | /** | ||
233 | * Calculate LPC coefficients for multiple orders | ||
234 | * | ||
235 | * @param lpc_type LPC method for determining coefficients, | ||
236 | * see #FFLPCType for details | ||
237 | */ | ||
238 | 9623 | int ff_lpc_calc_coefs(LPCContext *s, | |
239 | const int32_t *samples, int blocksize, int min_order, | ||
240 | int max_order, int precision, | ||
241 | int32_t coefs[][MAX_LPC_ORDER], int *shift, | ||
242 | enum FFLPCType lpc_type, int lpc_passes, | ||
243 | int omethod, int min_shift, int max_shift, int zero_shift) | ||
244 | { | ||
245 | double autoc[MAX_LPC_ORDER+1]; | ||
246 | 9623 | double ref[MAX_LPC_ORDER] = { 0 }; | |
247 | double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; | ||
248 | 9623 | int i, j, pass = 0; | |
249 | int opt_order; | ||
250 | |||
251 | av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && | ||
252 | lpc_type > FF_LPC_TYPE_FIXED); | ||
253 |
3/4✓ Branch 0 taken 9441 times.
✓ Branch 1 taken 182 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9441 times.
|
9623 | av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON); |
254 | |||
255 | /* reinit LPC context if parameters have changed */ | ||
256 |
3/4✓ Branch 0 taken 9565 times.
✓ Branch 1 taken 58 times.
✓ Branch 2 taken 9565 times.
✗ Branch 3 not taken.
|
9623 | if (blocksize != s->blocksize || max_order != s->max_order || |
257 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 9564 times.
|
9565 | lpc_type != s->lpc_type) { |
258 | 59 | ff_lpc_end(s); | |
259 | 59 | ff_lpc_init(s, blocksize, max_order, lpc_type); | |
260 | } | ||
261 | |||
262 |
2/2✓ Branch 0 taken 2589 times.
✓ Branch 1 taken 7034 times.
|
9623 | if(lpc_passes <= 0) |
263 | 2589 | lpc_passes = 2; | |
264 | |||
265 |
4/6✓ Branch 0 taken 182 times.
✓ Branch 1 taken 9441 times.
✓ Branch 2 taken 182 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 182 times.
✗ Branch 5 not taken.
|
9623 | if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) { |
266 | 9623 | s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples); | |
267 | |||
268 | 9623 | s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc); | |
269 | |||
270 | 9623 | compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); | |
271 | |||
272 |
2/2✓ Branch 0 taken 101410 times.
✓ Branch 1 taken 9623 times.
|
111033 | for(i=0; i<max_order; i++) |
273 | 101410 | ref[i] = fabs(lpc[i][i]); | |
274 | |||
275 | 9623 | pass++; | |
276 | } | ||
277 | |||
278 |
2/2✓ Branch 0 taken 182 times.
✓ Branch 1 taken 9441 times.
|
9623 | if (lpc_type == FF_LPC_TYPE_CHOLESKY) { |
279 | 182 | LLSModel *m = s->lls_models; | |
280 | 182 | LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]); | |
281 | 182 | double av_uninit(weight); | |
282 | 182 | memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var)); | |
283 | |||
284 | /* Avoids initializing with an unused value when lpc_passes == 1 */ | ||
285 |
1/2✓ Branch 0 taken 182 times.
✗ Branch 1 not taken.
|
182 | if (lpc_passes > 1) |
286 |
2/2✓ Branch 0 taken 1456 times.
✓ Branch 1 taken 182 times.
|
1638 | for(j=0; j<max_order; j++) |
287 | 1456 | m[0].coeff[max_order-1][j] = -lpc[max_order-1][j]; | |
288 | |||
289 |
2/2✓ Branch 0 taken 182 times.
✓ Branch 1 taken 182 times.
|
364 | for(; pass<lpc_passes; pass++){ |
290 | 182 | avpriv_init_lls(&m[pass&1], max_order); | |
291 | |||
292 | 182 | weight=0; | |
293 |
2/2✓ Branch 0 taken 836444 times.
✓ Branch 1 taken 182 times.
|
836626 | for(i=max_order; i<blocksize; i++){ |
294 |
2/2✓ Branch 0 taken 7527996 times.
✓ Branch 1 taken 836444 times.
|
8364440 | for(j=0; j<=max_order; j++) |
295 | 7527996 | var[j]= samples[i-j]; | |
296 | |||
297 |
1/2✓ Branch 0 taken 836444 times.
✗ Branch 1 not taken.
|
836444 | if(pass){ |
298 | double eval, inv, rinv; | ||
299 | 836444 | eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); | |
300 | 836444 | eval= (512>>pass) + fabs(eval - var[0]); | |
301 | 836444 | inv = 1/eval; | |
302 | 836444 | rinv = sqrt(inv); | |
303 |
2/2✓ Branch 0 taken 7527996 times.
✓ Branch 1 taken 836444 times.
|
8364440 | for(j=0; j<=max_order; j++) |
304 | 7527996 | var[j] *= rinv; | |
305 | 836444 | weight += inv; | |
306 | }else | ||
307 | ✗ | weight++; | |
308 | |||
309 | 836444 | m[pass&1].update_lls(&m[pass&1], var); | |
310 | } | ||
311 | 182 | avpriv_solve_lls(&m[pass&1], 0.001, 0); | |
312 | } | ||
313 | |||
314 |
2/2✓ Branch 0 taken 1456 times.
✓ Branch 1 taken 182 times.
|
1638 | for(i=0; i<max_order; i++){ |
315 |
2/2✓ Branch 0 taken 11648 times.
✓ Branch 1 taken 1456 times.
|
13104 | for(j=0; j<max_order; j++) |
316 | 11648 | lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; | |
317 | 1456 | ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; | |
318 | } | ||
319 |
2/2✓ Branch 0 taken 1274 times.
✓ Branch 1 taken 182 times.
|
1456 | for(i=max_order-1; i>0; i--) |
320 | 1274 | ref[i] = ref[i-1] - ref[i]; | |
321 | } | ||
322 | |||
323 | 9623 | opt_order = max_order; | |
324 | |||
325 |
2/2✓ Branch 0 taken 9387 times.
✓ Branch 1 taken 236 times.
|
9623 | if(omethod == ORDER_METHOD_EST) { |
326 | 9387 | opt_order = estimate_best_order(ref, min_order, max_order); | |
327 | 9387 | i = opt_order-1; | |
328 | 9387 | quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], | |
329 | min_shift, max_shift, zero_shift); | ||
330 | } else { | ||
331 |
2/2✓ Branch 0 taken 2832 times.
✓ Branch 1 taken 236 times.
|
3068 | for(i=min_order-1; i<max_order; i++) { |
332 | 2832 | quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], | |
333 | min_shift, max_shift, zero_shift); | ||
334 | } | ||
335 | } | ||
336 | |||
337 | 9623 | return opt_order; | |
338 | } | ||
339 | |||
340 | 195 | av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, | |
341 | enum FFLPCType lpc_type) | ||
342 | { | ||
343 | 195 | s->blocksize = blocksize; | |
344 | 195 | s->max_order = max_order; | |
345 | 195 | s->lpc_type = lpc_type; | |
346 | |||
347 | 195 | s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) * | |
348 | sizeof(*s->windowed_samples)); | ||
349 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 195 times.
|
195 | if (!s->windowed_buffer) |
350 | ✗ | return AVERROR(ENOMEM); | |
351 | 195 | s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4); | |
352 | |||
353 | 195 | s->lpc_apply_welch_window = lpc_apply_welch_window_c; | |
354 | 195 | s->lpc_compute_autocorr = lpc_compute_autocorr_c; | |
355 | |||
356 | #if ARCH_RISCV | ||
357 | ff_lpc_init_riscv(s); | ||
358 | #elif ARCH_X86 | ||
359 | 195 | ff_lpc_init_x86(s); | |
360 | #endif | ||
361 | |||
362 | 195 | return 0; | |
363 | } | ||
364 | |||
365 | 195 | av_cold void ff_lpc_end(LPCContext *s) | |
366 | { | ||
367 | 195 | av_freep(&s->windowed_buffer); | |
368 | 195 | } | |
369 |