/* ----------------------------------------------------------------- */
/*           The HMM-Based Speech Synthesis Engine "hts_engine API"  */
/*           developed by HTS Working Group                          */
/*           http://hts-engine.sourceforge.net/                      */
/* ----------------------------------------------------------------- */
/*                                                                   */
/*  Copyright (c) 2001-2013  Nagoya Institute of Technology          */
/*                           Department of Computer Science          */
/*                                                                   */
/*                2001-2008  Tokyo Institute of Technology           */
/*                           Interdisciplinary Graduate School of    */
/*                           Science and Engineering                 */
/*                                                                   */
/* All rights reserved.                                              */
/*                                                                   */
/* Redistribution and use in source and binary forms, with or        */
/* without modification, are permitted provided that the following   */
/* conditions are met:                                               */
/*                                                                   */
/* - Redistributions of source code must retain the above copyright  */
/*   notice, this list of conditions and the following disclaimer.   */
/* - Redistributions in binary form must reproduce the above         */
/*   copyright notice, this list of conditions and the following     */
/*   disclaimer in the documentation and/or other materials provided */
/*   with the distribution.                                          */
/* - Neither the name of the HTS working group nor the names of its  */
/*   contributors may be used to endorse or promote products derived */
/*   from this software without specific prior written permission.   */
/*                                                                   */
/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND            */
/* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,       */
/* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF          */
/* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE          */
/* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS */
/* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,          */
/* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED   */
/* TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,     */
/* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON */
/* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,   */
/* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    */
/* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           */
/* POSSIBILITY OF SUCH DAMAGE.                                       */
/* ----------------------------------------------------------------- */

#ifndef HTS_VOCODER_C
#define HTS_VOCODER_C

#ifdef __cplusplus
#define HTS_VOCODER_C_START extern "C" {
#define HTS_VOCODER_C_END   }
#else
#define HTS_VOCODER_C_START
#define HTS_VOCODER_C_END
#endif                          /* __CPLUSPLUS */

HTS_VOCODER_C_START;

#include <math.h>               /* for sqrt(),log(),exp(),pow(),cos() */

/* hts_engine libraries */
#include "HTS_hidden.h"

static const double HTS_pade[21] = {
   1.00000000000,
   1.00000000000,
   0.00000000000,
   1.00000000000,
   0.00000000000,
   0.00000000000,
   1.00000000000,
   0.00000000000,
   0.00000000000,
   0.00000000000,
   1.00000000000,
   0.49992730000,
   0.10670050000,
   0.01170221000,
   0.00056562790,
   1.00000000000,
   0.49993910000,
   0.11070980000,
   0.01369984000,
   0.00095648530,
   0.00003041721
};

/* HTS_movem: move memory */
static void HTS_movem(double *a, double *b, const int nitem)
{
   long i = (long) nitem;

   if (a > b)
      while (i--)
         *b++ = *a++;
   else {
      a += i;
      b += i;
      while (i--)
         *--b = *--a;
   }
}

/* HTS_mlsafir: sub functions for MLSA filter */
static double HTS_mlsafir(const double x, const double *b, const int m, const double a, const double aa, double *d)
{
   double y = 0.0;
   int i;

   d[0] = x;
   d[1] = aa * d[0] + a * d[1];

   for (i = 2; i <= m; i++)
      d[i] += a * (d[i + 1] - d[i - 1]);

   for (i = 2; i <= m; i++)
      y += d[i] * b[i];

   for (i = m + 1; i > 1; i--)
      d[i] = d[i - 1];

   return (y);
}

/* HTS_mlsadf1: sub functions for MLSA filter */
static double HTS_mlsadf1(double x, const double *b, const int m, const double a, const double aa, const int pd, double *d, const double *ppade)
{
   double v, out = 0.0, *pt;
   int i;

   pt = &d[pd + 1];

   for (i = pd; i >= 1; i--) {
      d[i] = aa * pt[i - 1] + a * d[i];
      pt[i] = d[i] * b[1];
      v = pt[i] * ppade[i];
      x += (1 & i) ? v : -v;
      out += v;
   }

   pt[0] = x;
   out += x;

   return (out);
}

