ffmpeg / libavcodec / lpc.c @ 7101b185
History  View  Annotate  Download (7.78 KB)
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 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 apply_welch_window_c(const int32_t *data, int len, double *w_data) 
32 
{ 
33 
int i, n2;

34 
double w;

35 
double c;

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

39  
40 
n2 = (len >> 1);

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

54 
* Calculate autocorrelation data from audio samples

55 
* A Welch window function is applied before calculation.

56 
*/

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

59 
{ 
60 
int i, j;

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

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

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

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

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

83 
* Quantize LPC coefficients

84 
*/

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

89 
double cmax, error;

90 
int32_t qmax; 
91 
int sh;

92  
93 
/* define maximum levels */

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

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

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

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

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

116 
coefficients instead */

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

125 
error=0;

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

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

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

142 
break;

143 
} 
144 
} 
145 
return est;

146 
} 
147  
148 
/**

149 
* Calculate LPC coefficients for multiple orders

150 
*

151 
* @param use_lpc LPC method for determining coefficients

152 
* 0 = LPC with fixed predefined coeffs

153 
* 1 = LPC with coeffs determined by LevinsonDurbin recursion

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

155 
*/

156 
int ff_lpc_calc_coefs(LPCContext *s,

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

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

165 
double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];

166 
int i, j, pass;

167 
int opt_order;

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

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

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

180 
double *windowed_samples = s>windowed_samples + max_order;

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

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

196  
197 
weight=0;

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

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

203 
double eval, inv, rinv;

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

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

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

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

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

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

233 
} else {

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

236 
} 
237 
} 
238  
239 
return opt_order;

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

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

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

251 
sizeof(*s>windowed_samples));

252 
if (!s>windowed_samples)

253 
return AVERROR(ENOMEM);

254 
} else {

255 
s>windowed_samples = NULL;

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

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

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