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.h"

25 #include "libavutil/mem_internal.h"

27 #define LPC_USE_DOUBLE

28 #include "lpc.h"

29 #include "lpc_functions.h"

30 #include "libavutil/avassert.h"

32 /**

33 * Schur recursion.

34 * Produces reflection coefficients from autocorrelation data.

35 */

36 static inline void compute_ref_coefs(const LPC_TYPE *autoc, int max_order,

37 LPC_TYPE *ref, LPC_TYPE *error)

38 {

39 LPC_TYPE err;

40 LPC_TYPE gen0[MAX_LPC_ORDER], gen1[MAX_LPC_ORDER];

42 for (int i = 0; i < max_order; i++)

43 gen0[i] = gen1[i] = autoc[i + 1];

45 err = autoc[0];

46 ref[0] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1);

47 err += gen1[0] * ref[0];

48 if (error)

49 error[0] = err;

50 for (int i = 1; i < max_order; i++) {

51 for (int j = 0; j < max_order - i; j++) {

52 gen1[j] = gen1[j + 1] + ref[i - 1] * gen0[j];

53 gen0[j] = gen1[j + 1] * ref[i - 1] + gen0[j];

54 }

55 ref[i] = -gen1[0] / ((LPC_USE_FIXED || err) ? err : 1);

56 err += gen1[0] * ref[i];

57 if (error)

58 error[i] = err;

59 }

60 }

63 /**

64 * Apply Welch window function to audio block

65 */

66 static void lpc_apply_welch_window_c(const int32_t *data, ptrdiff_t len,

67 double *w_data)

68 {

69 int i, n2;

70 double w;

71 double c;

73 if (len == 1) {

74 w_data[0] = 0.0;

75 return;

76 }

78 n2 = (len >> 1);

79 c = 2.0 / (len - 1.0);

81 if (len & 1) {

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

83 w = c - i - 1.0;

84 w = 1.0 - (w * w);

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

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

87 }

88 w_data[n2] = 0.0;

89 return;

90 }

92 w_data+=n2;

93 data+=n2;

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

95 w = c - n2 + i;

96 w = 1.0 - (w * w);

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

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

99 }

100 }

101

102 /**

103 * Calculate autocorrelation data from audio samples

104 * A Welch window function is applied before calculation.

105 */

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

107 double *autoc)

108 {

109 int i, j;

110

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

112 double sum0 = 1.0, sum1 = 1.0;

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

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

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

116 }

117 autoc[j ] = sum0;

118 autoc[j+1] = sum1;

119 }

120

121 if(j==lag){

122 double sum = 1.0;

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

124 sum += data[i] * data[i-j];

125 }

126 autoc[j] = sum;

127 }

128 }

129

130 /**

131 * Quantize LPC coefficients

132 */

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

134 int32_t *lpc_out, int *shift, int min_shift,

135 int max_shift, int zero_shift)

136 {

137 int i;

138 double cmax, error;

139 int32_t qmax;

140 int sh;

141

142 /* define maximum levels */

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

144

145 /* find maximum coefficient value */

146 cmax = 0.0;

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

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

149 }

150

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

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

153 *shift = zero_shift;

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

155 return;

156 }

157

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

159 sh = max_shift;

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

161 sh--;

162 }

163

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

165 coefficients instead */

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

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

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

169 lpc_in[i] *= scale;

170 }

171 }

172

173 /* output quantized coefficients and level shift */

174 error=0;

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

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

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

178 error -= lpc_out[i];

179 }

180 *shift = sh;

181 }

182

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

184 {

185 int i, est;

186

187 est = min_order;

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

189 if(ref[i] > 0.10) {

190 est = i+1;

191 break;

192 }

193 }

194 return est;

195 }

196

197 int ff_lpc_calc_ref_coefs(LPCContext *s,

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

199 {

200 double autoc[MAX_LPC_ORDER + 1];

201

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

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

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

205

206 return order;

207 }

208

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

210 int order, double *ref)

211 {

212 int i;

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

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

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

216

217 /* Apply windowing */

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

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

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

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

222 }

223

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

225 signal = autoc[0];

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

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

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

229 return avg_err ? signal/avg_err : NAN;

230 }

231

232 /**

233 * Calculate LPC coefficients for multiple orders

234 *

235 * @param lpc_type LPC method for determining coefficients,

236 * see #FFLPCType for details

237 */

238 int ff_lpc_calc_coefs(LPCContext *s,

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

240 int max_order, int precision,

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

242 enum FFLPCType lpc_type, int lpc_passes,

243 int omethod, int min_shift, int max_shift, int zero_shift)

244 {

245 double autoc[MAX_LPC_ORDER+1];

246 double ref[MAX_LPC_ORDER] = { 0 };

247 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];

248 int i, j, pass = 0;

249 int opt_order;

