1 /* ----------------------------------------------------------------- */
2 /*           The HMM-Based Speech Synthesis Engine "hts_engine API"  */
3 /*           developed by HTS Working Group                          */
4 /*           http://hts-engine.sourceforge.net/                      */
5 /* ----------------------------------------------------------------- */
6 /*                                                                   */
7 /*  Copyright (c) 2001-2011  Nagoya Institute of Technology          */
8 /*                           Department of Computer Science          */
9 /*                                                                   */
10 /*                2001-2008  Tokyo Institute of Technology           */
11 /*                           Interdisciplinary Graduate School of    */
12 /*                           Science and Engineering                 */
13 /*                                                                   */
14 /* All rights reserved.                                              */
15 /*                                                                   */
16 /* Redistribution and use in source and binary forms, with or        */
17 /* without modification, are permitted provided that the following   */
18 /* conditions are met:                                               */
19 /*                                                                   */
20 /* - Redistributions of source code must retain the above copyright  */
21 /*   notice, this list of conditions and the following disclaimer.   */
22 /* - Redistributions in binary form must reproduce the above         */
23 /*   copyright notice, this list of conditions and the following     */
24 /*   disclaimer in the documentation and/or other materials provided */
25 /*   with the distribution.                                          */
26 /* - Neither the name of the HTS working group nor the names of its  */
27 /*   contributors may be used to endorse or promote products derived */
28 /*   from this software without specific prior written permission.   */
29 /*                                                                   */
30 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND            */
31 /* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,       */
32 /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF          */
33 /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE          */
34 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS */
35 /* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,          */
36 /* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED   */
37 /* TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,     */
38 /* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON */
39 /* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,   */
40 /* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    */
41 /* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           */
42 /* POSSIBILITY OF SUCH DAMAGE.                                       */
43 /* ----------------------------------------------------------------- */
44 
45 #ifndef HTS106_VOCODER_C
46 #define HTS106_VOCODER_C
47 
48 #ifdef __cplusplus
49 #define HTS106_VOCODER_C_START extern "C" {
50 #define HTS106_VOCODER_C_END   }
51 #else
52 #define HTS106_VOCODER_C_START
53 #define HTS106_VOCODER_C_END
54 #endif                          /* __CPLUSPLUS */
55 
56 HTS106_VOCODER_C_START;
57 
58 #include <math.h>               /* for sqrt(),log(),exp(),pow(),cos() */
59 
60 /* hts_engine libraries */
61 #include "HTS106_hidden.h"
62 
63 static const double HTS106_pade[21] = {
64    1.00000000000,
65    1.00000000000,
66    0.00000000000,
67    1.00000000000,
68    0.00000000000,
69    0.00000000000,
70    1.00000000000,
71    0.00000000000,
72    0.00000000000,
73    0.00000000000,
74    1.00000000000,
75    0.49992730000,
76    0.10670050000,
77    0.01170221000,
78    0.00056562790,
79    1.00000000000,
80    0.49993910000,
81    0.11070980000,
82    0.01369984000,
83    0.00095648530,
84    0.00003041721
85 };
86 
87 /* HTS106_movem: move memory */
HTS106_movem(double * a,double * b,const int nitem)88 static void HTS106_movem(double *a, double *b, const int nitem)
89 {
90    long i = (long) nitem;
91 
92    if (a > b)
93       while (i--)
94          *b++ = *a++;
95    else {
96       a += i;
97       b += i;
98       while (i--)
99          *--b = *--a;
100    }
101 }
102 
103 /* HTS106_mlsafir: sub functions for MLSA filter */
HTS106_mlsafir(const double x,const double * b,const int m,const double a,const double aa,double * d)104 static double HTS106_mlsafir(const double x, const double *b, const int m, const double a, const double aa, double *d)
105 {
106    double y = 0.0;
107    int i;
108 
109    d[0] = x;
110    d[1] = aa * d[0] + a * d[1];
111 
112    for (i = 2; i <= m; i++)
113       d[i] += a * (d[i + 1] - d[i - 1]);
114 
115    for (i = 2; i <= m; i++)
116       y += d[i] * b[i];
117 
118    for (i = m + 1; i > 1; i--)
119       d[i] = d[i - 1];
120 
121    return (y);
122 }
123 
124 /* HTS106_mlsadf1: sub functions for MLSA filter */
HTS106_mlsadf1(double x,const double * b,const int m,const double a,const double aa,const int pd,double * d,const double * ppade)125 static double HTS106_mlsadf1(double x, const double *b, const int m, const double a, const double aa, const int pd, double *d, const double *ppade)
126 {
127    double v, out = 0.0, *pt;
128    int i;
129 
130    pt = &d[pd + 1];
131 
132    for (i = pd; i >= 1; i--) {
133       d[i] = aa * pt[i - 1] + a * d[i];
134       pt[i] = d[i] * b[1];
135       v = pt[i] * ppade[i];
136       x += (1 & i) ? v : -v;
137       out += v;
138    }
139 
140    pt[0] = x;
141    out += x;
142 
143    return (out);
144 }
145 
146 /* HTS106_mlsadf2: sub functions for MLSA filter */
HTS106_mlsadf2(double x,const double * b,const int m,const double a,const double aa,const int pd,double * d,const double * ppade)147 static double HTS106_mlsadf2(double x, const double *b, const int m, const double a, const double aa, const int pd, double *d, const double *ppade)
148 {
149    double v, out = 0.0, *pt;
150    int i;
151 
152    pt = &d[pd * (m + 2)];
153 
154    for (i = pd; i >= 1; i--) {
155       pt[i] = HTS106_mlsafir(pt[i - 1], b, m, a, aa, &d[(i - 1) * (m + 2)]);
156       v = pt[i] * ppade[i];
157 
158       x += (1 & i) ? v : -v;
159       out += v;
160    }
161 
162    pt[0] = x;
163    out += x;
164 
165    return (out);
166 }
167 
168 /* HTS106_mlsadf: functions for MLSA filter */
HTS106_mlsadf(double x,const double * b,const int m,const double a,const int pd,double * d)169 static double HTS106_mlsadf(double x, const double *b, const int m, const double a, const int pd, double *d)
170 {
171    const double aa = 1 - a * a;
172    const double *ppade = &(HTS106_pade[pd * (pd + 1) / 2]);
173 
174    x = HTS106_mlsadf1(x, b, m, a, aa, pd, d, ppade);
175    x = HTS106_mlsadf2(x, b, m, a, aa, pd, &d[2 * (pd + 1)], ppade);
176 
177    return (x);
178 }
179 
180 /* HTS106_rnd: functions for random noise generation */
HTS106_rnd(unsigned long * next)181 static double HTS106_rnd(unsigned long *next)
182 {
183    double r;
184 
185    *next = *next * 1103515245L + 12345;
186    r = (*next / 65536L) % 32768L;
187 
188    return (r / RANDMAX);
189 }
190 
191 /* HTS106_nrandom: functions for gaussian random noise generation */
HTS106_nrandom(HTS106_Vocoder * v)192 static double HTS106_nrandom(HTS106_Vocoder * v)
193 {
194    if (v->sw == 0) {
195       v->sw = 1;
196       do {
197          v->r1 = 2 * HTS106_rnd(&v->next) - 1;
198          v->r2 = 2 * HTS106_rnd(&v->next) - 1;
199          v->s = v->r1 * v->r1 + v->r2 * v->r2;
200       } while (v->s > 1 || v->s == 0);
201       v->s = sqrt(-2 * log(v->s) / v->s);
202       return (v->r1 * v->s);
203    } else {
204       v->sw = 0;
205       return (v->r2 * v->s);
206    }
207 }
208 
209 /* HTS106_srnd: functions for gaussian random noise generation */
HTS106_srnd(unsigned long seed)210 static unsigned long HTS106_srnd(unsigned long seed)
211 {
212    return (seed);
213 }
214 
215 /* HTS106_mceq: function for M-sequence random noise generation */
HTS106_mseq(HTS106_Vocoder * v)216 static int HTS106_mseq(HTS106_Vocoder * v)
217 {
218    int x0, x28;
219 
220    v->x >>= 1;
221    if (v->x & B0)
222       x0 = 1;
223    else
224       x0 = -1;
225    if (v->x & B28)
226       x28 = 1;
227    else
228       x28 = -1;
229    if (x0 + x28)
230       v->x &= B31_;
231    else
232       v->x |= B31;
233 
234    return (x0);
235 }
236 
237 /* HTS106_mc2b: transform mel-cepstrum to MLSA digital fillter coefficients */
HTS106_mc2b(double * mc,double * b,int m,const double a)238 static void HTS106_mc2b(double *mc, double *b, int m, const double a)
239 {
240    if (mc != b) {
241       if (a != 0.0) {
242          b[m] = mc[m];
243          for (m--; m >= 0; m--)
244             b[m] = mc[m] - a * b[m + 1];
245       } else
246          HTS106_movem(mc, b, m + 1);
247    } else if (a != 0.0)
248       for (m--; m >= 0; m--)
249          b[m] -= a * b[m + 1];
250 }
251 
252 /* HTS106_b2bc: transform MLSA digital filter coefficients to mel-cepstrum */
HTS106_b2mc(const double * b,double * mc,int m,const double a)253 static void HTS106_b2mc(const double *b, double *mc, int m, const double a)
254 {
255    double d, o;
256 
257    d = mc[m] = b[m];
258    for (m--; m >= 0; m--) {
259       o = b[m] + a * d;
260       d = b[m];
261       mc[m] = o;
262    }
263 }
264 
265 /* HTS106_freqt: frequency transformation */
HTS106_freqt(HTS106_Vocoder * v,const double * c1,const int m1,double * c2,const int m2,const double a)266 static void HTS106_freqt(HTS106_Vocoder * v, const double *c1, const int m1, double *c2, const int m2, const double a)
267 {
268    int i, j;
269    const double b = 1 - a * a;
270    double *g;
271 
272    if (m2 > v->freqt_size) {
273       if (v->freqt_buff != NULL)
274          HTS106_free(v->freqt_buff);
275       v->freqt_buff = (double *) HTS106_calloc(m2 + m2 + 2, sizeof(double));
276       v->freqt_size = m2;
277    }
278    g = v->freqt_buff + v->freqt_size + 1;
279 
280    for (i = 0; i < m2 + 1; i++)
281       g[i] = 0.0;
282 
283    for (i = -m1; i <= 0; i++) {
284       if (0 <= m2)
285          g[0] = c1[-i] + a * (v->freqt_buff[0] = g[0]);
286       if (1 <= m2)
287          g[1] = b * v->freqt_buff[0] + a * (v->freqt_buff[1] = g[1]);
288       for (j = 2; j <= m2; j++)
289          g[j] = v->freqt_buff[j - 1] + a * ((v->freqt_buff[j] = g[j]) - g[j - 1]);
290    }
291 
292    HTS106_movem(g, c2, m2 + 1);
293 }
294 
295 /* HTS106_c2ir: The minimum phase impulse response is evaluated from the minimum phase cepstrum */
HTS106_c2ir(const double * c,const int nc,double * h,const int leng)296 static void HTS106_c2ir(const double *c, const int nc, double *h, const int leng)
297 {
298    int n, k, upl;
299    double d;
300 
301    h[0] = exp(c[0]);
302    for (n = 1; n < leng; n++) {
303       d = 0;
304       upl = (n >= nc) ? nc - 1 : n;
305       for (k = 1; k <= upl; k++)
306          d += k * c[k] * h[n - k];
307       h[n] = d / n;
308    }
309 }
310 
311 /* HTS106_b2en: calculate frame energy */
HTS106_b2en(HTS106_Vocoder * v,const double * b,const int m,const double a)312 static double HTS106_b2en(HTS106_Vocoder * v, const double *b, const int m, const double a)
313 {
314    int i;
315    double en = 0.0;
316    double *cep;
317    double *ir;
318 
319    if (v->spectrum2en_size < m) {
320       if (v->spectrum2en_buff != NULL)
321          HTS106_free(v->spectrum2en_buff);
322       v->spectrum2en_buff = (double *) HTS106_calloc((m + 1) + 2 * IRLENG, sizeof(double));
323       v->spectrum2en_size = m;
324    }
325    cep = v->spectrum2en_buff + m + 1;
326    ir = cep + IRLENG;
327 
328    HTS106_b2mc(b, v->spectrum2en_buff, m, a);
329    HTS106_freqt(v, v->spectrum2en_buff, m, cep, IRLENG - 1, -a);
330    HTS106_c2ir(cep, IRLENG, ir, IRLENG);
331 
332    for (i = 0; i < IRLENG; i++)
333       en += ir[i] * ir[i];
334 
335    return (en);
336 }
337 
338 /* HTS106_ignorm: inverse gain normalization */
HTS106_ignorm(double * c1,double * c2,int m,const double g)339 static void HTS106_ignorm(double *c1, double *c2, int m, const double g)
340 {
341    double k;
342    if (g != 0.0) {
343       k = pow(c1[0], g);
344       for (; m >= 1; m--)
345          c2[m] = k * c1[m];
346       c2[0] = (k - 1.0) / g;
347    } else {
348       HTS106_movem(&c1[1], &c2[1], m);
349       c2[0] = log(c1[0]);
350    }
351 }
352 
353 /* HTS106_gnorm: gain normalization */
HTS106_gnorm(double * c1,double * c2,int m,const double g)354 static void HTS106_gnorm(double *c1, double *c2, int m, const double g)
355 {
356    double k;
357    if (g != 0.0) {
358       k = 1.0 + g * c1[0];
359       for (; m >= 1; m--)
360          c2[m] = c1[m] / k;
361       c2[0] = pow(k, 1.0 / g);
362    } else {
363       HTS106_movem(&c1[1], &c2[1], m);
364       c2[0] = exp(c1[0]);
365    }
366 }
367 
368 /* HTS106_lsp2lpc: transform LSP to LPC */
HTS106_lsp2lpc(HTS106_Vocoder * v,double * lsp,double * a,const int m)369 static void HTS106_lsp2lpc(HTS106_Vocoder * v, double *lsp, double *a, const int m)
370 {
371    int i, k, mh1, mh2, flag_odd;
372    double xx, xf, xff;
373    double *p, *q;
374    double *a0, *a1, *a2, *b0, *b1, *b2;
375 
376    flag_odd = 0;
377    if (m % 2 == 0)
378       mh1 = mh2 = m / 2;
379    else {
380       mh1 = (m + 1) / 2;
381       mh2 = (m - 1) / 2;
382       flag_odd = 1;
383    }
384 
385    if (m > v->lsp2lpc_size) {
386       if (v->lsp2lpc_buff != NULL)
387          HTS106_free(v->lsp2lpc_buff);
388       v->lsp2lpc_buff = (double *) HTS106_calloc(5 * m + 6, sizeof(double));
389       v->lsp2lpc_size = m;
390    }
391    p = v->lsp2lpc_buff + m;
392    q = p + mh1;
393    a0 = q + mh2;
394    a1 = a0 + (mh1 + 1);
395    a2 = a1 + (mh1 + 1);
396    b0 = a2 + (mh1 + 1);
397    b1 = b0 + (mh2 + 1);
398    b2 = b1 + (mh2 + 1);
399 
400    HTS106_movem(lsp, v->lsp2lpc_buff, m);
401 
402    for (i = 0; i < mh1 + 1; i++)
403       a0[i] = 0.0;
404    for (i = 0; i < mh1 + 1; i++)
405       a1[i] = 0.0;
406    for (i = 0; i < mh1 + 1; i++)
407       a2[i] = 0.0;
408    for (i = 0; i < mh2 + 1; i++)
409       b0[i] = 0.0;
410    for (i = 0; i < mh2 + 1; i++)
411       b1[i] = 0.0;
412    for (i = 0; i < mh2 + 1; i++)
413       b2[i] = 0.0;
414 
415    /* lsp filter parameters */
416    for (i = k = 0; i < mh1; i++, k += 2)
417       p[i] = -2.0 * cos(v->lsp2lpc_buff[k]);
418    for (i = k = 0; i < mh2; i++, k += 2)
419       q[i] = -2.0 * cos(v->lsp2lpc_buff[k + 1]);
420 
421    /* impulse response of analysis filter */
422    xx = 1.0;
423    xf = xff = 0.0;
424 
425    for (k = 0; k <= m; k++) {
426       if (flag_odd) {
427          a0[0] = xx;
428          b0[0] = xx - xff;
429          xff = xf;
430          xf = xx;
431       } else {
432          a0[0] = xx + xf;
433          b0[0] = xx - xf;
434          xf = xx;
435       }
436 
437       for (i = 0; i < mh1; i++) {
438          a0[i + 1] = a0[i] + p[i] * a1[i] + a2[i];
439          a2[i] = a1[i];
440          a1[i] = a0[i];
441       }
442 
443       for (i = 0; i < mh2; i++) {
444          b0[i + 1] = b0[i] + q[i] * b1[i] + b2[i];
445          b2[i] = b1[i];
446          b1[i] = b0[i];
447       }
448 
449       if (k != 0)
450          a[k - 1] = -0.5 * (a0[mh1] + b0[mh2]);
451       xx = 0.0;
452    }
453 
454    for (i = m - 1; i >= 0; i--)
455       a[i + 1] = -a[i];
456    a[0] = 1.0;
457 }
458 
459 /* HTS106_gc2gc: generalized cepstral transformation */
HTS106_gc2gc(HTS106_Vocoder * v,double * c1,const int m1,const double g1,double * c2,const int m2,const double g2)460 static void HTS106_gc2gc(HTS106_Vocoder * v, double *c1, const int m1, const double g1, double *c2, const int m2, const double g2)
461 {
462    int i, min, k, mk;
463    double ss1, ss2, cc;
464 
465    if (m1 > v->gc2gc_size) {
466       if (v->gc2gc_buff != NULL)
467          HTS106_free(v->gc2gc_buff);
468       v->gc2gc_buff = (double *) HTS106_calloc(m1 + 1, sizeof(double));
469       v->gc2gc_size = m1;
470    }
471 
472    HTS106_movem(c1, v->gc2gc_buff, m1 + 1);
473 
474    c2[0] = v->gc2gc_buff[0];
475    for (i = 1; i <= m2; i++) {
476       ss1 = ss2 = 0.0;
477       min = m1 < i ? m1 : i - 1;
478       for (k = 1; k <= min; k++) {
479          mk = i - k;
480          cc = v->gc2gc_buff[k] * c2[mk];
481          ss2 += k * cc;
482          ss1 += mk * cc;
483       }
484 
485       if (i <= m1)
486          c2[i] = v->gc2gc_buff[i] + (g2 * ss2 - g1 * ss1) / i;
487       else
488          c2[i] = (g2 * ss2 - g1 * ss1) / i;
489    }
490 }
491 
492 /* HTS106_mgc2mgc: frequency and generalized cepstral transformation */
HTS106_mgc2mgc(HTS106_Vocoder * v,double * c1,const int m1,const double a1,const double g1,double * c2,const int m2,const double a2,const double g2)493 static void HTS106_mgc2mgc(HTS106_Vocoder * v, double *c1, const int m1, const double a1, const double g1, double *c2, const int m2, const double a2, const double g2)
494 {
495    double a;
496 
497    if (a1 == a2) {
498       HTS106_gnorm(c1, c1, m1, g1);
499       HTS106_gc2gc(v, c1, m1, g1, c2, m2, g2);
500       HTS106_ignorm(c2, c2, m2, g2);
501    } else {
502       a = (a2 - a1) / (1 - a1 * a2);
503       HTS106_freqt(v, c1, m1, c2, m2, a);
504       HTS106_gnorm(c2, c2, m2, g1);
505       HTS106_gc2gc(v, c2, m2, g1, c2, m2, g2);
506       HTS106_ignorm(c2, c2, m2, g2);
507    }
508 }
509 
510 /* HTS106_lsp2mgc: transform LSP to MGC */
HTS106_lsp2mgc(HTS106_Vocoder * v,double * lsp,double * mgc,const int m,const double alpha)511 static void HTS106_lsp2mgc(HTS106_Vocoder * v, double *lsp, double *mgc, const int m, const double alpha)
512 {
513    int i;
514    /* lsp2lpc */
515    HTS106_lsp2lpc(v, lsp + 1, mgc, m);
516    if (v->use_log_gain)
517       mgc[0] = exp(lsp[0]);
518    else
519       mgc[0] = lsp[0];
520 
521    /* mgc2mgc */
522    if (NORMFLG1)
523       HTS106_ignorm(mgc, mgc, m, v->gamma);
524    else if (MULGFLG1)
525       mgc[0] = (1.0 - mgc[0]) * v->stage;
526    if (MULGFLG1)
527       for (i = m; i >= 1; i--)
528          mgc[i] *= -v->stage;
529    HTS106_mgc2mgc(v, mgc, m, alpha, v->gamma, mgc, m, alpha, v->gamma);
530    if (NORMFLG2)
531       HTS106_gnorm(mgc, mgc, m, v->gamma);
532    else if (MULGFLG2)
533       mgc[0] = mgc[0] * v->gamma + 1.0;
534    if (MULGFLG2)
535       for (i = m; i >= 1; i--)
536          mgc[i] *= v->gamma;
537 }
538 
539 /* HTS106_mglsadff: sub functions for MGLSA filter */
HTS106_mglsadff(double x,const double * b,const int m,const double a,double * d)540 static double HTS106_mglsadff(double x, const double *b, const int m, const double a, double *d)
541 {
542    int i;
543 
544    double y;
545    y = d[0] * b[1];
546    for (i = 1; i < m; i++) {
547       d[i] += a * (d[i + 1] - d[i - 1]);
548       y += d[i] * b[i + 1];
549    }
550    x -= y;
551 
552    for (i = m; i > 0; i--)
553       d[i] = d[i - 1];
554    d[0] = a * d[0] + (1 - a * a) * x;
555    return x;
556 }
557 
558 /* HTS106_mglsadf: sub functions for MGLSA filter */
HTS106_mglsadf(double x,const double * b,const int m,const double a,const int n,double * d)559 static double HTS106_mglsadf(double x, const double *b, const int m, const double a, const int n, double *d)
560 {
561    int i;
562 
563    for (i = 0; i < n; i++)
564       x = HTS106_mglsadff(x, b, m, a, &d[i * (m + 1)]);
565 
566    return x;
567 }
568 
569 /* HTS106_white_noise: return white noise */
HTS106_white_noise(HTS106_Vocoder * v)570 static double HTS106_white_noise(HTS106_Vocoder * v)
571 {
572    if (v->gauss)
573       return (double) HTS106_nrandom(v);
574    else
575       return HTS106_mseq(v);
576 }
577 
578 /* HTS106_ping_pulse: ping pulse using low-pass filter */
HTS106_ping_pulse(HTS106_Vocoder * v,const int ping_place,const double p,const int nlpf,const double * lpf)579 static void HTS106_ping_pulse(HTS106_Vocoder * v, const int ping_place, const double p, const int nlpf, const double *lpf)
580 {
581    int i, j;
582    const double power = sqrt(p);
583 
584    for (i = ping_place - nlpf, j = 0; i <= ping_place + nlpf; i++, j++)
585       if (0 <= i && i < PULSELISTSIZE)
586          v->pulse_list[i] += power * lpf[j];
587 }
588 
589 /* HTS106_ping_noise: ping noise using low-pass filter */
HTS106_ping_noise(HTS106_Vocoder * v,const int ping_place,const int nlpf,const double * lpf)590 static void HTS106_ping_noise(HTS106_Vocoder * v, const int ping_place, const int nlpf, const double *lpf)
591 {
592    int i, j;
593    const double power = HTS106_white_noise(v);
594 
595    for (i = ping_place - nlpf, j = 0; i <= ping_place + nlpf; i++, j++)
596       if (0 <= i && i < PULSELISTSIZE) {
597          if (j == nlpf)
598             v->pulse_list[i] += power * (1.0 - lpf[j]);
599          else
600             v->pulse_list[i] += power * (0.0 - lpf[j]);
601       }
602 }
603 
604 /* HTS106_Vocoder_initialize_excitation: initialize excitation */
HTS106_Vocoder_initialize_excitation(HTS106_Vocoder * v)605 static void HTS106_Vocoder_initialize_excitation(HTS106_Vocoder * v)
606 {
607    int i;
608 
609    for (i = 0; i < PULSELISTSIZE; i++)
610       v->pulse_list[i] = 0.0;
611    v->p1 = 0.0;
612    v->pc = 0.0;
613    v->p = 0.0;
614    v->inc = 0.0;
615 }
616 
617 /* HTS106_Vocoder_start_excitation: start excitation of each frame */
HTS106_Vocoder_start_excitation(HTS106_Vocoder * v,const double pitch,const int nlpf)618 static void HTS106_Vocoder_start_excitation(HTS106_Vocoder * v, const double pitch, const int nlpf)
619 {
620    if (v->p1 != 0.0 && pitch != 0.0)
621       v->inc = (pitch - v->p1) * v->iprd / v->fprd;
622    else {
623       v->inc = 0.0;
624       if (nlpf <= 0) {
625          v->p1 = 0.0;
626          v->pc = pitch;
627       }
628    }
629    v->p = pitch;
630 }
631 
632 /* HTS106_Vocoder_get_excitation: get excitation of each sample */
HTS106_Vocoder_get_excitation(HTS106_Vocoder * v,const int fprd_index,const int iprd_index,const int nlpf,const double * lpf)633 static double HTS106_Vocoder_get_excitation(HTS106_Vocoder * v, const int fprd_index, const int iprd_index, const int nlpf, const double *lpf)
634 {
635    double x;
636    int i, j;
637 
638    if (nlpf > 0) {
639       if (fprd_index == 0) {
640          if (v->p1 == 0.0) {
641             for (i = 0; i < v->fprd; i++)
642                v->pulse_list[i] += HTS106_white_noise(v);
643             if (v->p != 0.0) {
644                HTS106_ping_pulse(v, v->fprd + nlpf, v->p, nlpf, lpf);
645                v->pc = v->fprd + nlpf;
646             }
647          } else if (v->p == 0.0) {
648             for (i = nlpf; i < v->fprd; i++)
649                v->pulse_list[i] += HTS106_white_noise(v);
650          } else {
651             for (i = 0, j = (v->iprd + 1) / 2; i < v->fprd; i++) {
652                if ((v->pc + v->p1) <= i) {
653                   HTS106_ping_pulse(v, i, v->p1, nlpf, lpf);
654                   v->pc += v->p1;
655                }
656                if (!--j) {
657                   v->p1 += v->inc;
658                   j = v->iprd;
659                }
660             }
661             for (i = v->fprd; i < v->fprd + nlpf; i++) {
662                if ((v->pc + v->p) <= i) {
663                   HTS106_ping_pulse(v, i, v->p, nlpf, lpf);
664                   v->pc += v->p;
665                }
666             }
667             for (i = nlpf; i < v->fprd + nlpf; i++)
668                HTS106_ping_noise(v, i, nlpf, lpf);
669          }
670       }
671       x = v->pulse_list[fprd_index];
672    } else {
673       if (v->p1 == 0.0)
674          x = HTS106_white_noise(v);
675       else {
676          if ((v->pc += 1.0) >= v->p1) {
677             x = sqrt(v->p1);
678             v->pc -= v->p1;
679          } else
680             x = 0.0;
681       }
682       if (iprd_index <= 1)
683          v->p1 += v->inc;
684    }
685    return x;
686 }
687 
688 /* HTS106_Vocoder_end_excitation: end excitation of each frame */
HTS106_Vocoder_end_excitation(HTS106_Vocoder * v,const int nlpf)689 static void HTS106_Vocoder_end_excitation(HTS106_Vocoder * v, const int nlpf)
690 {
691    int i;
692 
693    if (nlpf > 0) {
694       v->pc -= v->fprd;
695       for (i = 0; i < PULSELISTSIZE; i++)
696          if (i < PULSELISTSIZE - v->fprd)
697             v->pulse_list[i] = v->pulse_list[i + v->fprd];
698          else
699             v->pulse_list[i] = 0.0;
700    }
701    v->p1 = v->p;
702 }
703 
704 /* HTS106_Vocoder_initialize: initialize vocoder */
HTS106_Vocoder_initialize(HTS106_Vocoder * v,const int m,const int stage,HTS106_Boolean use_log_gain,const int rate,const int fperiod)705 void HTS106_Vocoder_initialize(HTS106_Vocoder * v, const int m, const int stage, HTS106_Boolean use_log_gain, const int rate, const int fperiod)
706 {
707    /* set parameter */
708    v->stage = stage;
709    if (stage != 0)
710       v->gamma = -1.0 / v->stage;
711    else
712       v->gamma = 0.0;
713    v->use_log_gain = use_log_gain;
714    v->fprd = fperiod;
715    v->iprd = IPERIOD;
716    v->seed = SEED;
717    v->next = SEED;
718    v->gauss = GAUSS;
719    v->rate = rate;
720    v->p1 = -1.0;
721    v->sw = 0;
722    v->x = 0x55555555;
723    /* init buffer */
724    v->freqt_buff = NULL;
725    v->freqt_size = 0;
726    v->gc2gc_buff = NULL;
727    v->gc2gc_size = 0;
728    v->lsp2lpc_buff = NULL;
729    v->lsp2lpc_size = 0;
730    v->postfilter_buff = NULL;
731    v->postfilter_size = 0;
732    v->spectrum2en_buff = NULL;
733    v->spectrum2en_size = 0;
734    if (v->stage == 0) {         /* for MCP */
735       v->c = (double *) HTS106_calloc(m * (3 + PADEORDER) + 5 * PADEORDER + 6, sizeof(double));
736       v->cc = v->c + m + 1;
737       v->cinc = v->cc + m + 1;
738       v->d1 = v->cinc + m + 1;
739    } else {                     /* for LSP */
740       v->c = (double *) HTS106_calloc((m + 1) * (v->stage + 3), sizeof(double));
741       v->cc = v->c + m + 1;
742       v->cinc = v->cc + m + 1;
743       v->d1 = v->cinc + m + 1;
744    }
745    v->pulse_list = (double *) HTS106_calloc(PULSELISTSIZE, sizeof(double));
746 }
747 
748 /* HTS106_Vocoder_synthesize: pulse/noise excitation and MLSA/MGLSA filster based waveform synthesis */
HTS106_Vocoder_synthesize(HTS106_Vocoder * v,const int m,double lf0,double * spectrum,const int nlpf,double * lpf,double alpha,double beta,double volume,short * rawdata,HTS106_Audio * audio)749 void HTS106_Vocoder_synthesize(HTS106_Vocoder * v, const int m, double lf0, double *spectrum, const int nlpf, double *lpf, double alpha, double beta, double volume, short *rawdata, HTS106_Audio * audio)
750 {
751    double x;
752    int i, j;
753    short xs;
754    int rawidx = 0;
755    double p;
756 
757    /* lf0 -> pitch */
758    if (lf0 == LZERO)
759       p = 0.0;
760    else
761       p = v->rate / exp(lf0);
762 
763    /* first time */
764    if (v->p1 < 0.0) {
765       if (v->gauss & (v->seed != 1))
766          v->next = HTS106_srnd((unsigned) v->seed);
767       HTS106_Vocoder_initialize_excitation(v);
768       if (v->stage == 0) {      /* for MCP */
769          HTS106_mc2b(spectrum, v->c, m, alpha);
770       } else {                  /* for LSP */
771          if (v->use_log_gain)
772             v->c[0] = LZERO;
773          else
774             v->c[0] = ZERO;
775          for (i = 1; i <= m; i++)
776             v->c[i] = i * PI / (m + 1);
777          HTS106_lsp2mgc(v, v->c, v->c, m, alpha);
778          HTS106_mc2b(v->c, v->c, m, alpha);
779          HTS106_gnorm(v->c, v->c, m, v->gamma);
780          for (i = 1; i <= m; i++)
781             v->c[i] *= v->gamma;
782       }
783    }
784 
785    HTS106_Vocoder_start_excitation(v, p, nlpf);
786    if (v->stage == 0) {         /* for MCP */
787       HTS106_Vocoder_postfilter_mcp(v, spectrum, m, alpha, beta);
788       HTS106_mc2b(spectrum, v->cc, m, alpha);
789       for (i = 0; i <= m; i++)
790          v->cinc[i] = (v->cc[i] - v->c[i]) * v->iprd / v->fprd;
791    } else {                     /* for LSP */
792       HTS106_lsp2mgc(v, spectrum, v->cc, m, alpha);
793       HTS106_mc2b(v->cc, v->cc, m, alpha);
794       HTS106_gnorm(v->cc, v->cc, m, v->gamma);
795       for (i = 1; i <= m; i++)
796          v->cc[i] *= v->gamma;
797       for (i = 0; i <= m; i++)
798          v->cinc[i] = (v->cc[i] - v->c[i]) * v->iprd / v->fprd;
799    }
800 
801    for (j = 0, i = (v->iprd + 1) / 2; j < v->fprd; j++) {
802       x = HTS106_Vocoder_get_excitation(v, j, i, nlpf, lpf);
803       if (v->stage == 0) {      /* for MCP */
804          if (x != 0.0)
805             x *= exp(v->c[0]);
806          x = HTS106_mlsadf(x, v->c, m, alpha, PADEORDER, v->d1);
807       } else {                  /* for LSP */
808          if (!NGAIN)
809             x *= v->c[0];
810          x = HTS106_mglsadf(x, v->c, m, alpha, v->stage, v->d1);
811       }
812       x *= volume;
813 
814       /* output */
815       if (x > 32767.0)
816          xs = 32767;
817       else if (x < -32768.0)
818          xs = -32768;
819       else
820          xs = (short) x;
821       if (rawdata)
822          rawdata[rawidx++] = xs;
823       if (audio)
824          HTS106_Audio_write(audio, xs);
825 
826       if (!--i) {
827          for (i = 0; i <= m; i++)
828             v->c[i] += v->cinc[i];
829          i = v->iprd;
830       }
831    }
832 
833    HTS106_Vocoder_end_excitation(v, nlpf);
834    HTS106_movem(v->cc, v->c, m + 1);
835 }
836 
837 /* HTS106_Vocoder_postfilter_mcp: postfilter for MCP */
HTS106_Vocoder_postfilter_mcp(HTS106_Vocoder * v,double * mcp,const int m,double alpha,double beta)838 void HTS106_Vocoder_postfilter_mcp(HTS106_Vocoder * v, double *mcp, const int m, double alpha, double beta)
839 {
840    double e1, e2;
841    int k;
842 
843    if (beta > 0.0 && m > 1) {
844       if (v->postfilter_size < m) {
845          if (v->postfilter_buff != NULL)
846             HTS106_free(v->postfilter_buff);
847          v->postfilter_buff = (double *) HTS106_calloc(m + 1, sizeof(double));
848          v->postfilter_size = m;
849       }
850       HTS106_mc2b(mcp, v->postfilter_buff, m, alpha);
851       e1 = HTS106_b2en(v, v->postfilter_buff, m, alpha);
852 
853       v->postfilter_buff[1] -= beta * alpha * mcp[2];
854       for (k = 2; k <= m; k++)
855          v->postfilter_buff[k] *= (1.0 + beta);
856 
857       e2 = HTS106_b2en(v, v->postfilter_buff, m, alpha);
858       v->postfilter_buff[0] += log(e1 / e2) / 2;
859       HTS106_b2mc(v->postfilter_buff, mcp, m, alpha);
860    }
861 }
862 
863 /* HTS106_Vocoder_clear: clear vocoder */
HTS106_Vocoder_clear(HTS106_Vocoder * v)864 void HTS106_Vocoder_clear(HTS106_Vocoder * v)
865 {
866    if (v != NULL) {
867       /* free buffer */
868       if (v->freqt_buff != NULL) {
869          HTS106_free(v->freqt_buff);
870          v->freqt_buff = NULL;
871       }
872       v->freqt_size = 0;
873       if (v->gc2gc_buff != NULL) {
874          HTS106_free(v->gc2gc_buff);
875          v->gc2gc_buff = NULL;
876       }
877       v->gc2gc_size = 0;
878       if (v->lsp2lpc_buff != NULL) {
879          HTS106_free(v->lsp2lpc_buff);
880          v->lsp2lpc_buff = NULL;
881       }
882       v->lsp2lpc_size = 0;
883       if (v->postfilter_buff != NULL) {
884          HTS106_free(v->postfilter_buff);
885          v->postfilter_buff = NULL;
886       }
887       v->postfilter_size = 0;
888       if (v->spectrum2en_buff != NULL) {
889          HTS106_free(v->spectrum2en_buff);
890          v->spectrum2en_buff = NULL;
891       }
892       v->spectrum2en_size = 0;
893       if (v->c != NULL) {
894          HTS106_free(v->c);
895          v->c = NULL;
896       }
897       if (v->pulse_list != NULL)
898          HTS106_free(v->pulse_list);
899    }
900 }
901 
902 HTS106_VOCODER_C_END;
903 
904 #endif                          /* !HTS106_VOCODER_C */
905