/* HTS_mlsadf2: sub functions for MLSA filter */
static double HTS_mlsadf2(double x, const double *b, const int m, const double a, const double aa, const int pd, double *d, const double *ppade)
{
   double v, out = 0.0, *pt;
   int i;

   pt = &d[pd * (m + 2)];

   for (i = pd; i >= 1; i--) {
      pt[i] = HTS_mlsafir(pt[i - 1], b, m, a, aa, &d[(i - 1) * (m + 2)]);
      v = pt[i] * ppade[i];

      x += (1 & i) ? v : -v;
      out += v;
   }

   pt[0] = x;
   out += x;

   return (out);
}

/* HTS_mlsadf: functions for MLSA filter */
static double HTS_mlsadf(double x, const double *b, const int m, const double a, const int pd, double *d)
{
   const double aa = 1 - a * a;
   const double *ppade = &(HTS_pade[pd * (pd + 1) / 2]);

   x = HTS_mlsadf1(x, b, m, a, aa, pd, d, ppade);
   x = HTS_mlsadf2(x, b, m, a, aa, pd, &d[2 * (pd + 1)], ppade);

   return (x);
}

/* HTS_rnd: functions for random noise generation */
static double HTS_rnd(unsigned long *next)
{
   double r;

   *next = *next * 1103515245L + 12345;
   r = (*next / 65536L) % 32768L;

   return (r / RANDMAX);
}

/* HTS_nrandom: functions for gaussian random noise generation */
static double HTS_nrandom(HTS_Vocoder * v)
{
   if (v->sw == 0) {
      v->sw = 1;
      do {
         v->r1 = 2 * HTS_rnd(&v->next) - 1;
         v->r2 = 2 * HTS_rnd(&v->next) - 1;
         v->s = v->r1 * v->r1 + v->r2 * v->r2;
      } while (v->s > 1 || v->s == 0);
      v->s = sqrt(-2 * log(v->s) / v->s);
      return (v->r1 * v->s);
   } else {
      v->sw = 0;
      return (v->r2 * v->s);
   }
}

/* HTS_mceq: function for M-sequence random noise generation */
static int HTS_mseq(HTS_Vocoder * v)
{
   int x0, x28;

   v->x >>= 1;
   if (v->x & B0)
      x0 = 1;
   else
      x0 = -1;
   if (v->x & B28)
      x28 = 1;
   else
      x28 = -1;
   if (x0 + x28)
      v->x &= B31_;
   else
      v->x |= B31;

   return (x0);
}

/* HTS_mc2b: transform mel-cepstrum to MLSA digital fillter coefficients */
static void HTS_mc2b(double *mc, double *b, int m, const double a)
{
   if (mc != b) {
      if (a != 0.0) {
         b[m] = mc[m];
         for (m--; m >= 0; m--)
            b[m] = mc[m] - a * b[m + 1];
      } else
         HTS_movem(mc, b, m + 1);
   } else if (a != 0.0)
      for (m--; m >= 0; m--)
         b[m] -= a * b[m + 1];
}

/* HTS_b2bc: transform MLSA digital filter coefficients to mel-cepstrum */
static void HTS_b2mc(const double *b, double *mc, int m, const double a)
{
   double d, o;

   d = mc[m] = b[m];
   for (m--; m >= 0; m--) {
      o = b[m] + a * d;
      d = b[m];
      mc[m] = o;
   }
}

/* HTS_freqt: frequency transformation */
static void HTS_freqt(HTS_Vocoder * v, const double *c1, const int m1, double *c2, const int m2, const double a)
{
   int i, j;
   const double b = 1 - a * a;
   double *g;

   if (m2 > v->freqt_size) {
      if (v->freqt_buff != NULL)
         HTS_free(v->freqt_buff);
      v->freqt_buff = (double *) HTS_calloc(m2 + m2 + 2, sizeof(double));
      v->freqt_size = m2;
   }
   g = v->freqt_buff + v->freqt_size + 1;

   for (i = 0; i < m2 + 1; i++)
      g[i] = 0.0;

   for (i = -m1; i <= 0; i++) {
      if (0 <= m2)
         g[0] = c1[-i] + a * (v->freqt_buff[0] = g[0]);
      if (1 <= m2)
         g[1] = b * v->freqt_buff[0] + a * (v->freqt_buff[1] = g[1]);
      for (j = 2; j <= m2; j++)
         g[j] = v->freqt_buff[j - 1] + a * ((v->freqt_buff[j] = g[j]) - g[j - 1]);
   }

   HTS_movem(g, c2, m2 + 1);
}

