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/lls.h"
25 #define LPC_USE_DOUBLE
30 * Apply Welch window function to audio block
32 static void apply_welch_window(const int32_t *data
, int len
, double *w_data
)
38 assert(!(len
&1)); //the optimization in r11881 does not support odd len
39 //if someone wants odd len extend the change in r11881
42 c
= 2.0 / (len
- 1.0);
49 w_data
[-i
-1] = data
[-i
-1] * w
;
50 w_data
[+i
] = data
[+i
] * w
;
55 * Calculates autocorrelation data from audio samples
56 * A Welch window function is applied before calculation.
58 void ff_lpc_compute_autocorr(const int32_t *data
, int len
, int lag
,
62 double tmp
[len
+ lag
+ 1];
63 double *data1
= tmp
+ lag
;
65 apply_welch_window(data
, len
, data1
);
71 for(j
=0; j
<lag
; j
+=2){
72 double sum0
= 1.0, sum1
= 1.0;
74 sum0
+= data1
[i
] * data1
[i
-j
];
75 sum1
+= data1
[i
] * data1
[i
-j
-1];
83 for(i
=j
-1; i
<len
; i
+=2){
84 sum
+= data1
[i
] * data1
[i
-j
]
85 + data1
[i
+1] * data1
[i
-j
+1];
92 * Quantize LPC coefficients
94 static void quantize_lpc_coefs(double *lpc_in
, int order
, int precision
,
95 int32_t *lpc_out
, int *shift
, int max_shift
, int zero_shift
)
102 /* define maximum levels */
103 qmax
= (1 << (precision
- 1)) - 1;
105 /* find maximum coefficient value */
107 for(i
=0; i
<order
; i
++) {
108 cmax
= FFMAX(cmax
, fabs(lpc_in
[i
]));
111 /* if maximum value quantizes to zero, return all zeros */
112 if(cmax
* (1 << max_shift
) < 1.0) {
114 memset(lpc_out
, 0, sizeof(int32_t) * order
);
118 /* calculate level shift which scales max coeff to available bits */
120 while((cmax
* (1 << sh
) > qmax
) && (sh
> 0)) {
124 /* since negative shift values are unsupported in decoder, scale down
125 coefficients instead */
126 if(sh
== 0 && cmax
> qmax
) {
127 double scale
= ((double)qmax
) / cmax
;
128 for(i
=0; i
<order
; i
++) {
133 /* output quantized coefficients and level shift */
135 for(i
=0; i
<order
; i
++) {
136 error
-= lpc_in
[i
] * (1 << sh
);
137 lpc_out
[i
] = av_clip(lrintf(error
), -qmax
, qmax
);
143 static int estimate_best_order(double *ref
, int min_order
, int max_order
)
148 for(i
=max_order
-1; i
>=min_order
-1; i
--) {
158 * Calculate LPC coefficients for multiple orders
160 * @param use_lpc LPC method for determining coefficients
161 * 0 = LPC with fixed pre-defined coeffs
162 * 1 = LPC with coeffs determined by Levinson-Durbin recursion
163 * 2+ = LPC with coeffs determined by Cholesky factorization using (use_lpc-1) passes.
165 int ff_lpc_calc_coefs(DSPContext
*s
,
166 const int32_t *samples
, int blocksize
, int min_order
,
167 int max_order
, int precision
,
168 int32_t coefs
[][MAX_LPC_ORDER
], int *shift
, int use_lpc
,
169 int omethod
, int max_shift
, int zero_shift
)
171 double autoc
[MAX_LPC_ORDER
+1];
172 double ref
[MAX_LPC_ORDER
];
173 double lpc
[MAX_LPC_ORDER
][MAX_LPC_ORDER
];
177 assert(max_order
>= MIN_LPC_ORDER
&& max_order
<= MAX_LPC_ORDER
&& use_lpc
> 0);
180 s
->lpc_compute_autocorr(samples
, blocksize
, max_order
, autoc
);
182 compute_lpc_coefs(autoc
, max_order
, &lpc
[0][0], MAX_LPC_ORDER
, 0, 1);
184 for(i
=0; i
<max_order
; i
++)
185 ref
[i
] = fabs(lpc
[i
][i
]);
188 double var
[MAX_LPC_ORDER
+1], av_uninit(weight
);
190 for(pass
=0; pass
<use_lpc
-1; pass
++){
191 av_init_lls(&m
[pass
&1], max_order
);
194 for(i
=max_order
; i
<blocksize
; i
++){
195 for(j
=0; j
<=max_order
; j
++)
196 var
[j
]= samples
[i
-j
];
199 double eval
, inv
, rinv
;
200 eval
= av_evaluate_lls(&m
[(pass
-1)&1], var
+1, max_order
-1);
201 eval
= (512>>pass
) + fabs(eval
- var
[0]);
204 for(j
=0; j
<=max_order
; j
++)
210 av_update_lls(&m
[pass
&1], var
, 1.0);
212 av_solve_lls(&m
[pass
&1], 0.001, 0);
215 for(i
=0; i
<max_order
; i
++){
216 for(j
=0; j
<max_order
; j
++)
217 lpc
[i
][j
]=-m
[(pass
-1)&1].coeff
[i
][j
];
218 ref
[i
]= sqrt(m
[(pass
-1)&1].variance
[i
] / weight
) * (blocksize
- max_order
) / 4000;
220 for(i
=max_order
-1; i
>0; i
--)
221 ref
[i
] = ref
[i
-1] - ref
[i
];
223 opt_order
= max_order
;
225 if(omethod
== ORDER_METHOD_EST
) {
226 opt_order
= estimate_best_order(ref
, min_order
, max_order
);
228 quantize_lpc_coefs(lpc
[i
], i
+1, precision
, coefs
[i
], &shift
[i
], max_shift
, zero_shift
);
230 for(i
=min_order
-1; i
<max_order
; i
++) {
231 quantize_lpc_coefs(lpc
[i
], i
+1, precision
, coefs
[i
], &shift
[i
], max_shift
, zero_shift
);