250

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

252 lpc_type > FF_LPC_TYPE_FIXED);

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

254

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

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

257 lpc_type != s->lpc_type) {

258 ff_lpc_end(s);

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

260 }

261

262 if(lpc_passes <= 0)

263 lpc_passes = 2;

264

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

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

267

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

269

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

271

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

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

274

275 pass++;

276 }

277

278 if (lpc_type == FF_LPC_TYPE_CHOLESKY) {

279 LLSModel *m = s->lls_models;

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

281 double av_uninit(weight);

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

283

284 /* Avoids initializing with an unused value when lpc_passes == 1 */

285 if (lpc_passes > 1)

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

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

288

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

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

291

292 weight=0;

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

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

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

296

297 if(pass){

298 double eval, inv, rinv;

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

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

301 inv = 1/eval;

302 rinv = sqrt(inv);

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

304 var[j] *= rinv;

305 weight += inv;

306 }else

307 weight++;

308

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

310 }

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

312 }

313

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

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

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

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

318 }

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

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

321 }

322

323 opt_order = max_order;

324

325 if(omethod == ORDER_METHOD_EST) {

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

327 i = opt_order-1;

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

329 min_shift, max_shift, zero_shift);

330 } else {

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

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

333 min_shift, max_shift, zero_shift);

334 }

335 }

336

337 return opt_order;

338 }

339

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

341 enum FFLPCType lpc_type)

342 {

343 s->blocksize = blocksize;

344 s->max_order = max_order;

345 s->lpc_type = lpc_type;

346

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

348 sizeof(*s->windowed_samples));

349 if (!s->windowed_buffer)

350 return AVERROR(ENOMEM);

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

352

353 s->lpc_apply_welch_window = lpc_apply_welch_window_c;

354 s->lpc_compute_autocorr = lpc_compute_autocorr_c;

355

356 #if ARCH_RISCV

357 ff_lpc_init_riscv(s);

358 #elif ARCH_X86

359 ff_lpc_init_x86(s);

360 #endif

361

362 return 0;

363 }

364

365 av_cold void ff_lpc_end(LPCContext *s)

366 {

367 av_freep(&s->windowed_buffer);

368 }

error

static void error(const char *err)

Definition: target_bsf_fuzzer.c:32

avpriv_solve_lls

void avpriv_solve_lls(LLSModel *m, double threshold, unsigned short min_order)

Definition: lls.c:47

LLSModel

Linear least squares model.

Definition: lls.h:37

FFLPCType

LPC analysis type.

Definition: lpc.h:42

av_clip

#define av_clip

Definition: common.h:100

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

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:238

ff_lpc_init

av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, enum FFLPCType lpc_type)

Initialize LPCContext.

Definition: lpc.c:340

lpc_apply_welch_window_c

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

Apply Welch window function to audio block.

Definition: lpc.c:66

uint8_t w

Definition: llviddspenc.c:38

FF_LPC_TYPE_CHOLESKY

@ FF_LPC_TYPE_CHOLESKY

Cholesky factorization.

Definition: lpc.h:47

#define b

Definition: input.c:42

data

const char data[16]

Definition: mxf.c:149

FFMAX

#define FFMAX(a, b)

Definition: macros.h:47

lpc.h

LPCContext

Definition: lpc.h:51

weight

const h264_weight_func weight

Definition: h264dsp_init.c:33

avassert.h

av_cold

#define av_cold

Definition: attributes.h:106

LOCAL_ALIGNED

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

Definition: mem_internal.h:124

compute_ref_coefs

static void compute_ref_coefs(const LPC_TYPE *autoc, int max_order, LPC_TYPE *ref, LPC_TYPE *error)

Schur recursion.

Definition: lpc.c:36

#define s(width, name)

Definition: cbs_vp9.c:198

av_assert0

#define av_assert0(cond)

assert() equivalent, that is always enabled.

Definition: avassert.h:41

lls.h

NAN

#define NAN

Definition: mathematics.h:115

fabs

static __device__ float fabs(float a)

Definition: cuda_runtime.h:182

NULL

#define NULL

Definition: coverity.c:32

ff_lpc_end

av_cold void ff_lpc_end(LPCContext *s)

Uninitialize LPCContext.

Definition: lpc.c:365

double

Definition: af_crystalizer.c:132

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

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

shift

static int shift(int a, int b)

Definition: bonk.c:261

MAX_LPC_ORDER

#define MAX_LPC_ORDER

Definition: lpc.h:37

MIN_LPC_ORDER

#define MIN_LPC_ORDER

Definition: lpc.h:36

ff_lpc_init_x86

void ff_lpc_init_x86(LPCContext *s)

Definition: lpc_init.c:106

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:29

M_PI

#define M_PI

Definition: mathematics.h:67

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:209

LLSModel::evaluate_lls

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