Skip to content

Commit

Permalink
reduce calling levinson--durbin algorithm.
Browse files Browse the repository at this point in the history
  • Loading branch information
@aikiriao
aikiriao committed Jun 22, 2024
1 parent e6f1888 commit 7394a2c
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 113 deletions.
10 changes: 2 additions & 8 deletions libs/lpc/include/lpc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,16 +49,10 @@ LPCApiResult LPCCalculator_CalculateLPCCoefficients(
LPCWindowType window_type, double regular_term);

/* Levinson-Durbin再帰計算により与えられた次数まで全てのLPC係数を求める（倍精度） */
/* error_varsは0次の誤差分散（分散）からmax_coef_order次の分散まで求めるためerror_varsのサイズはmax_coef_order+1要する */
LPCApiResult LPCCalculator_CalculateMultipleLPCCoefficients(
struct LPCCalculator* lpcc,
const double* data, uint32_t num_samples, double **lpc_coefs, uint32_t max_coef_order,
LPCWindowType window_type, double regular_term);

/* Levinson-Durbin再帰計算により与えられた次数までの残差分散を求める（倍精度） */
/* 0次の誤差分散（分散）からmax_coef_order次の分散まで求めるためerror_varsのサイズはmax_coef_order+1要する */
LPCApiResult LPCCalculator_CalculateErrorVariances(
struct LPCCalculator* lpcc,
const double* data, uint32_t num_samples, double *error_vars, uint32_t max_coef_order,
const double* data, uint32_t num_samples, double **lpc_coefs, double *error_vars, uint32_t max_coef_order,
LPCWindowType window_type, double regular_term);

/* 補助関数法よりLPC係数を求める（倍精度） */
Expand Down
80 changes: 32 additions & 48 deletions libs/lpc/src/lpc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,24 @@ static LPCError LPC_ApplyWindow(
return LPC_ERROR_OK;
}

/* 窓の二乗和の逆数計算 */
static double LPC_ComputeWindowInverseSquaredSum(LPCWindowType window_type, const uint32_t num_samples)
{
switch (window_type) {
case LPC_WINDOWTYPE_RECTANGULAR:
return 1.0;
case LPC_WINDOWTYPE_WELCH:
{
const double n = num_samples - 1;
return (15 * (n - 1) * (n - 1) * (n - 1)) / (8 * n * (n - 2) * (n * n - 2 * n + 2));
}
default:
assert(0);
}

return 1.0;
}

/*（標本）自己相関の計算 */
static LPCError LPC_CalculateAutoCorrelation(
const double *data, uint32_t num_samples, double *auto_corr, uint32_t order)
Expand Down Expand Up @@ -461,6 +479,15 @@ static LPCError LPC_CalculateCoef(
return LPC_ERROR_NG;
}

/* 誤差分散で窓の影響を考慮 */
{
uint32_t i;
const double inverse_sqr = LPC_ComputeWindowInverseSquaredSum(window_type, num_samples);
for (i = 0; i < coef_order + 1; i++) {
lpcc->error_vars[i] *= inverse_sqr;
}
}

return LPC_ERROR_OK;
}

