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 | 4359 | 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 4359 times.
|
4359 | if (len == 1) { |
74 | ✗ | w_data[0] = 0.0; | |
75 | ✗ | return; | |
76 | } | ||
77 | |||
78 | 4359 | n2 = (len >> 1); | |
79 | 4359 | c = 2.0 / (len - 1.0); | |
80 | |||
81 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 4352 times.
|
4359 | if (len & 1) { |
82 |
2/2✓ Branch 0 taken 5099 times.
✓ Branch 1 taken 7 times.
|
5106 | for(i=0; i<n2; i++) { |
83 | 5099 | w = c - i - 1.0; | |
84 | 5099 | w = 1.0 - (w * w); | |
85 | 5099 | w_data[i] = data[i] * w; | |
86 | 5099 | w_data[len-1-i] = data[len-1-i] * w; | |
87 | } | ||
88 | 7 | w_data[n2] = 0.0; | |
89 | 7 | return; | |
90 | } | ||
91 | |||
92 | 4352 | w_data+=n2; | |
93 | 4352 | data+=n2; | |
94 |
2/2✓ Branch 0 taken 10094977 times.
✓ Branch 1 taken 4352 times.
|
10099329 | for(i=0; i<n2; i++) { |
95 | 10094977 | w = c - n2 + i; | |
96 | 10094977 | w = 1.0 - (w * w); | |
97 | 10094977 | w_data[-i-1] = data[-i-1] * w; | |
98 | 10094977 | 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 | 6500 | 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 41915 times.
✓ Branch 1 taken 6500 times.
|
48415 | for(j=0; j<lag; j+=2){ |
112 | 41915 | double sum0 = 1.0, sum1 = 1.0; | |
113 |
2/2✓ Branch 0 taken 142066738 times.
✓ Branch 1 taken 41915 times.
|
142108653 | for(i=j; i<len; i++){ |
114 | 142066738 | sum0 += data[i] * data[i-j]; | |
115 | 142066738 | sum1 += data[i] * data[i-j-1]; | |
116 | } | ||
117 | 41915 | autoc[j ] = sum0; | |
118 | 41915 | autoc[j+1] = sum1; | |
119 | } | ||
120 | |||
121 |
2/2✓ Branch 0 taken 6204 times.
✓ Branch 1 taken 296 times.
|
6500 | if(j==lag){ |
122 | 6204 | double sum = 1.0; | |
123 |
2/2✓ Branch 0 taken 10645539 times.
✓ Branch 1 taken 6204 times.
|
10651743 | for(i=j-1; i<len; i+=2){ |
124 | 10645539 | sum += data[i ] * data[i-j ] | |
125 | 10645539 | + data[i+1] * data[i-j+1]; | |
126 | } | ||
127 | 6204 | autoc[j] = sum; | |
128 | } | ||
129 | 6500 | } | |
130 | |||
131 | /** | ||
132 | * Quantize LPC coefficients | ||
133 | */ | ||
134 | 12149 | static void quantize_lpc_coefs(double *lpc_in, int order, int precision, | |
135 | int32_t *lpc_out, int *shift, int min_shift, | ||
136 | int max_shift, int zero_shift) | ||
137 | { | ||
138 | int i; | ||
139 | double cmax, error; | ||
140 | int32_t qmax; | ||
141 | int sh; | ||
142 | |||
143 | /* define maximum levels */ | ||
144 | 12149 | qmax = (1 << (precision - 1)) - 1; | |
145 | |||
146 | /* find maximum coefficient value */ | ||
147 | 12149 | cmax = 0.0; | |
148 |
2/2✓ Branch 0 taken 60925 times.
✓ Branch 1 taken 12149 times.
|
73074 | for(i=0; i<order; i++) { |
149 |
2/2✓ Branch 0 taken 42903 times.
✓ Branch 1 taken 18022 times.
|
60925 | cmax= FFMAX(cmax, fabs(lpc_in[i])); |
150 | } | ||
151 | |||
152 | /* if maximum value quantizes to zero, return all zeros */ | ||
153 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12149 times.
|
12149 | if(cmax * (1 << max_shift) < 1.0) { |
154 | ✗ | *shift = zero_shift; | |
155 | ✗ | memset(lpc_out, 0, sizeof(int32_t) * order); | |
156 | ✗ | return; | |
157 | } | ||
158 | |||
159 | /* calculate level shift which scales max coeff to available bits */ | ||
160 | 12149 | sh = max_shift; | |
161 |
3/4✓ Branch 0 taken 21467 times.
✓ Branch 1 taken 12149 times.
✓ Branch 2 taken 21467 times.
✗ Branch 3 not taken.
|
33616 | while((cmax * (1 << sh) > qmax) && (sh > min_shift)) { |
162 | 21467 | sh--; | |
163 | } | ||
164 | |||
165 | /* since negative shift values are unsupported in decoder, scale down | ||
166 | coefficients instead */ | ||
167 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 12149 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12149 | if(sh == 0 && cmax > qmax) { |
168 | ✗ | double scale = ((double)qmax) / cmax; | |
169 | ✗ | for(i=0; i<order; i++) { | |
170 | ✗ | lpc_in[i] *= scale; | |
171 | } | ||
172 | } | ||
173 | |||
174 | /* output quantized coefficients and level shift */ | ||
175 | 12149 | error=0; | |
176 |
2/2✓ Branch 0 taken 60925 times.
✓ Branch 1 taken 12149 times.
|
73074 | for(i=0; i<order; i++) { |
177 | 60925 | error -= lpc_in[i] * (1 << sh); | |
178 | 60925 | lpc_out[i] = av_clip(lrintf(error), -qmax, qmax); | |
179 | 60925 | error -= lpc_out[i]; | |
180 | } | ||
181 | 12149 | *shift = sh; | |
182 | } | ||
183 | |||
184 | 9317 | static int estimate_best_order(double *ref, int min_order, int max_order) | |
185 | { | ||
186 | int i, est; | ||
187 | |||
188 | 9317 | est = min_order; | |
189 |
2/2✓ Branch 0 taken 64818 times.
✓ Branch 1 taken 174 times.
|
64992 | for(i=max_order-1; i>=min_order-1; i--) { |
190 |
2/2✓ Branch 0 taken 9143 times.
✓ Branch 1 taken 55675 times.
|
64818 | if(ref[i] > 0.10) { |
191 | 9143 | est = i+1; | |
192 | 9143 | break; | |
193 | } | ||
194 | } | ||
195 | 9317 | return est; | |
196 | } | ||
197 | |||
198 | ✗ | int ff_lpc_calc_ref_coefs(LPCContext *s, | |
199 | const int32_t *samples, int order, double *ref) | ||
200 | { | ||
201 | double autoc[MAX_LPC_ORDER + 1]; | ||
202 | |||
203 | ✗ | s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples); | |
204 | ✗ | s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc); | |
205 | ✗ | compute_ref_coefs(autoc, order, ref, NULL); | |
206 | |||
207 | ✗ | return order; | |
208 | } | ||
209 | |||
210 | 2147 | double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, | |
211 | int order, double *ref) | ||
212 | { | ||
213 | int i; | ||
214 | 2147 | double signal = 0.0f, avg_err = 0.0f; | |
215 | 2147 | double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0}; | |
216 | 2147 | const double a = 0.5f, b = 1.0f - a; | |
217 | |||
218 | /* Apply windowing */ | ||
219 |
2/2✓ Branch 0 taken 600615 times.
✓ Branch 1 taken 2147 times.
|
602762 | for (i = 0; i <= len / 2; i++) { |
220 | 600615 | double weight = a - b*cos((2*M_PI*i)/(len - 1)); | |
221 | 600615 | s->windowed_samples[i] = weight*samples[i]; | |
222 | 600615 | s->windowed_samples[len-1-i] = weight*samples[len-1-i]; | |
223 | } | ||
224 | |||
225 | 2147 | s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc); | |
226 | 2147 | signal = autoc[0]; | |
227 | 2147 | compute_ref_coefs(autoc, order, ref, error); | |
228 |
2/2✓ Branch 0 taken 24284 times.
✓ Branch 1 taken 2147 times.
|
26431 | for (i = 0; i < order; i++) |
229 | 24284 | avg_err = (avg_err + error[i])/2.0f; | |
230 |
2/2✓ Branch 0 taken 2114 times.
✓ Branch 1 taken 33 times.
|
2147 | return avg_err ? signal/avg_err : NAN; |
231 | } | ||
232 | |||
233 | /** | ||
234 | * Calculate LPC coefficients for multiple orders | ||
235 | * | ||
236 | * @param lpc_type LPC method for determining coefficients, | ||
237 | * see #FFLPCType for details | ||
238 | */ | ||
239 | 9553 | int ff_lpc_calc_coefs(LPCContext *s, | |
240 | const int32_t *samples, int blocksize, int min_order, | ||
241 | int max_order, int precision, | ||
242 | int32_t coefs[][MAX_LPC_ORDER], int *shift, | ||
243 | enum FFLPCType lpc_type, int lpc_passes, | ||
244 | int omethod, int min_shift, int max_shift, int zero_shift) | ||
245 | { | ||
246 | double autoc[MAX_LPC_ORDER+1]; | ||
247 | 9553 | double ref[MAX_LPC_ORDER] = { 0 }; | |
248 | double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; | ||
249 | 9553 | int i, j, pass = 0; | |
250 | int opt_order; | ||
251 | |||
252 | av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && | ||
253 | lpc_type > FF_LPC_TYPE_FIXED); | ||
254 |
3/4✓ Branch 0 taken 9371 times.
✓ Branch 1 taken 182 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9371 times.
|
9553 | av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON); |
255 | |||
256 | /* reinit LPC context if parameters have changed */ | ||
257 |
3/4✓ Branch 0 taken 9502 times.
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 9502 times.
✗ Branch 3 not taken.
|
9553 | if (blocksize != s->blocksize || max_order != s->max_order || |
258 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 9501 times.
|
9502 | lpc_type != s->lpc_type) { |
259 | 52 | ff_lpc_end(s); | |
260 | 52 | ff_lpc_init(s, blocksize, max_order, lpc_type); | |
261 | } | ||
262 | |||
263 |
2/2✓ Branch 0 taken 2589 times.
✓ Branch 1 taken 6964 times.
|
9553 | if(lpc_passes <= 0) |
264 | 2589 | lpc_passes = 2; | |
265 | |||
266 |
4/6✓ Branch 0 taken 182 times.
✓ Branch 1 taken 9371 times.
✓ Branch 2 taken 182 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 182 times.
✗ Branch 5 not taken.
|
9553 | if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) { |
267 | 9553 | s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples); | |
268 | |||
269 | 9553 | s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc); | |
270 | |||
271 | 9553 | compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); | |
272 | |||
273 |
2/2✓ Branch 0 taken 100850 times.
✓ Branch 1 taken 9553 times.
|
110403 | for(i=0; i<max_order; i++) |
274 | 100850 | ref[i] = fabs(lpc[i][i]); | |
275 | |||
276 | 9553 | pass++; | |
277 | } | ||
278 | |||
279 |
2/2✓ Branch 0 taken 182 times.
✓ Branch 1 taken 9371 times.
|
9553 | if (lpc_type == FF_LPC_TYPE_CHOLESKY) { |
280 | 182 | LLSModel *m = s->lls_models; | |
281 | 182 | LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]); | |
282 | 182 | double av_uninit(weight); | |
283 | 182 | memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var)); | |
284 | |||
285 |
2/2✓ Branch 0 taken 1456 times.
✓ Branch 1 taken 182 times.
|
1638 | for(j=0; j<max_order; j++) |
286 | 1456 | m[0].coeff[max_order-1][j] = -lpc[max_order-1][j]; | |
287 | |||
288 |
2/2✓ Branch 0 taken 182 times.
✓ Branch 1 taken 182 times.
|
364 | for(; pass<lpc_passes; pass++){ |
289 | 182 | avpriv_init_lls(&m[pass&1], max_order); | |
290 | |||
291 | 182 | weight=0; | |
292 |
2/2✓ Branch 0 taken 836444 times.
✓ Branch 1 taken 182 times.
|
836626 | for(i=max_order; i<blocksize; i++){ |
293 |
2/2✓ Branch 0 taken 7527996 times.
✓ Branch 1 taken 836444 times.
|
8364440 | for(j=0; j<=max_order; j++) |
294 | 7527996 | var[j]= samples[i-j]; | |
295 | |||
296 |
1/2✓ Branch 0 taken 836444 times.
✗ Branch 1 not taken.
|
836444 | if(pass){ |
297 | double eval, inv, rinv; | ||
298 | 836444 | eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); | |
299 | 836444 | eval= (512>>pass) + fabs(eval - var[0]); | |
300 | 836444 | inv = 1/eval; | |
301 | 836444 | rinv = sqrt(inv); | |
302 |
2/2✓ Branch 0 taken 7527996 times.
✓ Branch 1 taken 836444 times.
|
8364440 | for(j=0; j<=max_order; j++) |
303 | 7527996 | var[j] *= rinv; | |
304 | 836444 | weight += inv; | |
305 | }else | ||
306 | ✗ | weight++; | |
307 | |||
308 | 836444 | m[pass&1].update_lls(&m[pass&1], var); | |
309 | } | ||
310 | 182 | avpriv_solve_lls(&m[pass&1], 0.001, 0); | |
311 | } | ||
312 | |||
313 |
2/2✓ Branch 0 taken 1456 times.
✓ Branch 1 taken 182 times.
|
1638 | for(i=0; i<max_order; i++){ |
314 |
2/2✓ Branch 0 taken 11648 times.
✓ Branch 1 taken 1456 times.
|
13104 | for(j=0; j<max_order; j++) |
315 | 11648 | lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; | |
316 | 1456 | ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; | |
317 | } | ||
318 |
2/2✓ Branch 0 taken 1274 times.
✓ Branch 1 taken 182 times.
|
1456 | for(i=max_order-1; i>0; i--) |
319 | 1274 | ref[i] = ref[i-1] - ref[i]; | |
320 | } | ||
321 | |||
322 | 9553 | opt_order = max_order; | |
323 | |||
324 |
2/2✓ Branch 0 taken 9317 times.
✓ Branch 1 taken 236 times.
|
9553 | if(omethod == ORDER_METHOD_EST) { |
325 | 9317 | opt_order = estimate_best_order(ref, min_order, max_order); | |
326 | 9317 | i = opt_order-1; | |
327 | 9317 | quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], | |
328 | min_shift, max_shift, zero_shift); | ||
329 | } else { | ||
330 |
2/2✓ Branch 0 taken 2832 times.
✓ Branch 1 taken 236 times.
|
3068 | for(i=min_order-1; i<max_order; i++) { |
331 | 2832 | quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], | |
332 | min_shift, max_shift, zero_shift); | ||
333 | } | ||
334 | } | ||
335 | |||
336 | 9553 | return opt_order; | |
337 | } | ||
338 | |||
339 | 155 | av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, | |
340 | enum FFLPCType lpc_type) | ||
341 | { | ||
342 | 155 | s->blocksize = blocksize; | |
343 | 155 | s->max_order = max_order; | |
344 | 155 | s->lpc_type = lpc_type; | |
345 | |||
346 | 155 | s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) * | |
347 | sizeof(*s->windowed_samples)); | ||
348 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 155 times.
|
155 | if (!s->windowed_buffer) |
349 | ✗ | return AVERROR(ENOMEM); | |
350 | 155 | s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4); | |
351 | |||
352 | 155 | s->lpc_apply_welch_window = lpc_apply_welch_window_c; | |
353 | 155 | s->lpc_compute_autocorr = lpc_compute_autocorr_c; | |
354 | |||
355 | #if ARCH_RISCV | ||
356 | ff_lpc_init_riscv(s); | ||
357 | #elif ARCH_X86 | ||
358 | 155 | ff_lpc_init_x86(s); | |
359 | #endif | ||
360 | |||
361 | 155 | return 0; | |
362 | } | ||
363 | |||
364 | 155 | av_cold void ff_lpc_end(LPCContext *s) | |
365 | { | ||
366 | 155 | av_freep(&s->windowed_buffer); | |
367 | 155 | } | |
368 |