/* HTS_c2ir: The minimum phase impulse response is evaluated from the minimum phase cepstrum */
static void HTS_c2ir(const double *c, const int nc, double *h, const int leng)
{
   int n, k, upl;
   double d;

   h[0] = exp(c[0]);
   for (n = 1; n < leng; n++) {
      d = 0;
      upl = (n >= nc) ? nc - 1 : n;
      for (k = 1; k <= upl; k++)
         d += k * c[k] * h[n - k];
      h[n] = d / n;
   }
}

/* HTS_b2en: calculate frame energy */
static double HTS_b2en(HTS_Vocoder * v, const double *b, const int m, const double a)
{
   int i;
   double en = 0.0;
   double *cep;
   double *ir;

   if (v->spectrum2en_size < m) {
      if (v->spectrum2en_buff != NULL)
         HTS_free(v->spectrum2en_buff);
      v->spectrum2en_buff = (double *) HTS_calloc((m + 1) + 2 * IRLENG, sizeof(double));
      v->spectrum2en_size = m;
   }
   cep = v->spectrum2en_buff + m + 1;
   ir = cep + IRLENG;

   HTS_b2mc(b, v->spectrum2en_buff, m, a);
   HTS_freqt(v, v->spectrum2en_buff, m, cep, IRLENG - 1, -a);
   HTS_c2ir(cep, IRLENG, ir, IRLENG);

   for (i = 0; i < IRLENG; i++)
      en += ir[i] * ir[i];

   return (en);
}

/* HTS_ignorm: inverse gain normalization */
static void HTS_ignorm(double *c1, double *c2, int m, const double g)
{
   double k;
   if (g != 0.0) {
      k = pow(c1[0], g);
      for (; m >= 1; m--)
         c2[m] = k * c1[m];
      c2[0] = (k - 1.0) / g;
   } else {
      HTS_movem(&c1[1], &c2[1], m);
      c2[0] = log(c1[0]);
   }
}

/* HTS_gnorm: gain normalization */
static void HTS_gnorm(double *c1, double *c2, int m, const double g)
{
   double k;
   if (g != 0.0) {
      k = 1.0 + g * c1[0];
      for (; m >= 1; m--)
         c2[m] = c1[m] / k;
      c2[0] = pow(k, 1.0 / g);
   } else {
      HTS_movem(&c1[1], &c2[1], m);
      c2[0] = exp(c1[0]);
   }
}

/* HTS_lsp2lpc: transform LSP to LPC */
static void HTS_lsp2lpc(HTS_Vocoder * v, double *lsp, double *a, const int m)
{
   int i, k, mh1, mh2, flag_odd;
   double xx, xf, xff;
   double *p, *q;
   double *a0, *a1, *a2, *b0, *b1, *b2;

   flag_odd = 0;
   if (m % 2 == 0)
      mh1 = mh2 = m / 2;
   else {
      mh1 = (m + 1) / 2;
      mh2 = (m - 1) / 2;
      flag_odd = 1;
   }

   if (m > v->lsp2lpc_size) {
      if (v->lsp2lpc_buff != NULL)
         HTS_free(v->lsp2lpc_buff);
      v->lsp2lpc_buff = (double *) HTS_calloc(5 * m + 6, sizeof(double));
      v->lsp2lpc_size = m;
   }
   p = v->lsp2lpc_buff + m;
   q = p + mh1;
   a0 = q + mh2;
   a1 = a0 + (mh1 + 1);
   a2 = a1 + (mh1 + 1);
   b0 = a2 + (mh1 + 1);
   b1 = b0 + (mh2 + 1);
   b2 = b1 + (mh2 + 1);

   HTS_movem(lsp, v->lsp2lpc_buff, m);

   for (i = 0; i < mh1 + 1; i++)
      a0[i] = 0.0;
   for (i = 0; i < mh1 + 1; i++)
      a1[i] = 0.0;
   for (i = 0; i < mh1 + 1; i++)
      a2[i] = 0.0;
   for (i = 0; i < mh2 + 1; i++)
      b0[i] = 0.0;
   for (i = 0; i < mh2 + 1; i++)
      b1[i] = 0.0;
   for (i = 0; i < mh2 + 1; i++)
      b2[i] = 0.0;

   /* lsp filter parameters */
   for (i = k = 0; i < mh1; i++, k += 2)
      p[i] = -2.0 * cos(v->lsp2lpc_buff[k]);
   for (i = k = 0; i < mh2; i++, k += 2)
      q[i] = -2.0 * cos(v->lsp2lpc_buff[k + 1]);

   /* impulse response of analysis filter */
   xx = 1.0;
   xf = xff = 0.0;

   for (k = 0; k <= m; k++) {
      if (flag_odd) {
         a0[0] = xx;
         b0[0] = xx - xff;
         xff = xf;
         xf = xx;
      } else {
         a0[0] = xx + xf;
         b0[0] = xx - xf;
         xf = xx;
      }

      for (i = 0; i < mh1; i++) {
         a0[i + 1] = a0[i] + p[i] * a1[i] + a2[i];
         a2[i] = a1[i];
         a1[i] = a0[i];
      }

      for (i = 0; i < mh2; i++) {
         b0[i + 1] = b0[i] + q[i] * b1[i] + b2[i];
         b2[i] = b1[i];
         b1[i] = b0[i];
      }

      if (k != 0)
         a[k - 1] = -0.5 * (a0[mh1] + b0[mh2]);
      xx = 0.0;
   }

   for (i = m - 1; i >= 0; i--)
      a[i + 1] = -a[i];
   a[0] = 1.0;
}