Expand Down Expand Up @@ -499,7 +526,7 @@ LPCApiResult LPCCalculator_CalculateLPCCoefficients(
/* Levinson-Durbin再帰計算により与えられた次数まで全てのLPC係数を求める（倍精度） */
LPCApiResult LPCCalculator_CalculateMultipleLPCCoefficients(
struct LPCCalculator* lpcc,
const double* data, uint32_t num_samples, double **lpc_coefs, uint32_t max_coef_order,
const double* data, uint32_t num_samples, double **lpc_coefs, double *error_vars, uint32_t max_coef_order,
LPCWindowType window_type, double regular_term)
{
uint32_t k;
Expand Down Expand Up @@ -528,36 +555,6 @@ LPCApiResult LPCCalculator_CalculateMultipleLPCCoefficients(
for (k = 0; k < max_coef_order; k++) {
memmove(lpc_coefs[k], &lpcc->a_vecs[k][1], sizeof(double) * max_coef_order);
}

return LPC_APIRESULT_OK;
}

/* Levinson-Durbin再帰計算により与えられた次数までの残差分散を求める（倍精度） */
LPCApiResult LPCCalculator_CalculateErrorVariances(
struct LPCCalculator *lpcc,
const double *data, uint32_t num_samples, double *error_vars, uint32_t max_coef_order,
LPCWindowType window_type, double regular_term)
{
/* 引数チェック */
if ((data == NULL) || (error_vars == NULL)) {
return LPC_APIRESULT_INVALID_ARGUMENT;
}

/* 次数チェック */
if (max_coef_order > lpcc->max_order) {
return LPC_APIRESULT_EXCEED_MAX_ORDER;
}

/* 入力サンプル数チェック */
if (num_samples > lpcc->max_num_buffer_samples) {
return LPC_APIRESULT_EXCEED_MAX_NUM_SAMPLES;
}

/* 係数計算 */
if (LPC_CalculateCoef(lpcc, data, num_samples, max_coef_order, window_type, regular_term) != LPC_ERROR_OK) {
return LPC_APIRESULT_FAILED_TO_CALCULATION;
}

/* 計算成功時は結果をコピー */
memmove(error_vars, lpcc->error_vars, sizeof(double) * (max_coef_order + 1));

Expand Down Expand Up @@ -1050,22 +1047,8 @@ static LPCError LPC_CalculateCoefSVR(
return LPC_ERROR_INVALID_ARGUMENT;
}

/* Levinson-Durbin法で係数を求める */
if ((err = LPC_CalculateCoef(lpcc, data, num_samples, coef_order, window_type, regular_term)) != LPC_ERROR_OK) {
return err;
}
/* 0次自己相関（信号の二乗和）が小さい場合
* => 係数は全て0として無音出力システムを予測 */
if (fabs(lpcc->auto_corr[0]) < FLT_EPSILON) {
for (i = 0; i < coef_order; i++) {
coef[i] = 0.0;
}
return LPC_ERROR_OK;
}

/* 学習しない場合はLevinson-Durbin法の結果をそのまま採用 */
/* 学習しない場合は即時終了 */
if (max_num_iteration == 0) {
memcpy(coef, &lpcc->a_vecs[coef_order - 1][1], sizeof(double) * coef_order);
return LPC_ERROR_OK;
}

Expand All @@ -1086,8 +1069,8 @@ static LPCError LPC_CalculateCoefSVR(
return LPC_ERROR_OK;
}

/* 初期値をLevinson-Durbin法の係数に設定 */
memcpy(init_coef, &lpcc->a_vecs[coef_order - 1][1], sizeof(double) * coef_order);
/* 初期値を引数の係数に設定 */
memcpy(init_coef, coef, sizeof(double) * coef_order);
memcpy(best_coef, init_coef, sizeof(double) * coef_order);

/* TODO: 係数は順序反転した方がresidualの計算が早そう（要検証） */
Expand Down Expand Up @@ -1137,6 +1120,7 @@ static LPCError LPC_CalculateCoefSVR(
}
}

/* 最良係数を記録 */
memcpy(coef, best_coef, sizeof(double) * coef_order);

return LPC_ERROR_OK;
Expand Down
111 changes: 54 additions & 57 deletions libs/srla_encoder/src/srla_encoder.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ struct SRLAEncoder {
int32_t **buffer_int; /* 信号バッファ(int) */
int32_t **residual; /* 残差信号 */
double *buffer_double; /* 信号バッファ(double) */
double *error_vars; /* 各予測係数の残差分散列 */
double **multiple_lpc_coefs; /* 各次数の予測係数 */
uint32_t *partitions_buffer; /* 最適な分割設定の記録領域 */
struct StaticHuffmanCodes param_codes; /* パラメータ符号化用Huffman符号 */
const struct SRLAParameterPreset *parameter_preset; /* パラメータプリセット */
Expand Down Expand Up @@ -499,6 +501,10 @@ int32_t SRLAEncoder_CalculateWorkSize(const struct SRLAEncoderConfig *config)
work_size += (int32_t)(config->max_num_samples_per_block * sizeof(double) + SRLA_MEMORY_ALIGNMENT);
/* 残差信号のサイズ */
work_size += (int32_t)SRLA_CALCULATE_2DIMARRAY_WORKSIZE(int32_t, config->max_num_channels, config->max_num_samples_per_block);
/* 残差分散領域のサイズ */
work_size += (int32_t)(SRLA_MEMORY_ALIGNMENT + sizeof(double) * (config->max_num_parameters + 1));
/* 係数領域のサイズ */
work_size += (int32_t)SRLA_CALCULATE_2DIMARRAY_WORKSIZE(double, config->max_num_parameters, config->max_num_parameters);
/* 分割設定記録領域のサイズ */
work_size += (int32_t)(SRLAENCODER_CALCULATE_NUM_NODES(config->max_num_samples_per_block, config->min_num_samples_per_block) * sizeof(uint32_t) + SRLA_MEMORY_ALIGNMENT);

Expand Down Expand Up @@ -620,6 +626,15 @@ struct SRLAEncoder* SRLAEncoder_Create(const struct SRLAEncoderConfig *config, v
SRLA_ALLOCATE_2DIMARRAY(encoder->residual,
work_ptr, int32_t, config->max_num_channels, config->max_num_samples_per_block);

/* 残差分散領域 */
work_ptr = (uint8_t *)SRLAUTILITY_ROUNDUP((uintptr_t)work_ptr, SRLA_MEMORY_ALIGNMENT);
encoder->error_vars = (double *)work_ptr;
work_ptr += (config->max_num_parameters + 1) * sizeof(double);

/* 前次数の係数 */
SRLA_ALLOCATE_2DIMARRAY(encoder->multiple_lpc_coefs,
work_ptr, double, config->max_num_parameters, config->max_num_parameters);

/* doubleバッファ */
work_ptr = (uint8_t *)SRLAUTILITY_ROUNDUP((uintptr_t)work_ptr, SRLA_MEMORY_ALIGNMENT);
encoder->buffer_double = (double *)work_ptr;
Expand Down Expand Up @@ -843,21 +858,15 @@ static double SRLAEncoder_CalculateRGRMeanCodeLength(double mean_abs_error, uint
return (1.0 + k1) * (1.0 - k1factor) + (1.0 + k2 + (1.0 / (1.0 - k2factor))) * k1factor;
}

static double SRLAEncoder_ComputeInverseSquaredSumWelchWindow(const uint32_t window_size)
{
const double n = window_size - 1;
return (15 * (n - 1) * (n - 1) * (n - 1)) / (8 * n * (n - 2) * (n * n - 2 * n + 2));
}

/* 最適なLPC次数の選択 */
static SRLAError SRLAEncoder_SelectBestLPCOrder(
struct SRLAEncoder *encoder,
const double *input, uint32_t num_samples, SRLAChannelLPCOrderDecisionTactics tactics,
const struct SRLAHeader *header, SRLAChannelLPCOrderDecisionTactics tactics,
const double *input, uint32_t num_samples, const double **coefs, const double *error_vars,
uint32_t max_coef_order, uint32_t *best_coef_order)
{
SRLA_ASSERT(encoder != NULL);
SRLA_ASSERT(input != NULL);
SRLA_ASSERT(input == encoder->buffer_double); /* 現状エンコーダハンドルの領域の使用を想定 */
SRLA_ASSERT(coefs != NULL);
SRLA_ASSERT(error_vars != NULL);
SRLA_ASSERT(best_coef_order != NULL);

switch (tactics) {
Expand All @@ -868,22 +877,12 @@ static SRLAError SRLAEncoder_SelectBestLPCOrder(
case SRLA_LPC_ORDER_DECISION_TACTICS_BRUTEFORCE_SEARCH:
/* 網羅探索 */
{
LPCApiResult ret;
double minlen, len, mabse;
uint32_t i, order, smpl, tmp_best_order = 0;
double coefs[SRLA_MAX_COEFFICIENT_ORDER][SRLA_MAX_COEFFICIENT_ORDER];
double *pcoefs[SRLA_MAX_COEFFICIENT_ORDER];
for (i = 0; i < SRLA_MAX_COEFFICIENT_ORDER; i++) {
pcoefs[i] = &coefs[i][0];
}
/* 次数選択のため係数計算 */
ret = LPCCalculator_CalculateMultipleLPCCoefficients(encoder->lpcc,
input, num_samples, pcoefs, max_coef_order, LPC_WINDOWTYPE_WELCH, SRLA_LPC_RIDGE_REGULARIZATION_PARAMETER);
SRLA_ASSERT(ret == LPC_APIRESULT_OK);

minlen = FLT_MAX;
for (order = 1; order <= max_coef_order; order++) {
const double *coef = pcoefs[order - 1];
const double *coef = coefs[order - 1];
mabse = 0.0;
for (smpl = order; smpl < num_samples; smpl++) {
double residual = input[smpl];
Expand All @@ -893,7 +892,7 @@ static SRLAError SRLAEncoder_SelectBestLPCOrder(
mabse += SRLAUTILITY_ABS(residual);
}
/* 残差符号のサイズ */
len = SRLAEncoder_CalculateRGRMeanCodeLength(mabse / num_samples, encoder->header.bits_per_sample) * num_samples;
len = SRLAEncoder_CalculateRGRMeanCodeLength(mabse / num_samples, header->bits_per_sample) * num_samples;
/* 係数のサイズ */
len += SRLA_LPC_COEFFICIENT_BITWIDTH * order;
if (minlen > len) {
Expand All @@ -907,40 +906,25 @@ static SRLAError SRLAEncoder_SelectBestLPCOrder(
return SRLA_ERROR_OK;
}
case SRLA_LPC_ORDER_DECISION_TACTICS_BRUTEFORCE_ESTIMATION:
/* 最小推定符号長を与える係数次数の探索 */
{
LPCApiResult ret;
uint32_t i, tmp_best_order = 0;
double error_vars[SRLA_MAX_COEFFICIENT_ORDER + 1];
const double var_norm_factor = SRLAEncoder_ComputeInverseSquaredSumWelchWindow(num_samples);

/* 残差分散の計算 */
ret = LPCCalculator_CalculateErrorVariances(encoder->lpcc,
input, num_samples, error_vars, max_coef_order, LPC_WINDOWTYPE_WELCH, SRLA_LPC_RIDGE_REGULARIZATION_PARAMETER);
SRLA_ASSERT(ret == LPC_APIRESULT_OK);

for (i = 0; i <= max_coef_order; i++) {
error_vars[i] *= var_norm_factor;
}
uint32_t order, tmp_best_order = 0;
double len, mabse, minlen = FLT_MAX;

/* 最小推定符号長を与える係数次数の探索 */
{
uint32_t order;
double len, mabse, minlen = FLT_MAX;
/* 次数あたり係数量子化誤差分散 */
const double qerr_var_unit = pow(2.0, -2.0 * SRLA_LPC_COEFFICIENT_BITWIDTH) * error_vars[0] / 12.0;
for (order = 1; order <= max_coef_order; order++) {
/* 係数量子化誤差分散を加えた誤差分散 */
const double err_var = error_vars[order] + order * qerr_var_unit;
/* Laplace分布の仮定で残差分散から平均絶対値を推定 */
mabse = 2.0 * sqrt(err_var / 2.0); /* 符号化で非負整数化するため2倍 */
/* 残差符号のサイズ */
len = SRLAEncoder_CalculateRGRMeanCodeLength(mabse, encoder->header.bits_per_sample) * num_samples;
/* 係数のサイズ */
len += SRLA_LPC_COEFFICIENT_BITWIDTH * order;
if (minlen > len) {
minlen = len;
tmp_best_order = order;
}
/* 次数あたり係数量子化誤差分散 */
const double qerr_var_unit = pow(2.0, -2.0 * SRLA_LPC_COEFFICIENT_BITWIDTH) * error_vars[0] / 12.0;
for (order = 1; order <= max_coef_order; order++) {
/* 係数量子化誤差分散を加えた誤差分散 */
const double err_var = error_vars[order] + order * qerr_var_unit;
/* Laplace分布の仮定で残差分散から平均絶対値を推定 */
mabse = 2.0 * sqrt(err_var / 2.0); /* 符号化で非負整数化するため2倍 */
/* 残差符号のサイズ */
len = SRLAEncoder_CalculateRGRMeanCodeLength(mabse, header->bits_per_sample) * num_samples;
/* 係数のサイズ */
len += SRLA_LPC_COEFFICIENT_BITWIDTH * order;
if (minlen > len) {
minlen = len;
tmp_best_order = order;
}
}

Expand Down Expand Up @@ -1090,16 +1074,28 @@ static SRLAApiResult SRLAEncoder_EncodeCompressData(
for (ch = 0; ch < header->num_channels; ch++) {
uint32_t smpl, p;
LPCApiResult ret;
SRLAError err;
/* double精度の信号に変換（[-1,1]の範囲に正規化） */
const double norm_const = pow(2.0, -(int32_t)(header->bits_per_sample - 1));
for (smpl = 0; smpl < num_samples; smpl++) {
encoder->buffer_double[smpl] = encoder->buffer_int[ch][smpl] * norm_const;
}
/* 最大次数まで係数と誤差分散を計算 */
ret = LPCCalculator_CalculateMultipleLPCCoefficients(encoder->lpcc,
encoder->buffer_double, num_samples,
encoder->multiple_lpc_coefs, encoder->error_vars, encoder->parameter_preset->max_num_parameters,
LPC_WINDOWTYPE_WELCH, SRLA_LPC_RIDGE_REGULARIZATION_PARAMETER);
SRLA_ASSERT(ret == LPC_APIRESULT_OK);
/* 次数選択 */
SRLAEncoder_SelectBestLPCOrder(encoder,
encoder->buffer_double, num_samples, encoder->parameter_preset->lpc_order_tactics,
err = SRLAEncoder_SelectBestLPCOrder(&encoder->header,
encoder->parameter_preset->lpc_order_tactics,
encoder->buffer_double, num_samples, (const double **)encoder->multiple_lpc_coefs, encoder->error_vars,
encoder->parameter_preset->max_num_parameters, &encoder->coef_order[ch]);
/* LPC係数計算 */
SRLA_ASSERT(err == SRLA_ERROR_OK);
/* 最良の次数をパラメータに設定 */
memcpy(encoder->params_double[ch],
encoder->multiple_lpc_coefs[encoder->coef_order[ch] - 1], sizeof(double) * encoder->coef_order[ch]);
/* SVRによるLPC係数計算 */
ret = LPCCalculator_CalculateLPCCoefficientsSVR(encoder->lpcc,
encoder->buffer_double, num_samples,
encoder->params_double[ch], encoder->coef_order[ch], encoder->parameter_preset->svr_max_num_iterations,
Expand All @@ -1111,6 +1107,7 @@ static SRLAApiResult SRLAEncoder_EncodeCompressData(
encoder->params_double[ch][p] = encoder->params_double[ch][encoder->coef_order[ch] - p - 1];
encoder->params_double[ch][encoder->coef_order[ch] - p - 1] = tmp;
}
/* LPC係数量子化 */
ret = LPC_QuantizeCoefficients(encoder->params_double[ch], encoder->coef_order[ch],
SRLA_LPC_COEFFICIENT_BITWIDTH, (1 << SRLA_RSHIFT_LPC_COEFFICIENT_BITWIDTH),
encoder->params_int[ch], &encoder->rshifts[ch]);
Expand Down

0 comments on commit 7394a2c

Please sign in to comment.