LPC（Linear Predictive Coding, 線形予測符号化）

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>

/* （標本）自己相関の計算 */
static int
calc_auto_correlation(double *auto_corr, const double *data, const size_t num_sample, const size_t max_order);

/* Levinson-Durbin再帰計算 */
static int
levinson_durbin_recursion(double *lpc_coef, const double *auto_corr, const size_t max_order);

int main(void)
{
  int num_sample = 300; /* サンプル数 */
  int max_delay  = 10;  /* LPC係数の数 */

  int i_sample, i_delay;
  double *data     = (double *)malloc(sizeof(double) * num_sample);
  double *predict  = (double *)malloc(sizeof(double) * num_sample);
  double *auto_cor = (double *)malloc(sizeof(double) * (max_delay + 1));
  double *coff     = (double *)malloc(sizeof(double) * (max_delay + 1));
  double error;

  /* 波形の生成 */
  for (i_sample = 0; i_sample < num_sample; i_sample++) {
    data[i_sample] = sin(i_sample * 0.01) + 0.5 * cos(4.0f * sin(i_sample * 0.05));
  }

  /* 自己相関・Levinson-Durbin再帰計算 */
  calc_auto_correlation(auto_cor, data, num_sample, max_delay + 1);
  levinson_durbin_recursion(coff, auto_cor, max_delay);

  /* 予測テスト */
  for (i_sample = 0; i_sample < num_sample; i_sample++) {
    if (i_sample < max_delay) {
      /* 最初のmax_delayステップ分は元信号を単純コピー */
      predict[i_sample] = data[i_sample];
    } else {
      /* 以降は予測 */
      predict[i_sample] = 0.0f;
      for (i_delay = 1; i_delay <= max_delay; i_delay++) {
        predict[i_sample] -= (coff[i_delay] * data[i_sample-i_delay]);
      }
    }
  }

  /* 誤差計算・結果表示 */
  error = 0.0f;
  for (i_sample = 0; i_sample < num_sample; i_sample++) {
    error += pow(predict[i_sample] - data[i_sample], 2);
    printf("No:%d Data: %f Predict: %f \n", i_sample, data[i_sample], predict[i_sample]);
  }
  printf("Error : %f \n", sqrt(error / num_sample));

  free(data); free(predict);
  free(auto_cor); free(coff);

  return 0;
}

static int levinson_durbin_recursion(double *lpc_coef, const double *auto_corr, const size_t max_order)
{
  int delay, i_delay;
  double lambda;
  double *u_vec, *v_vec, *a_vec, *e_vec;

  if (lpc_coef == NULL || auto_corr == NULL) {
    fprintf(stderr, "Data or result pointer point to NULL. \n");
    return -1;
  }

  /*
   * 0次自己相関（信号の二乗和）が0に近い場合、入力信号は無音と判定
   * => 予測誤差, LPC係数は全て0として無音出力システムを予測.
   */
  if (fabs(auto_corr[0]) < FLT_EPSILON) {
    for (i_delay = 0; i_delay < max_order+1; ++i_delay) {
      lpc_coef[i_delay] = 0.0f;
    }
    return 0;
  }

  /* 初期化 */
  a_vec = (double *)malloc(sizeof(double) * (max_order + 2)); /* a_0, a_k+1を含めるとmax_order+2 */
  e_vec = (double *)malloc(sizeof(double) * (max_order + 2)); /* e_0, e_k+1を含めるとmax_order+2 */
  u_vec = (double *)malloc(sizeof(double) * (max_order + 2));
  v_vec = (double *)malloc(sizeof(double) * (max_order + 2));
  for (i_delay = 0; i_delay < max_order + 2; i_delay++) {
    u_vec[i_delay] = v_vec[i_delay] = a_vec[i_delay] = 0.0f;
  }

  /* 最初のステップの係数をセット */
  a_vec[0] = 1.0f;
  e_vec[0] = auto_corr[0];
  a_vec[1] = - auto_corr[1] / auto_corr[0];
  e_vec[1] = auto_corr[0] + auto_corr[1] * a_vec[1];
  u_vec[0] = 1.0f; u_vec[1] = 0.0f;
  v_vec[0] = 0.0f; v_vec[1] = 1.0f;

  /* 再帰処理 */
  for (delay = 1; delay < max_order; delay++) {
    lambda = 0.0f;
    for (i_delay = 0; i_delay < delay+1; i_delay++) {
      lambda += a_vec[i_delay] * auto_corr[delay+1-i_delay];
    }
    lambda /= (-e_vec[delay]);
    e_vec[delay+1] = (1 - lambda * lambda) * e_vec[delay];

    /* u_vec, v_vecの更新 */
    for (i_delay = 0; i_delay < delay; i_delay++) {
      u_vec[i_delay+1] = v_vec[delay-i_delay] = a_vec[i_delay+1];
    }
    u_vec[0] = 1.0f; u_vec[delay+1] = 0.0f;
    v_vec[0] = 0.0f; v_vec[delay+1] = 1.0f;

    /* resultの更新 */
    for (i_delay = 0; i_delay < delay+2; i_delay++) {
       a_vec[i_delay] = u_vec[i_delay] + lambda * v_vec[i_delay];
    }

  }

  /* 結果の取得 */
  memcpy(lpc_coef, a_vec, sizeof(double) * (max_order + 1));

  free(u_vec); free(v_vec);
  free(a_vec); free(e_vec);

  return 0;
}

static int calc_auto_correlation(double *auto_corr, const double *data, const size_t num_sample, const size_t max_order)
{
  int i_sample, delay_time;

  if (max_order > num_sample) {
    fprintf(stderr, "Max order(%zu) is larger than number of samples(%zu). \n", max_order, num_sample);
    return -1;
  }

  if (auto_corr == NULL || data == NULL) {
    fprintf(stderr, "Data or result pointer point to NULL. \n");
    return -2;
  }

  /* （標本）自己相関の計算 */
  for (delay_time = 0; delay_time < max_order; delay_time++) {
    auto_corr[delay_time] = 0.0f;
    for (i_sample = delay_time; i_sample < num_sample; i_sample++) {
      auto_corr[delay_time] += data[i_sample] * data[i_sample - delay_time];
    }
  }

  return 0;
}