/* HTS_gc2gc: generalized cepstral transformation */
static void HTS_gc2gc(HTS_Vocoder * v, double *c1, const int m1, const double g1, double *c2, const int m2, const double g2)
{
   int i, min, k, mk;
   double ss1, ss2, cc;

   if (m1 > v->gc2gc_size) {
      if (v->gc2gc_buff != NULL)
         HTS_free(v->gc2gc_buff);
      v->gc2gc_buff = (double *) HTS_calloc(m1 + 1, sizeof(double));
      v->gc2gc_size = m1;
   }

   HTS_movem(c1, v->gc2gc_buff, m1 + 1);

   c2[0] = v->gc2gc_buff[0];
   for (i = 1; i <= m2; i++) {
      ss1 = ss2 = 0.0;
      min = m1 < i ? m1 : i - 1;
      for (k = 1; k <= min; k++) {
         mk = i - k;
         cc = v->gc2gc_buff[k] * c2[mk];
         ss2 += k * cc;
         ss1 += mk * cc;
      }

      if (i <= m1)
         c2[i] = v->gc2gc_buff[i] + (g2 * ss2 - g1 * ss1) / i;
      else
         c2[i] = (g2 * ss2 - g1 * ss1) / i;
   }
}

/* HTS_mgc2mgc: frequency and generalized cepstral transformation */
static void HTS_mgc2mgc(HTS_Vocoder * v, double *c1, const int m1, const double a1, const double g1, double *c2, const int m2, const double a2, const double g2)
{
   double a;

   if (a1 == a2) {
      HTS_gnorm(c1, c1, m1, g1);
      HTS_gc2gc(v, c1, m1, g1, c2, m2, g2);
      HTS_ignorm(c2, c2, m2, g2);
   } else {
      a = (a2 - a1) / (1 - a1 * a2);
      HTS_freqt(v, c1, m1, c2, m2, a);
      HTS_gnorm(c2, c2, m2, g1);
      HTS_gc2gc(v, c2, m2, g1, c2, m2, g2);
      HTS_ignorm(c2, c2, m2, g2);
   }
}

/* HTS_lsp2mgc: transform LSP to MGC */
static void HTS_lsp2mgc(HTS_Vocoder * v, double *lsp, double *mgc, const int m, const double alpha)
{
   int i;
   /* lsp2lpc */
   HTS_lsp2lpc(v, lsp + 1, mgc, m);
   if (v->use_log_gain)
      mgc[0] = exp(lsp[0]);
   else
      mgc[0] = lsp[0];

   /* mgc2mgc */
   if (NORMFLG1)
      HTS_ignorm(mgc, mgc, m, v->gamma);
   else if (MULGFLG1)
      mgc[0] = (1.0 - mgc[0]) * ((double) v->stage);
   if (MULGFLG1)
      for (i = m; i >= 1; i--)
         mgc[i] *= -((double) v->stage);
   HTS_mgc2mgc(v, mgc, m, alpha, v->gamma, mgc, m, alpha, v->gamma);
   if (NORMFLG2)
      HTS_gnorm(mgc, mgc, m, v->gamma);
   else if (MULGFLG2)
      mgc[0] = mgc[0] * v->gamma + 1.0;
   if (MULGFLG2)
      for (i = m; i >= 1; i--)
         mgc[i] *= v->gamma;
}

