3 * Copyright (c) 2006 Justin Ruggles <justin.ruggles@gmail.com>
5 * This file is part of FFmpeg.
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.
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.
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
22 #include "libavutil/common.h"
23 #include "libavutil/lls.h"
25 #define LPC_USE_DOUBLE
27 #include "libavutil/avassert.h"
31 * Apply Welch window function to audio block
33 static void lpc_apply_welch_window_c(const int32_t *data
, int len
,
40 /* The optimization in commit fa4ed8c does not support odd len.
41 * If someone wants odd len extend that change. */
42 av_assert2(!(len
& 1));
45 c
= 2.0 / (len
- 1.0);
52 w_data
[-i
-1] = data
[-i
-1] * w
;
53 w_data
[+i
] = data
[+i
] * w
;
58 * Calculate autocorrelation data from audio samples
59 * A Welch window function is applied before calculation.
61 static void lpc_compute_autocorr_c(const double *data
, int len
, int lag
,
66 for(j
=0; j
<lag
; j
+=2){
67 double sum0
= 1.0, sum1
= 1.0;
69 sum0
+= data
[i
] * data
[i
-j
];
70 sum1
+= data
[i
] * data
[i
-j
-1];
78 for(i
=j
-1; i
<len
; i
+=2){
79 sum
+= data
[i
] * data
[i
-j
]
80 + data
[i
+1] * data
[i
-j
+1];
87 * Quantize LPC coefficients
89 static void quantize_lpc_coefs(double *lpc_in
, int order
, int precision
,
90 int32_t *lpc_out
, int *shift
, int max_shift
, int zero_shift
)
97 /* define maximum levels */
98 qmax
= (1 << (precision
- 1)) - 1;
100 /* find maximum coefficient value */
102 for(i
=0; i
<order
; i
++) {
103 cmax
= FFMAX(cmax
, fabs(lpc_in
[i
]));
106 /* if maximum value quantizes to zero, return all zeros */
107 if(cmax
* (1 << max_shift
) < 1.0) {
109 memset(lpc_out
, 0, sizeof(int32_t) * order
);
113 /* calculate level shift which scales max coeff to available bits */
115 while((cmax
* (1 << sh
) > qmax
) && (sh
> 0)) {
119 /* since negative shift values are unsupported in decoder, scale down
120 coefficients instead */
121 if(sh
== 0 && cmax
> qmax
) {
122 double scale
= ((double)qmax
) / cmax
;
123 for(i
=0; i
<order
; i
++) {
128 /* output quantized coefficients and level shift */
130 for(i
=0; i
<order
; i
++) {
131 error
-= lpc_in
[i
] * (1 << sh
);
132 lpc_out
[i
] = av_clip(lrintf(error
), -qmax
, qmax
);
138 static int estimate_best_order(double *ref
, int min_order
, int max_order
)
143 for(i
=max_order
-1; i
>=min_order
-1; i
--) {
152 int ff_lpc_calc_ref_coefs(LPCContext
*s
,
153 const int32_t *samples
, int order
, double *ref
)
155 double autoc
[MAX_LPC_ORDER
+ 1];
157 s
->lpc_apply_welch_window(samples
, s
->blocksize
, s
->windowed_samples
);
158 s
->lpc_compute_autocorr(s
->windowed_samples
, s
->blocksize
, order
, autoc
);
159 compute_ref_coefs(autoc
, order
, ref
, NULL
);
165 * Calculate LPC coefficients for multiple orders
167 * @param lpc_type LPC method for determining coefficients,
168 * see #FFLPCType for details
170 int ff_lpc_calc_coefs(LPCContext
*s
,
171 const int32_t *samples
, int blocksize
, int min_order
,
172 int max_order
, int precision
,
173 int32_t coefs
[][MAX_LPC_ORDER
], int *shift
,
174 enum FFLPCType lpc_type
, int lpc_passes
,
175 int omethod
, int max_shift
, int zero_shift
)
177 double autoc
[MAX_LPC_ORDER
+1];
178 double ref
[MAX_LPC_ORDER
] = { 0 };
179 double lpc
[MAX_LPC_ORDER
][MAX_LPC_ORDER
];
183 av_assert2(max_order
>= MIN_LPC_ORDER
&& max_order
<= MAX_LPC_ORDER
&&
184 lpc_type
> FF_LPC_TYPE_FIXED
);
185 av_assert0(lpc_type
== FF_LPC_TYPE_CHOLESKY
|| lpc_type
== FF_LPC_TYPE_LEVINSON
);
187 /* reinit LPC context if parameters have changed */
188 if (blocksize
!= s
->blocksize
|| max_order
!= s
->max_order
||
189 lpc_type
!= s
->lpc_type
) {
191 ff_lpc_init(s
, blocksize
, max_order
, lpc_type
);
197 if (lpc_type
== FF_LPC_TYPE_LEVINSON
|| (lpc_type
== FF_LPC_TYPE_CHOLESKY
&& lpc_passes
> 1)) {
198 s
->lpc_apply_welch_window(samples
, blocksize
, s
->windowed_samples
);
200 s
->lpc_compute_autocorr(s
->windowed_samples
, blocksize
, max_order
, autoc
);
202 compute_lpc_coefs(autoc
, max_order
, &lpc
[0][0], MAX_LPC_ORDER
, 0, 1);
204 for(i
=0; i
<max_order
; i
++)
205 ref
[i
] = fabs(lpc
[i
][i
]);
210 if (lpc_type
== FF_LPC_TYPE_CHOLESKY
) {
211 LLSModel
*m
= s
->lls_models
;
212 LOCAL_ALIGNED(32, double, var
, [FFALIGN(MAX_LPC_ORDER
+1,4)]);
213 double av_uninit(weight
);
214 memset(var
, 0, FFALIGN(MAX_LPC_ORDER
+1,4)*sizeof(*var
));
216 for(j
=0; j
<max_order
; j
++)
217 m
[0].coeff
[max_order
-1][j
] = -lpc
[max_order
-1][j
];
219 for(; pass
<lpc_passes
; pass
++){
220 avpriv_init_lls(&m
[pass
&1], max_order
);
223 for(i
=max_order
; i
<blocksize
; i
++){
224 for(j
=0; j
<=max_order
; j
++)
225 var
[j
]= samples
[i
-j
];
228 double eval
, inv
, rinv
;
229 eval
= m
[pass
&1].evaluate_lls(&m
[(pass
-1)&1], var
+1, max_order
-1);
230 eval
= (512>>pass
) + fabs(eval
- var
[0]);
233 for(j
=0; j
<=max_order
; j
++)
239 m
[pass
&1].update_lls(&m
[pass
&1], var
);
241 avpriv_solve_lls(&m
[pass
&1], 0.001, 0);
244 for(i
=0; i
<max_order
; i
++){
245 for(j
=0; j
<max_order
; j
++)
246 lpc
[i
][j
]=-m
[(pass
-1)&1].coeff
[i
][j
];
247 ref
[i
]= sqrt(m
[(pass
-1)&1].variance
[i
] / weight
) * (blocksize
- max_order
) / 4000;
249 for(i
=max_order
-1; i
>0; i
--)
250 ref
[i
] = ref
[i
-1] - ref
[i
];
253 opt_order
= max_order
;
255 if(omethod
== ORDER_METHOD_EST
) {
256 opt_order
= estimate_best_order(ref
, min_order
, max_order
);
258 quantize_lpc_coefs(lpc
[i
], i
+1, precision
, coefs
[i
], &shift
[i
], max_shift
, zero_shift
);
260 for(i
=min_order
-1; i
<max_order
; i
++) {
261 quantize_lpc_coefs(lpc
[i
], i
+1, precision
, coefs
[i
], &shift
[i
], max_shift
, zero_shift
);
268 av_cold
int ff_lpc_init(LPCContext
*s
, int blocksize
, int max_order
,
269 enum FFLPCType lpc_type
)
271 s
->blocksize
= blocksize
;
272 s
->max_order
= max_order
;
273 s
->lpc_type
= lpc_type
;
275 s
->windowed_buffer
= av_mallocz((blocksize
+ 2 + FFALIGN(max_order
, 4)) *
276 sizeof(*s
->windowed_samples
));
277 if (!s
->windowed_buffer
)
278 return AVERROR(ENOMEM
);
279 s
->windowed_samples
= s
->windowed_buffer
+ FFALIGN(max_order
, 4);
281 s
->lpc_apply_welch_window
= lpc_apply_welch_window_c
;
282 s
->lpc_compute_autocorr
= lpc_compute_autocorr_c
;
290 av_cold
void ff_lpc_end(LPCContext
*s
)
292 av_freep(&s
->windowed_buffer
);