ffmpeg / libavcodec / lpc.c @ 0d8837bd
History  View  Annotate  Download (6.96 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(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 int32_t *data, int len, int lag, 
58 
double *autoc)

59 
{ 
60 
int i, j;

61 
double tmp[len + lag + 1]; 
62 
double *data1= tmp + lag;

63  
64 
apply_welch_window(data, len, data1); 
65  
66 
for(j=0; j<lag; j++) 
67 
data1[jlag]= 0.0; 
68 
data1[len] = 0.0; 
69  
70 
for(j=0; j<lag; j+=2){ 
71 
double sum0 = 1.0, sum1 = 1.0; 
72 
for(i=j; i<len; i++){

73 
sum0 += data1[i] * data1[ij]; 
74 
sum1 += data1[i] * data1[ij1];

75 
} 
76 
autoc[j ] = sum0; 
77 
autoc[j+1] = sum1;

78 
} 
79  
80 
if(j==lag){

81 
double sum = 1.0; 
82 
for(i=j1; i<len; i+=2){ 
83 
sum += data1[i ] * data1[ij ] 
84 
+ data1[i+1] * data1[ij+1]; 
85 
} 
86 
autoc[j] = sum; 
87 
} 
88 
} 
89  
90 
/**

91 
* Quantize LPC coefficients

92 
*/

93 
static void quantize_lpc_coefs(double *lpc_in, int order, int precision, 
94 
int32_t *lpc_out, int *shift, int max_shift, int zero_shift) 
95 
{ 
96 
int i;

97 
double cmax, error;

98 
int32_t qmax; 
99 
int sh;

100  
101 
/* define maximum levels */

102 
qmax = (1 << (precision  1))  1; 
103  
104 
/* find maximum coefficient value */

105 
cmax = 0.0; 
106 
for(i=0; i<order; i++) { 
107 
cmax= FFMAX(cmax, fabs(lpc_in[i])); 
108 
} 
109  
110 
/* if maximum value quantizes to zero, return all zeros */

111 
if(cmax * (1 << max_shift) < 1.0) { 
112 
*shift = zero_shift; 
113 
memset(lpc_out, 0, sizeof(int32_t) * order); 
114 
return;

115 
} 
116  
117 
/* calculate level shift which scales max coeff to available bits */

118 
sh = max_shift; 
119 
while((cmax * (1 << sh) > qmax) && (sh > 0)) { 
120 
sh; 
121 
} 
122  
123 
/* since negative shift values are unsupported in decoder, scale down

124 
coefficients instead */

125 
if(sh == 0 && cmax > qmax) { 
126 
double scale = ((double)qmax) / cmax; 
127 
for(i=0; i<order; i++) { 
128 
lpc_in[i] *= scale; 
129 
} 
130 
} 
131  
132 
/* output quantized coefficients and level shift */

133 
error=0;

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

136 
lpc_out[i] = av_clip(lrintf(error), qmax, qmax); 
137 
error = lpc_out[i]; 
138 
} 
139 
*shift = sh; 
140 
} 
141  
142 
static int estimate_best_order(double *ref, int min_order, int max_order) 
143 
{ 
144 
int i, est;

145  
146 
est = min_order; 
147 
for(i=max_order1; i>=min_order1; i) { 
148 
if(ref[i] > 0.10) { 
149 
est = i+1;

150 
break;

151 
} 
152 
} 
153 
return est;

154 
} 
155  
156 
/**

157 
* Calculate LPC coefficients for multiple orders

158 
*

159 
* @param use_lpc LPC method for determining coefficients

160 
* 0 = LPC with fixed predefined coeffs

161 
* 1 = LPC with coeffs determined by LevinsonDurbin recursion

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

163 
*/

164 
int ff_lpc_calc_coefs(LPCContext *s,

165 
const int32_t *samples, int blocksize, int min_order, 
166 
int max_order, int precision, 
167 
int32_t coefs[][MAX_LPC_ORDER], int *shift,

168 
enum AVLPCType lpc_type, int lpc_passes, 
169 
int omethod, int max_shift, int zero_shift) 
170 
{ 
171 
double autoc[MAX_LPC_ORDER+1]; 
172 
double ref[MAX_LPC_ORDER];

173 
double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];

174 
int i, j, pass;

175 
int opt_order;

176  
177 
assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && 
178 
lpc_type > AV_LPC_TYPE_FIXED); 
179  
180 
if (lpc_type == AV_LPC_TYPE_LEVINSON) {

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

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

193  
194 
weight=0;

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

196 
for(j=0; j<=max_order; j++) 
197 
var[j]= samples[ij]; 
198  
199 
if(pass){

200 
double eval, inv, rinv;

201 
eval= av_evaluate_lls(&m[(pass1)&1], var+1, max_order1); 
202 
eval= (512>>pass) + fabs(eval  var[0]); 
203 
inv = 1/eval;

204 
rinv = sqrt(inv); 
205 
for(j=0; j<=max_order; j++) 
206 
var[j] *= rinv; 
207 
weight += inv; 
208 
}else

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

223 
} 
224 
opt_order = max_order; 
225  
226 
if(omethod == ORDER_METHOD_EST) {

227 
opt_order = estimate_best_order(ref, min_order, max_order); 
228 
i = opt_order1;

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

230 
} else {

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

233 
} 
234 
} 
235  
236 
return opt_order;

237 
} 
238  
239 
av_cold void ff_lpc_init(LPCContext *s)

240 
{ 
241 
s>lpc_compute_autocorr = lpc_compute_autocorr_c; 
242  
243 
if (HAVE_MMX)

244 
ff_lpc_init_x86(s); 
245 
} 