/* HTS_mglsadff: sub functions for MGLSA filter */
static double HTS_mglsadff(double x, const double *b, const int m, const double a, double *d)
{
   int i;

   double y;
   y = d[0] * b[1];
   for (i = 1; i < m; i++) {
      d[i] += a * (d[i + 1] - d[i - 1]);
      y += d[i] * b[i + 1];
   }
   x -= y;

   for (i = m; i > 0; i--)
      d[i] = d[i - 1];
   d[0] = a * d[0] + (1 - a * a) * x;
   return x;
}

/* HTS_mglsadf: sub functions for MGLSA filter */
static double HTS_mglsadf(double x, const double *b, const int m, const double a, const int n, double *d)
{
   int i;

   for (i = 0; i < n; i++)
      x = HTS_mglsadff(x, b, m, a, &d[i * (m + 1)]);

   return x;
}

/* THS_check_lsp_stability: check LSP stability */
static void HTS_check_lsp_stability(double *lsp, size_t m)
{
   size_t i, j;
   double tmp;
   double min = (CHECK_LSP_STABILITY_MIN * PI) / (m + 1);
   HTS_Boolean find;

   for (i = 0; i < CHECK_LSP_STABILITY_NUM; i++) {
      find = FALSE;

      for (j = 1; j < m; j++) {
         tmp = lsp[j + 1] - lsp[j];
         if (tmp < min) {
            lsp[j] -= 0.5 * (min - tmp);
            lsp[j + 1] += 0.5 * (min - tmp);
            find = TRUE;
         }
      }

      if (lsp[1] < min) {
         lsp[1] = min;
         find = TRUE;
      }
      if (lsp[m] > PI - min) {
         lsp[m] = PI - min;
         find = TRUE;
      }

      if (find == FALSE)
         break;
   }
}

/* HTS_lsp2en: calculate frame energy */
static double HTS_lsp2en(HTS_Vocoder * v, double *lsp, size_t m, double alpha)
{
   size_t i;
   double en = 0.0;
   double *buff;

   if (v->spectrum2en_size < m) {
      if (v->spectrum2en_buff != NULL)
         HTS_free(v->spectrum2en_buff);
      v->spectrum2en_buff = (double *) HTS_calloc(m + 1 + IRLENG, sizeof(double));
      v->spectrum2en_size = m;
   }
   buff = v->spectrum2en_buff + m + 1;

   /* lsp2lpc */
   HTS_lsp2lpc(v, lsp + 1, v->spectrum2en_buff, m);
   if (v->use_log_gain)
      v->spectrum2en_buff[0] = exp(lsp[0]);
   else
      v->spectrum2en_buff[0] = lsp[0];

   /* mgc2mgc */
   if (NORMFLG1)
      HTS_ignorm(v->spectrum2en_buff, v->spectrum2en_buff, m, v->gamma);
   else if (MULGFLG1)
      v->spectrum2en_buff[0] = (1.0 - v->spectrum2en_buff[0]) * ((double) v->stage);
   if (MULGFLG1)
      for (i = 1; i <= m; i++)
         v->spectrum2en_buff[i] *= -((double) v->stage);
   HTS_mgc2mgc(v, v->spectrum2en_buff, m, alpha, v->gamma, buff, IRLENG - 1, 0.0, 1);

   for (i = 0; i < IRLENG; i++)
      en += buff[i] * buff[i];
   return en;
}

/* HTS_white_noise: return white noise */
static double HTS_white_noise(HTS_Vocoder * v)
{
   if (v->gauss)
      return (double) HTS_nrandom(v);
   else
      return (double) HTS_mseq(v);
}

/* HTS_Vocoder_initialize_excitation: initialize excitation */
static void HTS_Vocoder_initialize_excitation(HTS_Vocoder * v, double pitch, size_t nlpf)
{
   size_t i;

   v->pitch_of_curr_point = pitch;
   v->pitch_counter = pitch;
   v->pitch_inc_per_point = 0.0;
   if (nlpf > 0) {
      v->excite_buff_size = nlpf;
      v->excite_ring_buff = (double *) HTS_calloc(v->excite_buff_size, sizeof(double));
      for (i = 0; i < v->excite_buff_size; i++)
         v->excite_ring_buff[i] = 0.0;
      v->excite_buff_index = 0;
   } else {
      v->excite_buff_size = 0;
      v->excite_ring_buff = NULL;
      v->excite_buff_index = 0;
   }
}

