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 |
|
3969 |
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 |
|
3969 |
n2 = (len >> 1); |
42 |
|
3969 |
c = 2.0 / (len - 1.0); |
43 |
|
|
|
44 |
✓✓ |
3969 |
if (len & 1) { |
45 |
✓✓ |
1050 |
for(i=0; i<n2; i++) { |
46 |
|
1048 |
w = c - i - 1.0; |
47 |
|
1048 |
w = 1.0 - (w * w); |
48 |
|
1048 |
w_data[i] = data[i] * w; |
49 |
|
1048 |
w_data[len-1-i] = data[len-1-i] * w; |
50 |
|
|
} |
51 |
|
2 |
return; |
52 |
|
|
} |
53 |
|
|
|
54 |
|
3967 |
w_data+=n2; |
55 |
|
3967 |
data+=n2; |
56 |
✓✓ |
9253857 |
for(i=0; i<n2; i++) { |
57 |
|
9249890 |
w = c - n2 + i; |
58 |
|
9249890 |
w = 1.0 - (w * w); |
59 |
|
9249890 |
w_data[-i-1] = data[-i-1] * w; |
60 |
|
9249890 |
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 |
|
7733 |
static void lpc_compute_autocorr_c(const double *data, int len, int lag, |
69 |
|
|
double *autoc) |
70 |
|
|
{ |
71 |
|
|
int i, j; |
72 |
|
|
|
73 |
✓✓ |
57574 |
for(j=0; j<lag; j+=2){ |
74 |
|
49841 |
double sum0 = 1.0, sum1 = 1.0; |
75 |
✓✓ |
141092269 |
for(i=j; i<len; i++){ |
76 |
|
141042428 |
sum0 += data[i] * data[i-j]; |
77 |
|
141042428 |
sum1 += data[i] * data[i-j-1]; |
78 |
|
|
} |
79 |
|
49841 |
autoc[j ] = sum0; |
80 |
|
49841 |
autoc[j+1] = sum1; |
81 |
|
|
} |
82 |
|
|
|
83 |
✓✓ |
7733 |
if(j==lag){ |
84 |
|
7317 |
double sum = 1.0; |
85 |
✓✓ |
10279207 |
for(i=j-1; i<len; i+=2){ |
86 |
|
10271890 |
sum += data[i ] * data[i-j ] |
87 |
|
10271890 |
+ data[i+1] * data[i-j+1]; |
88 |
|
|
} |
89 |
|
7317 |
autoc[j] = sum; |
90 |
|
|
} |
91 |
|
7733 |
} |
92 |
|
|
|
93 |
|
|
/** |
94 |
|
|
* Quantize LPC coefficients |
95 |
|
|
*/ |
96 |
|
11765 |
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 |
|
11765 |
qmax = (1 << (precision - 1)) - 1; |
107 |
|
|
|
108 |
|
|
/* find maximum coefficient value */ |
109 |
|
11765 |
cmax = 0.0; |
110 |
✓✓ |
69630 |
for(i=0; i<order; i++) { |
111 |
✓✓ |
57865 |
cmax= FFMAX(cmax, fabs(lpc_in[i])); |
112 |
|
|
} |
113 |
|
|
|
114 |
|
|
/* if maximum value quantizes to zero, return all zeros */ |
115 |
✗✓ |
11765 |
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 |
|
11765 |
sh = max_shift; |
123 |
✓✓✓✗
|
32239 |
while((cmax * (1 << sh) > qmax) && (sh > min_shift)) { |
124 |
|
20474 |
sh--; |
125 |
|
|
} |
126 |
|
|
|
127 |
|
|
/* since negative shift values are unsupported in decoder, scale down |
128 |
|
|
coefficients instead */ |
129 |
✗✓✗✗
|
11765 |
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 |
|
11765 |
error=0; |
138 |
✓✓ |
69630 |
for(i=0; i<order; i++) { |
139 |
|
57865 |
error -= lpc_in[i] * (1 << sh); |
140 |
|
57865 |
lpc_out[i] = av_clip(lrintf(error), -qmax, qmax); |
141 |
|
57865 |
error -= lpc_out[i]; |
142 |
|
|
} |
143 |
|
11765 |
*shift = sh; |
144 |
|
|
} |
145 |
|
|
|
146 |
|
8933 |
static int estimate_best_order(double *ref, int min_order, int max_order) |
147 |
|
|
{ |
148 |
|
|
int i, est; |
149 |
|
|
|
150 |
|
8933 |
est = min_order; |
151 |
✓✓ |
64596 |
for(i=max_order-1; i>=min_order-1; i--) { |
152 |
✓✓ |
64422 |
if(ref[i] > 0.10) { |
153 |
|
8759 |
est = i+1; |
154 |
|
8759 |
break; |
155 |
|
|
} |
156 |
|
|
} |
157 |
|
8933 |
return est; |
158 |
|
|
} |
159 |
|
|
|
160 |
|
|
int ff_lpc_calc_ref_coefs(LPCContext *s, |
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 |
|
3764 |
double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, |
173 |
|
|
int order, double *ref) |
174 |
|
|
{ |
175 |
|
|
int i; |
176 |
|
3764 |
double signal = 0.0f, avg_err = 0.0f; |
177 |
|
3764 |
double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0}; |
178 |
|
3764 |
const double a = 0.5f, b = 1.0f - a; |
179 |
|
|
|
180 |
|
|
/* Apply windowing */ |
181 |
✓✓ |
1090012 |
for (i = 0; i <= len / 2; i++) { |
182 |
|
1086248 |
double weight = a - b*cos((2*M_PI*i)/(len - 1)); |
183 |
|
1086248 |
s->windowed_samples[i] = weight*samples[i]; |
184 |
|
1086248 |
s->windowed_samples[len-1-i] = weight*samples[len-1-i]; |
185 |
|
|
} |
186 |
|
|
|
187 |
|
3764 |
s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc); |
188 |
|
3764 |
signal = autoc[0]; |
189 |
|
3764 |
compute_ref_coefs(autoc, order, ref, error); |
190 |
✓✓ |
46852 |
for (i = 0; i < order; i++) |
191 |
|
43088 |
avg_err = (avg_err + error[i])/2.0f; |
192 |
|
3764 |
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 |
|
|
*/ |
201 |
|
9169 |
int ff_lpc_calc_coefs(LPCContext *s, |
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 |
|
9169 |
double ref[MAX_LPC_ORDER] = { 0 }; |
210 |
|
|
double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; |
211 |
|
9169 |
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 |
✓✓✗✓
|
9169 |
av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON); |
217 |
|
|
|
218 |
|
|
/* reinit LPC context if parameters have changed */ |
219 |
✓✓✓✗
|
9169 |
if (blocksize != s->blocksize || max_order != s->max_order || |
220 |
✓✓ |
9156 |
lpc_type != s->lpc_type) { |
221 |
|
14 |
ff_lpc_end(s); |
222 |
|
14 |
ff_lpc_init(s, blocksize, max_order, lpc_type); |
223 |
|
|
} |
224 |
|
|
|
225 |
✓✓ |
9169 |
if(lpc_passes <= 0) |
226 |
|
2589 |
lpc_passes = 2; |
227 |
|
|
|
228 |
✓✓✓✗ ✓✗ |
9169 |
if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) { |
229 |
|
9169 |
s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples); |
230 |
|
|
|
231 |
|
9169 |
s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc); |
232 |
|
|
|
233 |
|
9169 |
compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); |
234 |
|
|
|
235 |
✓✓ |
106947 |
for(i=0; i<max_order; i++) |
236 |
|
97778 |
ref[i] = fabs(lpc[i][i]); |
237 |
|
|
|
238 |
|
9169 |
pass++; |
239 |
|
|
} |
240 |
|
|
|
241 |
✓✓ |
9169 |
if (lpc_type == FF_LPC_TYPE_CHOLESKY) { |
242 |
|
182 |
LLSModel *m = s->lls_models; |
243 |
|
182 |
LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]); |
244 |
|
182 |
double av_uninit(weight); |
245 |
|
182 |
memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var)); |
246 |
|
|
|
247 |
✓✓ |
1638 |
for(j=0; j<max_order; j++) |
248 |
|
1456 |
m[0].coeff[max_order-1][j] = -lpc[max_order-1][j]; |
249 |
|
|
|
250 |
✓✓ |
364 |
for(; pass<lpc_passes; pass++){ |
251 |
|
182 |
avpriv_init_lls(&m[pass&1], max_order); |
252 |
|
|
|
253 |
|
182 |
weight=0; |
254 |
✓✓ |
836626 |
for(i=max_order; i<blocksize; i++){ |
255 |
✓✓ |
8364440 |
for(j=0; j<=max_order; j++) |
256 |
|
7527996 |
var[j]= samples[i-j]; |
257 |
|
|
|
258 |
✓✗ |
836444 |
if(pass){ |
259 |
|
|
double eval, inv, rinv; |
260 |
|
836444 |
eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); |
261 |
|
836444 |
eval= (512>>pass) + fabs(eval - var[0]); |
262 |
|
836444 |
inv = 1/eval; |
263 |
|
836444 |
rinv = sqrt(inv); |
264 |
✓✓ |
8364440 |
for(j=0; j<=max_order; j++) |
265 |
|
7527996 |
var[j] *= rinv; |
266 |
|
836444 |
weight += inv; |
267 |
|
|
}else |
268 |
|
|
weight++; |
269 |
|
|
|
270 |
|
836444 |
m[pass&1].update_lls(&m[pass&1], var); |
271 |
|
|
} |
272 |
|
182 |
avpriv_solve_lls(&m[pass&1], 0.001, 0); |
273 |
|
|
} |
274 |
|
|
|
275 |
✓✓ |
1638 |
for(i=0; i<max_order; i++){ |
276 |
✓✓ |
13104 |
for(j=0; j<max_order; j++) |
277 |
|
11648 |
lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; |
278 |
|
1456 |
ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; |
279 |
|
|
} |
280 |
✓✓ |
1456 |
for(i=max_order-1; i>0; i--) |
281 |
|
1274 |
ref[i] = ref[i-1] - ref[i]; |
282 |
|
|
} |
283 |
|
|
|
284 |
|
9169 |
opt_order = max_order; |
285 |
|
|
|
286 |
✓✓ |
9169 |
if(omethod == ORDER_METHOD_EST) { |
287 |
|
8933 |
opt_order = estimate_best_order(ref, min_order, max_order); |
288 |
|
8933 |
i = opt_order-1; |
289 |
|
8933 |
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], |
290 |
|
|
min_shift, max_shift, zero_shift); |
291 |
|
|
} else { |
292 |
✓✓ |
3068 |
for(i=min_order-1; i<max_order; i++) { |
293 |
|
2832 |
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], |
294 |
|
|
min_shift, max_shift, zero_shift); |
295 |
|
|
} |
296 |
|
|
} |
297 |
|
|
|
298 |
|
9169 |
return opt_order; |
299 |
|
|
} |
300 |
|
|
|
301 |
|
66 |
av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, |
302 |
|
|
enum FFLPCType lpc_type) |
303 |
|
|
{ |
304 |
|
66 |
s->blocksize = blocksize; |
305 |
|
66 |
s->max_order = max_order; |
306 |
|
66 |
s->lpc_type = lpc_type; |
307 |
|
|
|
308 |
|
66 |
s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) * |
309 |
|
|
sizeof(*s->windowed_samples)); |
310 |
✗✓ |
66 |
if (!s->windowed_buffer) |
311 |
|
|
return AVERROR(ENOMEM); |
312 |
|
66 |
s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4); |
313 |
|
|
|
314 |
|
66 |
s->lpc_apply_welch_window = lpc_apply_welch_window_c; |
315 |
|
66 |
s->lpc_compute_autocorr = lpc_compute_autocorr_c; |
316 |
|
|
|
317 |
|
|
if (ARCH_X86) |
318 |
|
66 |
ff_lpc_init_x86(s); |
319 |
|
|
|
320 |
|
66 |
return 0; |
321 |
|
|
} |
322 |
|
|
|
323 |
|
66 |
av_cold void ff_lpc_end(LPCContext *s) |
324 |
|
|
{ |
325 |
|
66 |
av_freep(&s->windowed_buffer); |
326 |
|
66 |
} |