FFmpeg: libavcodec/lpc.c Source File

FFmpeg

[フレーム]

libavcodec

lpc.c

Go to the documentation of this file.

1 /*

2 * LPC utility code

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 02110-1301 USA

20 */

22 #include "libavutil/common.h"

23 #include "libavutil/lls.h"

24 #include "libavutil/mem_internal.h"

26 #define LPC_USE_DOUBLE

27 #include "lpc.h"

28 #include "libavutil/avassert.h"

31 /**

32 * Apply Welch window function to audio block

33 */

34 static void lpc_apply_welch_window_c(const int32_t *data, int len,

35 double *w_data)

36 {

37 int i, n2;

38 double w;

39 double c;

41 n2 = (len >> 1);

42 c = 2.0 / (len - 1.0);

44 if (len & 1) {

45 for(i=0; i<n2; i++) {

46 w = c - i - 1.0;

47 w = 1.0 - (w * w);

48 w_data[i] = data[i] * w;

49 w_data[len-1-i] = data[len-1-i] * w;

50 }

51 return;

52 }

54 w_data+=n2;

55 data+=n2;

56 for(i=0; i<n2; i++) {

57 w = c - n2 + i;

58 w = 1.0 - (w * w);

59 w_data[-i-1] = data[-i-1] * w;

60 w_data[+i ] = data[+i ] * w;

61 }

62 }

64 /**

65 * Calculate autocorrelation data from audio samples

66 * A Welch window function is applied before calculation.

67 */

68 static void lpc_compute_autocorr_c(const double *data, int len, int lag,

69 double *autoc)

70 {

71 int i, j;

73 for(j=0; j<lag; j+=2){

74 double sum0 = 1.0, sum1 = 1.0;

75 for(i=j; i<len; i++){

76 sum0 += data[i] * data[i-j];

77 sum1 += data[i] * data[i-j-1];

78 }

79 autoc[j ] = sum0;

80 autoc[j+1] = sum1;

81 }

83 if(j==lag){

84 double sum = 1.0;

85 for(i=j-1; i<len; i+=2){

86 sum += data[i ] * data[i-j ]

87 + data[i+1] * data[i-j+1];

88 }

89 autoc[j] = sum;

90 }

91 }

93 /**

94 * Quantize LPC coefficients

95 */

96 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,

97 int32_t *lpc_out, int *shift, int min_shift,

98 int max_shift, int zero_shift)

99 {

100 int i;

101 double cmax, error;

102 int32_t qmax;

103 int sh;

104

105 /* define maximum levels */

106 qmax = (1 << (precision - 1)) - 1;

107

108 /* find maximum coefficient value */

109 cmax = 0.0;

110 for(i=0; i<order; i++) {

111 cmax= FFMAX(cmax, fabs(lpc_in[i]));

112 }

113

114 /* if maximum value quantizes to zero, return all zeros */

115 if(cmax * (1 << max_shift) < 1.0) {

116 *shift = zero_shift;

117 memset(lpc_out, 0, sizeof(int32_t) * order);

118 return;

119 }

120

121 /* calculate level shift which scales max coeff to available bits */

122 sh = max_shift;

123 while((cmax * (1 << sh) > qmax) && (sh > min_shift)) {

124 sh--;

125 }

126

127 /* since negative shift values are unsupported in decoder, scale down

128 coefficients instead */

129 if(sh == 0 && cmax > qmax) {

130 double scale = ((double)qmax) / cmax;

131 for(i=0; i<order; i++) {

132 lpc_in[i] *= scale;

133 }

134 }

135

136 /* output quantized coefficients and level shift */

137 error=0;

138 for(i=0; i<order; i++) {

139 error -= lpc_in[i] * (1 << sh);

140 lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);

141 error -= lpc_out[i];

142 }

143 *shift = sh;

144 }

145

146 static int estimate_best_order(double *ref, int min_order, int max_order)

147 {

148 int i, est;

149

150 est = min_order;

151 for(i=max_order-1; i>=min_order-1; i--) {

152 if(ref[i] > 0.10) {

153 est = i+1;

154 break;

155 }

156 }

157 return est;

158 }

159

160 int ff_lpc_calc_ref_coefs(LPCContext *s,

161 const int32_t *samples, int order, double *ref)

162 {

163 double autoc[MAX_LPC_ORDER + 1];

164

165 s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples);

166 s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc);

167 compute_ref_coefs(autoc, order, ref, NULL);

168

169 return order;

170 }

171

172 double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len,

173 int order, double *ref)