/* HTS_Vocoder_start_excitation: start excitation of each frame */
static void HTS_Vocoder_start_excitation(HTS_Vocoder * v, double pitch)
{
   if (v->pitch_of_curr_point != 0.0 && pitch != 0.0) {
      v->pitch_inc_per_point = (pitch - v->pitch_of_curr_point) / v->fprd;
   } else {
      v->pitch_inc_per_point = 0.0;
      v->pitch_of_curr_point = pitch;
      v->pitch_counter = pitch;
   }
}

/* HTS_Vocoder_excite_unvoiced_frame: ping noise to ring buffer */
static void HTS_Vocoder_excite_unvoiced_frame(HTS_Vocoder * v, double noise)
{
   size_t center = (v->excite_buff_size - 1) / 2;
   v->excite_ring_buff[(v->excite_buff_index + center) % v->excite_buff_size] += noise;
}

/* HTS_Vocoder_excite_vooiced_frame: ping noise and pulse to ring buffer */
static void HTS_Vocoder_excite_voiced_frame(HTS_Vocoder * v, double noise, double pulse, const double *lpf)
{
   size_t i;
   size_t center = (v->excite_buff_size - 1) / 2;

   if (noise != 0.0) {
      for (i = 0; i < v->excite_buff_size; i++) {
         if (i == center)
            v->excite_ring_buff[(v->excite_buff_index + i) % v->excite_buff_size] += noise * (1.0 - lpf[i]);
         else
            v->excite_ring_buff[(v->excite_buff_index + i) % v->excite_buff_size] += noise * (0.0 - lpf[i]);
      }
   }
   if (pulse != 0.0) {
      for (i = 0; i < v->excite_buff_size; i++)
         v->excite_ring_buff[(v->excite_buff_index + i) % v->excite_buff_size] += pulse * lpf[i];
   }
}

/* HTS_Vocoder_get_excitation: get excitation of each sample */
static double HTS_Vocoder_get_excitation(HTS_Vocoder * v, const double *lpf)
{
   double x;
   double noise, pulse = 0.0;

   if (v->excite_buff_size > 0) {
      noise = HTS_white_noise(v);
      pulse = 0.0;
      if (v->pitch_of_curr_point == 0.0) {
         HTS_Vocoder_excite_unvoiced_frame(v, noise);
      } else {
         v->pitch_counter += 1.0;
         if (v->pitch_counter >= v->pitch_of_curr_point) {
            pulse = sqrt(v->pitch_of_curr_point);
            v->pitch_counter -= v->pitch_of_curr_point;
         }
         HTS_Vocoder_excite_voiced_frame(v, noise, pulse, lpf);
         v->pitch_of_curr_point += v->pitch_inc_per_point;
      }
      x = v->excite_ring_buff[v->excite_buff_index];
      v->excite_ring_buff[v->excite_buff_index] = 0.0;
      v->excite_buff_index++;
      if (v->excite_buff_index >= v->excite_buff_size)
         v->excite_buff_index = 0;
   } else {
      if (v->pitch_of_curr_point == 0.0) {
         x = HTS_white_noise(v);
      } else {
         v->pitch_counter += 1.0;
         if (v->pitch_counter >= v->pitch_of_curr_point) {
            x = sqrt(v->pitch_of_curr_point);
            v->pitch_counter -= v->pitch_of_curr_point;
         } else {
            x = 0.0;
         }
         v->pitch_of_curr_point += v->pitch_inc_per_point;
      }
   }

   return x;
}

/* HTS_Vocoder_end_excitation: end excitation of each frame */
static void HTS_Vocoder_end_excitation(HTS_Vocoder * v, double pitch)
{
   v->pitch_of_curr_point = pitch;
}

/* HTS_Vocoder_postfilter_mcp: postfilter for MCP */
static void HTS_Vocoder_postfilter_mcp(HTS_Vocoder * v, double *mcp, const int m, double alpha, double beta)
{
   double e1, e2;
   int k;

   if (beta > 0.0 && m > 1) {
      if (v->postfilter_size < m) {
         if (v->postfilter_buff != NULL)
            HTS_free(v->postfilter_buff);
         v->postfilter_buff = (double *) HTS_calloc(m + 1, sizeof(double));
         v->postfilter_size = m;
      }
      HTS_mc2b(mcp, v->postfilter_buff, m, alpha);
      e1 = HTS_b2en(v, v->postfilter_buff, m, alpha);

      v->postfilter_buff[1] -= beta * alpha * v->postfilter_buff[2];
      for (k = 2; k <= m; k++)
         v->postfilter_buff[k] *= (1.0 + beta);

      e2 = HTS_b2en(v, v->postfilter_buff, m, alpha);
      v->postfilter_buff[0] += log(e1 / e2) / 2;
      HTS_b2mc(v->postfilter_buff, mcp, m, alpha);
   }
}

