ffmpeg / libavcodec / lpc.c @ 2912e87a
History  View  Annotate  Download (7.83 KB)
1 
/**


2 
* LPC utility code

3 
* Copyright (c) 2006 Justin Ruggles <justin.ruggles@gmail.com>

4 
*

5 
* This file is part of Libav.

6 
*

7 
* Libav 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 
* Libav 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 Libav; if not, write to the Free Software

19 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

20 
*/

21  
22 
#include "libavutil/lls.h" 
23  
24 
#define LPC_USE_DOUBLE

25 
#include "lpc.h" 
26  
27  
28 
/**

29 
* Apply Welch window function to audio block

30 
*/

31 
static void lpc_apply_welch_window_c(const int32_t *data, int len, 
32 
double *w_data)

33 
{ 
34 
int i, n2;

35 
double w;

36 
double c;

37  
38 
assert(!(len&1)); //the optimization in r11881 does not support odd len 
39 
//if someone wants odd len extend the change in r11881

40  
41 
n2 = (len >> 1);

42 
c = 2.0 / (len  1.0); 
43  
44 
w_data+=n2; 
45 
data+=n2; 
46 
for(i=0; i<n2; i++) { 
47 
w = c  n2 + i; 
48 
w = 1.0  (w * w); 
49 
w_data[i1] = data[i1] * w; 
50 
w_data[+i ] = data[+i ] * w; 
51 
} 
52 
} 
53  
54 
/**

55 
* Calculate autocorrelation data from audio samples

56 
* A Welch window function is applied before calculation.

57 
*/

58 
static void lpc_compute_autocorr_c(const double *data, int len, int lag, 
59 
double *autoc)

60 
{ 
61 
int i, j;

62  
63 
for(j=0; j<lag; j+=2){ 
64 
double sum0 = 1.0, sum1 = 1.0; 
65 
for(i=j; i<len; i++){

66 
sum0 += data[i] * data[ij]; 
67 
sum1 += data[i] * data[ij1];

68 
} 
69 
autoc[j ] = sum0; 
70 
autoc[j+1] = sum1;

71 
} 
72  
73 
if(j==lag){

74 
double sum = 1.0; 
75 
for(i=j1; i<len; i+=2){ 
76 
sum += data[i ] * data[ij ] 
77 
+ data[i+1] * data[ij+1]; 
78 
} 
79 
autoc[j] = sum; 
80 
} 
81 
} 
82  
83 
/**

84 
* Quantize LPC coefficients

85 
*/

86 
static void quantize_lpc_coefs(double *lpc_in, int order, int precision, 
87 
int32_t *lpc_out, int *shift, int max_shift, int zero_shift) 
88 
{ 
89 
int i;

90 
double cmax, error;

91 
int32_t qmax; 
92 
int sh;

93  
94 
/* define maximum levels */

95 
qmax = (1 << (precision  1))  1; 
96  
97 
/* find maximum coefficient value */

98 
cmax = 0.0; 
99 
for(i=0; i<order; i++) { 
100 
cmax= FFMAX(cmax, fabs(lpc_in[i])); 
101 
} 
102  
103 
/* if maximum value quantizes to zero, return all zeros */

104 
if(cmax * (1 << max_shift) < 1.0) { 
105 
*shift = zero_shift; 
106 
memset(lpc_out, 0, sizeof(int32_t) * order); 
107 
return;

108 
} 
109  
110 
/* calculate level shift which scales max coeff to available bits */

111 
sh = max_shift; 
112 
while((cmax * (1 << sh) > qmax) && (sh > 0)) { 
113 
sh; 
114 
} 
115  
116 
/* since negative shift values are unsupported in decoder, scale down

117 
coefficients instead */

118 
if(sh == 0 && cmax > qmax) { 
119 
double scale = ((double)qmax) / cmax; 
120 
for(i=0; i<order; i++) { 
121 
lpc_in[i] *= scale; 
122 
} 
123 
} 
124  
125 
/* output quantized coefficients and level shift */

126 
error=0;

127 
for(i=0; i<order; i++) { 
128 
error = lpc_in[i] * (1 << sh);

129 
lpc_out[i] = av_clip(lrintf(error), qmax, qmax); 
130 
error = lpc_out[i]; 
131 
} 
132 
*shift = sh; 
133 
} 
134  
135 
static int estimate_best_order(double *ref, int min_order, int max_order) 
136 
{ 
137 
int i, est;

138  
139 
est = min_order; 
140 
for(i=max_order1; i>=min_order1; i) { 
141 
if(ref[i] > 0.10) { 
142 
est = i+1;

143 
break;

144 
} 
145 
} 
146 
return est;

147 
} 
148  
149 
/**

150 
* Calculate LPC coefficients for multiple orders

151 
*

152 
* @param use_lpc LPC method for determining coefficients

153 
* 0 = LPC with fixed predefined coeffs

154 
* 1 = LPC with coeffs determined by LevinsonDurbin recursion

155 
* 2+ = LPC with coeffs determined by Cholesky factorization using (use_lpc1) passes.

156 
*/

