MLSA digital filter for speech synthesis in C++

mlsa_filter.h

#pragma once

#include <cmath>
#include <memory>
#include <vector>
#include <cassert>

namespace sp {

/**
 * MLSA BASE digital filter (Mel-log Spectrum Approximate digital filter)
 */
class mlsa_base_filter {
public:
  mlsa_base_filter(const int order, const double alpha);

  template <class Vector>
  double filter(const double x, const Vector& b);

 private:
  mlsa_base_filter();

  double alpha_;
  std::vector<double> delay_;
};

mlsa_base_filter::mlsa_base_filter(const int order, const double alpha)
: alpha_(alpha),
  delay_(order+1)
{
}

template <class Vector>
double mlsa_base_filter::filter(const double x, const Vector& b)
{
  double result = 0.0;
  
  delay_[0] = x;
  delay_[1] = (1.0-alpha_*alpha_)*delay_[0] + alpha_*delay_[1];

  for (size_t i = 2; i < b.size(); ++i) {
    delay_[i] = delay_[i] + alpha_*(delay_[i+1]-delay_[i-1]);
    result += delay_[i] * b[i];
  }

  // special case
  // TODO: other solution?
  if (b.size() == 2) {
    result += delay_[1] * b[1];
  } 
  
  // t <- t+1 in time
  for (size_t i = delay_.size()-1; i > 1; --i) {
    delay_[i] = delay_[i-1];
  }
  
  return result;
}

/**
 * MLSA digital filter cascaded
 */
class mlsa_base_cascaded_filter {
 public:
  mlsa_base_cascaded_filter(const int order,
			    const double alpha,
			    const int n_pade);

  template <class Vector>
  double filter(const double x, const Vector& b);
  
 private:
  mlsa_base_cascaded_filter();

  std::vector<std::unique_ptr<mlsa_base_filter>> base_f_; // cascadad filters
  std::vector<double> delay_;
  std::vector<double> pade_coef_;
};

mlsa_base_cascaded_filter::mlsa_base_cascaded_filter(const int order, 
						     const double alpha,
						     const int n_pade)
  : delay_(n_pade + 1),
  pade_coef_(n_pade + 1)
{ 
  using std::unique_ptr;

  if (n_pade != 4 && n_pade != 5) {
    std::cerr << "The number of pade approximations must be 4 or 5."
	      << std::endl;
  }
  assert(n_pade == 4 || n_pade == 5);

  for (int i = 0; i <= n_pade; ++i) {
    mlsa_base_filter* p = new mlsa_base_filter(order, alpha);
    base_f_.push_back(unique_ptr<mlsa_base_filter>(p));
  }
  
  if (n_pade == 4) {
    pade_coef_[0] = 1.0;
    pade_coef_[1] = 4.999273e-1;
    pade_coef_[2] = 1.067005e-1;
    pade_coef_[3] = 1.170221e-2;
    pade_coef_[4] = 5.656279e-4;
  }
  
  if (n_pade == 5) {
    pade_coef_[0] = 1.0;
    pade_coef_[1] = 4.999391e-1;
    pade_coef_[2] = 1.107098e-1;
    pade_coef_[3] = 1.369984e-2;
    pade_coef_[4] = 9.564853e-4;
    pade_coef_[5] = 3.041721e-5;
  }   
}

template <class Vector>
double mlsa_base_cascaded_filter::filter(const double x, const Vector& b)
{
  double result = 0.0;  
  double feed_back = 0.0;

  for (size_t i = pade_coef_.size()-1; i >= 1; --i) {
    delay_[i] = base_f_[i]->filter(delay_[i-1], b);
    double v = delay_[i] * pade_coef_[i];
    if (i % 2 == 1) {
      feed_back += v;
    } else {
      feed_back -= v;
    }
    result += v;
  }

  delay_[0] = feed_back + x;
  result += delay_[0];

  return result;
}

/**
 * MLSA digital filter (Mel-log Spectrum Approximate digital filter)
 * The filter consists of two stage cascade filters
 */
class mlsa_filter {
 public:
  mlsa_filter(const int order, const double alpha, const int n_pade);
 ~mlsa_filter();

 template <class Vector>
 double filter(const double x, const Vector& b);

 private:
 mlsa_filter();

  double alpha_;
  std::unique_ptr<mlsa_base_cascaded_filter> f1_; // first stage
  std::unique_ptr<mlsa_base_cascaded_filter> f2_; // second stage
};

mlsa_filter::mlsa_filter(const int order,
			 const double alpha,
			 const int n_pade)
  : alpha_(alpha),
  f1_(new mlsa_base_cascaded_filter(2, alpha, n_pade)),
  f2_(new mlsa_base_cascaded_filter(order, alpha, n_pade))
{
}

mlsa_filter::~mlsa_filter()
{
}

template <class Vector>
double mlsa_filter::filter(const double x, const Vector& b)
{
  // 1. First stage filtering
  Vector b1 = {0, b[1]};
  double y = f1_->filter(x, b1);
  
  // 2. Second stage filtering
  double result = f2_->filter(y, b);

  return result;
}

} // end namespace sp