174 {

175 int i;

176 double signal = 0.0f, avg_err = 0.0f;

177 double autoc[MAX_LPC_ORDER+1] = {0}, error[MAX_LPC_ORDER+1] = {0};

178 const double a = 0.5f, b = 1.0f - a;

179

180 /* Apply windowing */

181 for (i = 0; i <= len / 2; i++) {

182 double weight = a - b*cos((2*M_PI*i)/(len - 1));

183 s->windowed_samples[i] = weight*samples[i];

184 s->windowed_samples[len-1-i] = weight*samples[len-1-i];

185 }

186

187 s->lpc_compute_autocorr(s->windowed_samples, len, order, autoc);

188 signal = autoc[0];

189 compute_ref_coefs(autoc, order, ref, error);

190 for (i = 0; i < order; i++)

191 avg_err = (avg_err + error[i])/2.0f;

192 return avg_err ? signal/avg_err : NAN;

193 }

194

195 /**

196 * Calculate LPC coefficients for multiple orders

197 *

198 * @param lpc_type LPC method for determining coefficients,

199 * see #FFLPCType for details

200 */

201 int ff_lpc_calc_coefs(LPCContext *s,

202 const int32_t *samples, int blocksize, int min_order,

203 int max_order, int precision,

204 int32_t coefs[][MAX_LPC_ORDER], int *shift,

205 enum FFLPCType lpc_type, int lpc_passes,

206 int omethod, int min_shift, int max_shift, int zero_shift)

207 {

208 double autoc[MAX_LPC_ORDER+1];

209 double ref[MAX_LPC_ORDER] = { 0 };

210 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];

211 int i, j, pass = 0;

212 int opt_order;

213

214 av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER &&

215 lpc_type > FF_LPC_TYPE_FIXED);

216 av_assert0(lpc_type == FF_LPC_TYPE_CHOLESKY || lpc_type == FF_LPC_TYPE_LEVINSON);

217

218 /* reinit LPC context if parameters have changed */

219 if (blocksize != s->blocksize || max_order != s->max_order ||

220 lpc_type != s->lpc_type) {

221 ff_lpc_end(s);

222 ff_lpc_init(s, blocksize, max_order, lpc_type);

223 }

224

225 if(lpc_passes <= 0)

226 lpc_passes = 2;

227

228 if (lpc_type == FF_LPC_TYPE_LEVINSON || (lpc_type == FF_LPC_TYPE_CHOLESKY && lpc_passes > 1)) {

229 s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples);

230

231 s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc);

232

233 compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);

234

235 for(i=0; i<max_order; i++)

236 ref[i] = fabs(lpc[i][i]);

237

238 pass++;

239 }

240

241 if (lpc_type == FF_LPC_TYPE_CHOLESKY) {

242 LLSModel *m = s->lls_models;

243 LOCAL_ALIGNED(32, double, var, [FFALIGN(MAX_LPC_ORDER+1,4)]);

244 double av_uninit(weight);

245 memset(var, 0, FFALIGN(MAX_LPC_ORDER+1,4)*sizeof(*var));

246

247 for(j=0; j<max_order; j++)

248 m[0].coeff[max_order-1][j] = -lpc[max_order-1][j];

249

250 for(; pass<lpc_passes; pass++){

251 avpriv_init_lls(&m[pass&1], max_order);

252

253 weight=0;

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

255 for(j=0; j<=max_order; j++)

256 var[j]= samples[i-j];

257

258 if(pass){

259 double eval, inv, rinv;

260 eval= m[pass&1].evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);

261 eval= (512>>pass) + fabs(eval - var[0]);

262 inv = 1/eval;

263 rinv = sqrt(inv);

264 for(j=0; j<=max_order; j++)

265 var[j] *= rinv;

266 weight += inv;

267 }else

268 weight++;

269

270 m[pass&1].update_lls(&m[pass&1], var);

271 }

272 avpriv_solve_lls(&m[pass&1], 0.001, 0);

273 }

274

275 for(i=0; i<max_order; i++){

276 for(j=0; j<max_order; j++)

277 lpc[i][j]=-m[(pass-1)&1].coeff[i][j];

278 ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;

279 }

280 for(i=max_order-1; i>0; i--)

281 ref[i] = ref[i-1] - ref[i];

282 }

283

284 opt_order = max_order;

285

286 if(omethod == ORDER_METHOD_EST) {

287 opt_order = estimate_best_order(ref, min_order, max_order);

288 i = opt_order-1;

289 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],

290 min_shift, max_shift, zero_shift);

291 } else {

292 for(i=min_order-1; i<max_order; i++) {

293 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i],

294 min_shift, max_shift, zero_shift);

295 }

296 }

297

298 return opt_order;

299 }

300

301 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order,

302 enum FFLPCType lpc_type)

303 {

304 s->blocksize = blocksize;

305 s->max_order = max_order;

306 s->lpc_type = lpc_type;

307

308 s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) *

309 sizeof(*s->windowed_samples));

310 if (!s->windowed_buffer)

311 return AVERROR(ENOMEM);

312 s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4);

313

314 s->lpc_apply_welch_window = lpc_apply_welch_window_c;