157 
int ff_lpc_calc_coefs(LPCContext *s,

158 
const int32_t *samples, int blocksize, int min_order, 
159 
int max_order, int precision, 
160 
int32_t coefs[][MAX_LPC_ORDER], int *shift,

161 
enum AVLPCType lpc_type, int lpc_passes, 
162 
int omethod, int max_shift, int zero_shift) 
163 
{ 
164 
double autoc[MAX_LPC_ORDER+1]; 
165 
double ref[MAX_LPC_ORDER];

166 
double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];

167 
int i, j, pass;

168 
int opt_order;

169  
170 
assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && 
171 
lpc_type > AV_LPC_TYPE_FIXED); 
172  
173 
/* reinit LPC context if parameters have changed */

174 
if (blocksize != s>blocksize  max_order != s>max_order 

175 
lpc_type != s>lpc_type) { 
176 
ff_lpc_end(s); 
177 
ff_lpc_init(s, blocksize, max_order, lpc_type); 
178 
} 
179  
180 
if (lpc_type == AV_LPC_TYPE_LEVINSON) {

181 
double *windowed_samples = s>windowed_samples + max_order;

182  
183 
s>lpc_apply_welch_window(samples, blocksize, windowed_samples); 
184  
185 
s>lpc_compute_autocorr(windowed_samples, blocksize, max_order, autoc); 
186  
187 
compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); 
188  
189 
for(i=0; i<max_order; i++) 
190 
ref[i] = fabs(lpc[i][i]); 
191 
} else if (lpc_type == AV_LPC_TYPE_CHOLESKY) { 
192 
LLSModel m[2];

193 
double var[MAX_LPC_ORDER+1], av_uninit(weight); 
194  
195 
for(pass=0; pass<lpc_passes; pass++){ 
196 
av_init_lls(&m[pass&1], max_order);

197  
198 
weight=0;

199 
for(i=max_order; i<blocksize; i++){

200 
for(j=0; j<=max_order; j++) 
201 
var[j]= samples[ij]; 
202  
203 
if(pass){

204 
double eval, inv, rinv;

205 
eval= av_evaluate_lls(&m[(pass1)&1], var+1, max_order1); 
206 
eval= (512>>pass) + fabs(eval  var[0]); 
207 
inv = 1/eval;

208 
rinv = sqrt(inv); 
209 
for(j=0; j<=max_order; j++) 
210 
var[j] *= rinv; 
211 
weight += inv; 
212 
}else

213 
weight++; 
214  
215 
av_update_lls(&m[pass&1], var, 1.0); 
216 
} 
217 
av_solve_lls(&m[pass&1], 0.001, 0); 
218 
} 
219  
220 
for(i=0; i<max_order; i++){ 
221 
for(j=0; j<max_order; j++) 
222 
lpc[i][j]=m[(pass1)&1].coeff[i][j]; 
223 
ref[i]= sqrt(m[(pass1)&1].variance[i] / weight) * (blocksize  max_order) / 4000; 
224 
} 
225 
for(i=max_order1; i>0; i) 
226 
ref[i] = ref[i1]  ref[i];

227 
} 
228 
opt_order = max_order; 
229  
230 
if(omethod == ORDER_METHOD_EST) {

231 
opt_order = estimate_best_order(ref, min_order, max_order); 
232 
i = opt_order1;

233 
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);

234 
} else {

235 
for(i=min_order1; i<max_order; i++) { 
236 
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);

237 
} 
238 
} 
239  
240 
return opt_order;

241 
} 
242  
243 
av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, 
244 
enum AVLPCType lpc_type)

245 
{ 
246 
s>blocksize = blocksize; 
247 
s>max_order = max_order; 
248 
s>lpc_type = lpc_type; 
249  
250 
if (lpc_type == AV_LPC_TYPE_LEVINSON) {

251 
s>windowed_samples = av_mallocz((blocksize + max_order + 2) *

252 
sizeof(*s>windowed_samples));

253 
if (!s>windowed_samples)

254 
return AVERROR(ENOMEM);

255 
} else {

256 
s>windowed_samples = NULL;

257 
} 
258  
259 
s>lpc_apply_welch_window = lpc_apply_welch_window_c; 
260 
s>lpc_compute_autocorr = lpc_compute_autocorr_c; 
261  
262 
if (HAVE_MMX)

263 
ff_lpc_init_x86(s); 
264  
265 
return 0; 
266 
} 
267  
268 
av_cold void ff_lpc_end(LPCContext *s)

269 
{ 
270 
av_freep(&s>windowed_samples); 
271 
} 