/* HTS_Vocoder_postfilter_lsp: postfilter for LSP */
static void HTS_Vocoder_postfilter_lsp(HTS_Vocoder * v, double *lsp, size_t m, double alpha, double beta)
{
   double e1, e2;
   size_t i;
   double d1, d2;

   if (beta > 0.0 && m > 1) {
      if (v->postfilter_size < m) {
         if (v->postfilter_buff != NULL)
            HTS_free(v->postfilter_buff);
         v->postfilter_buff = (double *) HTS_calloc(m + 1, sizeof(double));
         v->postfilter_size = m;
      }

      e1 = HTS_lsp2en(v, lsp, m, alpha);

      /* postfiltering */
      for (i = 0; i <= m; i++) {
         if (i > 1 && i < m) {
            d1 = beta * (lsp[i + 1] - lsp[i]);
            d2 = beta * (lsp[i] - lsp[i - 1]);
            v->postfilter_buff[i] = lsp[i - 1] + d2 + (d2 * d2 * ((lsp[i + 1] - lsp[i - 1]) - (d1 + d2))) / ((d2 * d2) + (d1 * d1));
         } else {
            v->postfilter_buff[i] = lsp[i];
         }
      }
      HTS_movem(v->postfilter_buff, lsp, m + 1);

      e2 = HTS_lsp2en(v, lsp, m, alpha);

      if (e1 != e2) {
         if (v->use_log_gain)
            lsp[0] += 0.5 * log(e1 / e2);
         else
            lsp[0] *= sqrt(e1 / e2);
      }
   }
}

/* HTS_Vocoder_initialize: initialize vocoder */
void HTS_Vocoder_initialize(HTS_Vocoder * v, size_t m, size_t stage, HTS_Boolean use_log_gain, size_t rate, size_t fperiod)
{
   /* set parameter */
   v->is_first = TRUE;
   v->stage = stage;
   if (stage != 0)
      v->gamma = -1.0 / v->stage;
   else
      v->gamma = 0.0;
   v->use_log_gain = use_log_gain;
   v->fprd = fperiod;
   v->next = SEED;
   v->gauss = GAUSS;
   v->rate = rate;
   v->pitch_of_curr_point = 0.0;
   v->pitch_counter = 0.0;
   v->pitch_inc_per_point = 0.0;
   v->excite_ring_buff = NULL;
   v->excite_buff_size = 0;
   v->excite_buff_index = 0;
   v->sw = 0;
   v->x = 0x55555555;
   /* init buffer */
   v->freqt_buff = NULL;
   v->freqt_size = 0;
   v->gc2gc_buff = NULL;
   v->gc2gc_size = 0;
   v->lsp2lpc_buff = NULL;
   v->lsp2lpc_size = 0;
   v->postfilter_buff = NULL;
   v->postfilter_size = 0;
   v->spectrum2en_buff = NULL;
   v->spectrum2en_size = 0;
   if (v->stage == 0) {         /* for MCP */
      v->c = (double *) HTS_calloc(m * (3 + PADEORDER) + 5 * PADEORDER + 6, sizeof(double));
      v->cc = v->c + m + 1;
      v->cinc = v->cc + m + 1;
      v->d1 = v->cinc + m + 1;
   } else {                     /* for LSP */
      v->c = (double *) HTS_calloc((m + 1) * (v->stage + 3), sizeof(double));
      v->cc = v->c + m + 1;
      v->cinc = v->cc + m + 1;
      v->d1 = v->cinc + m + 1;
   }
}