315 s->lpc_compute_autocorr = lpc_compute_autocorr_c;

316

317 #if ARCH_X86

318 ff_lpc_init_x86(s);

319 #endif

320

321 return 0;

322 }

323

324 av_cold void ff_lpc_end(LPCContext *s)

325 {

326 av_freep(&s->windowed_buffer);

327 }

error

static void error(const char *err)

Definition: target_bsf_fuzzer.c:31

LLSModel

Linear least squares model.

Definition: lls.h:37

FFLPCType

LPC analysis type.

Definition: lpc.h:43

av_clip

#define av_clip

Definition: common.h:95

AVERROR

Filter the word "frame" indicates either a video frame or a group of audio as stored in an AVFrame structure Format for each input and each output the list of supported formats For video that means pixel format For audio that means channel sample they are references to shared objects When the negotiation mechanism computes the intersection of the formats supported at each end of a all references to both lists are replaced with a reference to the intersection And when a single format is eventually chosen for a link amongst the remaining all references to the list are updated That means that if a filter requires that its input and output have the same format amongst a supported all it has to do is use a reference to the same list of formats query_formats can leave some formats unset and return AVERROR(EAGAIN) to cause the negotiation mechanism toagain later. That can be used by filters with complex requirements to use the format negotiated on one link to set the formats supported on another. Frame references ownership and permissions

mem_internal.h

uint8_t w

Definition: llviddspenc.c:38

FF_LPC_TYPE_CHOLESKY

@ FF_LPC_TYPE_CHOLESKY

Cholesky factorization.

Definition: lpc.h:48

compute_lpc_coefs

static int AAC_RENAME() compute_lpc_coefs(const LPC_TYPE *autoc, int max_order, LPC_TYPE *lpc, int lpc_stride, int fail, int normalize)

Levinson-Durbin recursion.

Definition: lpc.h:166

#define b

Definition: input.c:34

data

const char data[16]

Definition: mxf.c:143

estimate_best_order

static int estimate_best_order(double *ref, int min_order, int max_order)

Definition: lpc.c:146

FFMAX

#define FFMAX(a, b)

Definition: macros.h:47

lpc.h

lpc_apply_welch_window_c

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

Apply Welch window function to audio block.

Definition: lpc.c:34

LPCContext

Definition: lpc.h:52

scale

static av_always_inline float scale(float x, float s)

Definition: vf_v360.c:1389

avassert.h

av_cold

#define av_cold

Definition: attributes.h:90

LOCAL_ALIGNED

#define LOCAL_ALIGNED(a, t, v,...)

Definition: mem_internal.h:113

#define s(width, name)

Definition: cbs_vp9.c:256

av_assert0

#define av_assert0(cond)

assert() equivalent, that is always enabled.

Definition: avassert.h:37

lls.h

NAN

#define NAN

Definition: mathematics.h:64

pass

#define pass

Definition: fft_template.c:608

fabs

static __device__ float fabs(float a)

Definition: cuda_runtime.h:182

NULL

#define NULL

Definition: coverity.c:32

double

Definition: af_crystalizer.c:132

avpriv_init_lls

av_cold void avpriv_init_lls(LLSModel *m, int indep_count)

Definition: lls.c:114

ff_lpc_calc_coefs

int ff_lpc_calc_coefs(LPCContext *s, const int32_t *samples, int blocksize, int min_order, int max_order, int precision, int32_t coefs[][MAX_LPC_ORDER], int *shift, enum FFLPCType lpc_type, int lpc_passes, int omethod, int min_shift, int max_shift, int zero_shift)

Calculate LPC coefficients for multiple orders.

Definition: lpc.c:201

Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c

Definition: undefined.txt:32

weight

static int weight(int i, int blen, int offset)

Definition: diracdec.c:1562

ff_lpc_calc_ref_coefs_f

double ff_lpc_calc_ref_coefs_f(LPCContext *s, const float *samples, int len, int order, double *ref)

Definition: lpc.c:172

LLSModel::update_lls

void(* update_lls)(struct LLSModel *m, const double *var)

Take the outer-product of var[] with itself, and add to the covariance matrix.

Definition: lls.h:49

MAX_LPC_ORDER

#define MAX_LPC_ORDER

Definition: lpc.h:38

MIN_LPC_ORDER

#define MIN_LPC_ORDER

Definition: lpc.h:37

The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a

Definition: undefined.txt:41

ORDER_METHOD_EST

#define ORDER_METHOD_EST

Definition: lpc.h:30

M_PI

#define M_PI

Definition: mathematics.h:52

ff_lpc_calc_ref_coefs

int ff_lpc_calc_ref_coefs(LPCContext *s, const int32_t *samples, int order, double *ref)

Definition: lpc.c:160

LLSModel::evaluate_lls

double(* evaluate_lls)(struct LLSModel *m, const double *var, int order)