/* HTS_Vocoder_synthesize: pulse/noise excitation and MLSA/MGLSA filster based waveform synthesis */
void HTS_Vocoder_synthesize(HTS_Vocoder * v, size_t m, double lf0, double *spectrum, size_t nlpf, double *lpf, double alpha, double beta, double volume, double *rawdata, HTS_Audio * audio)
{
   double x;
   int i, j;
   short xs;
   int rawidx = 0;
   double p;

   /* lf0 -> pitch */
   if (lf0 == LZERO)
      p = 0.0;
   else if (lf0 <= MIN_LF0)
      p = v->rate / MIN_F0;
   else if (lf0 >= MAX_LF0)
      p = v->rate / MAX_F0;
   else
      p = v->rate / exp(lf0);

   /* first time */
   if (v->is_first == TRUE) {
      HTS_Vocoder_initialize_excitation(v, p, nlpf);
      if (v->stage == 0) {      /* for MCP */
         HTS_mc2b(spectrum, v->c, m, alpha);
      } else {                  /* for LSP */
         HTS_movem(spectrum, v->c, m + 1);
         HTS_lsp2mgc(v, v->c, v->c, m, alpha);
         HTS_mc2b(v->c, v->c, m, alpha);
         HTS_gnorm(v->c, v->c, m, v->gamma);
         for (i = 1; i <= m; i++)
            v->c[i] *= v->gamma;
      }
      v->is_first = FALSE;
   }

   HTS_Vocoder_start_excitation(v, p);
   if (v->stage == 0) {         /* for MCP */
      HTS_Vocoder_postfilter_mcp(v, spectrum, m, alpha, beta);
      HTS_mc2b(spectrum, v->cc, m, alpha);
      for (i = 0; i <= m; i++)
         v->cinc[i] = (v->cc[i] - v->c[i]) / v->fprd;
   } else {                     /* for LSP */
      HTS_Vocoder_postfilter_lsp(v, spectrum, m, alpha, beta);
      HTS_check_lsp_stability(spectrum, m);
      HTS_lsp2mgc(v, spectrum, v->cc, m, alpha);
      HTS_mc2b(v->cc, v->cc, m, alpha);
      HTS_gnorm(v->cc, v->cc, m, v->gamma);
      for (i = 1; i <= m; i++)
         v->cc[i] *= v->gamma;
      for (i = 0; i <= m; i++)
         v->cinc[i] = (v->cc[i] - v->c[i]) / v->fprd;
   }

   for (j = 0; j < v->fprd; j++) {
      x = HTS_Vocoder_get_excitation(v, lpf);
      if (v->stage == 0) {      /* for MCP */
         if (x != 0.0)
            x *= exp(v->c[0]);
         x = HTS_mlsadf(x, v->c, m, alpha, PADEORDER, v->d1);
      } else {                  /* for LSP */
         if (!NGAIN)
            x *= v->c[0];
         x = HTS_mglsadf(x, v->c, m, alpha, v->stage, v->d1);
      }
      x *= volume;

      /* output */
      if (rawdata)
         rawdata[rawidx++] = x;
      if (audio) {
         if (x > 32767.0)
            xs = 32767;
         else if (x < -32768.0)
            xs = -32768;
         else
            xs = (short) x;
         HTS_Audio_write(audio, xs);
      }

      for (i = 0; i <= m; i++)
         v->c[i] += v->cinc[i];
   }

   HTS_Vocoder_end_excitation(v, p);
   HTS_movem(v->cc, v->c, m + 1);
}

/* HTS_Vocoder_clear: clear vocoder */
void HTS_Vocoder_clear(HTS_Vocoder * v)
{
   if (v != NULL) {
      /* free buffer */
      if (v->freqt_buff != NULL) {
         HTS_free(v->freqt_buff);
         v->freqt_buff = NULL;
      }
      v->freqt_size = 0;
      if (v->gc2gc_buff != NULL) {
         HTS_free(v->gc2gc_buff);
         v->gc2gc_buff = NULL;
      }
      v->gc2gc_size = 0;
      if (v->lsp2lpc_buff != NULL) {
         HTS_free(v->lsp2lpc_buff);
         v->lsp2lpc_buff = NULL;
      }
      v->lsp2lpc_size = 0;
      if (v->postfilter_buff != NULL) {
         HTS_free(v->postfilter_buff);
         v->postfilter_buff = NULL;
      }
      v->postfilter_size = 0;
      if (v->spectrum2en_buff != NULL) {
         HTS_free(v->spectrum2en_buff);
         v->spectrum2en_buff = NULL;
      }
      v->spectrum2en_size = 0;
      if (v->c != NULL) {
         HTS_free(v->c);
         v->c = NULL;
      }
      v->excite_buff_size = 0;
      v->excite_buff_index = 0;
      if (v->excite_ring_buff != NULL) {
         HTS_free(v->excite_ring_buff);
         v->excite_ring_buff = NULL;
      }
   }
}

HTS_VOCODER_C_END;

#endif                          /* !HTS_VOCODER_C */