Bug Summary

File:modules/clustergen/mlsa_resynthesis.cc
Location:line 530, column 4
Description:Value stored to 'aa' is never read

Annotated Source Code

1/* --------------------------------------------------------------- */
2/* The HMM-Based Speech Synthesis System (HTS): version 1.1b */
3/* HTS Working Group */
4/* */
5/* Department of Computer Science */
6/* Nagoya Institute of Technology */
7/* and */
8/* Interdisciplinary Graduate School of Science and Engineering */
9/* Tokyo Institute of Technology */
10/* Copyright (c) 2001-2003 */
11/* All Rights Reserved. */
12/* */
13/* Permission is hereby granted, free of charge, to use and */
14/* distribute this software and its documentation without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of this work, and to permit persons to whom this */
18/* work is furnished to do so, subject to the following conditions: */
19/* */
20/* 1. The code must retain the above copyright notice, this list */
21/* of conditions and the following disclaimer. */
22/* */
23/* 2. Any modifications must be clearly marked as such. */
24/* */
25/* NAGOYA INSTITUTE OF TECHNOLOGY, TOKYO INSITITUTE OF TECHNOLOGY, */
26/* HTS WORKING GROUP, AND THE CONTRIBUTORS TO THIS WORK DISCLAIM */
27/* ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL */
28/* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
29/* SHALL NAGOYA INSTITUTE OF TECHNOLOGY, TOKYO INSITITUTE OF */
30/* TECHNOLOGY, HTS WORKING GROUP, NOR THE CONTRIBUTORS BE LIABLE */
31/* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY */
32/* DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, */
33/* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTUOUS */
34/* ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR */
35/* PERFORMANCE OF THIS SOFTWARE. */
36/* */
37/* --------------------------------------------------------------- */
38/* This is Zen's MLSA filter as ported by Toda to festvox vc */
39/* and back ported into hts/festival so we can do MLSA filtering */
40/* If I took more time I could probably make this use the same as */
41/* as the other code in this directory -- awb@cs.cmu.edu 03JAN06 */
42/* --------------------------------------------------------------- */
43
44/*********************************************************************/
45/* */
46/* Mel-cepstral vocoder (pulse/noise excitation & MLSA filter) */
47/* 2003/12/26 by Heiga Zen */
48/* */
49/* Extracted from HTS and slightly modified */
50/* by Tomoki Toda (tomoki@ics.nitech.ac.jp) */
51/* June 2004 */
52/* Integrate as a Voice Conversion module */
53/* */
54/*-------------------------------------------------------------------*/
55
56#include <stdio.h>
57#include <stdlib.h>
58#include <string.h>
59#include <math.h>
60#include <EST_walloc.h>
61#include "festival.h"
62
63#include "mlsa_resynthesis.h"
64
65static void wavecompressor(DVECTOR wav);
66static DVECTOR xdvalloc(long length);
67static DVECTOR xdvcut(DVECTOR x, long offset, long length);
68static void xdvfree(DVECTOR vector);
69static double dvmax(DVECTOR x, long *index);
70static double dvmin(DVECTOR x, long *index);
71static DMATRIX xdmalloc(long row, long col);
72static void xdmfree(DMATRIX matrix);
73
74static void waveampcheck(DVECTOR wav, XBOOLint msg_flag);
75
76static void init_vocoder(double fs, int framel, int m, VocoderSetup *vs);
77static void vocoder(double p, double *mc, EST_Track *str,
78 int t,
79 int m, double a, double beta,
80 VocoderSetup *vs, double *wav, long *pos);
81static double mlsadf(double x, double *b, int m, double a, int pd, double *d,
82 VocoderSetup *vs);
83static double mlsadf1(double x, double *b, int m, double a, int pd, double *d,
84 VocoderSetup *vs);
85static double mlsadf2(double x, double *b, int m, double a, int pd, double *d,
86 VocoderSetup *vs);
87static double mlsafir (double x, double *b, int m, double a, double *d);
88static double nrandom (VocoderSetup *vs);
89static double rnd (unsigned long *next);
90static unsigned long srnd (unsigned long seed);
91static int mseq (VocoderSetup *vs);
92static void mc2b (double *mc, double *b, int m, double a);
93static double b2en (double *b, int m, double a, VocoderSetup *vs);
94static void b2mc (double *b, double *mc, int m, double a);
95static void freqt (double *c1, int m1, double *c2, int m2, double a,
96 VocoderSetup *vs);
97static void c2ir (double *c, int nc, double *h, int leng);
98
99
100#if 0
101static DVECTOR get_dpowvec(DMATRIX rmcep, DMATRIX cmcep);
102static double get_dpow(double *rmcep, double *cmcep, int m, double a,
103 VocoderSetup *vs);
104#endif
105static void free_vocoder(VocoderSetup *vs);
106
107LISP mlsa_resynthesis(LISP ltrack, LISP strtrack)
108{
109 /* Resynthesizes a wave from given track */
110 EST_Track *t;
111 EST_Track *str = 0;
112 EST_Wave *wave = 0;
113 DVECTOR w;
114 DMATRIX mcep;
115 DVECTOR f0v;
116 int sr = 16000;
117 int i,j;
118 double shift;
119 double ALPHA = 0.42;
120 double BETA = 0.0;
121
122 if ((ltrack == NULL__null) ||
123 (TYPEP(ltrack,tc_string)( (ltrack != __null) && ((((ltrack) == ((struct obj *
) 0)) ? 0 : ((*(ltrack)).type)) == (13)) ) &&
124 (streq(get_c_string(ltrack),"nil")(strcmp(get_c_string(ltrack),"nil")==0))))
125 return siod(new EST_Wave(0,1,sr));
126
127 t = track(ltrack);
128
129 if (strtrack != NULL__null)
130 { /* We have to do mixed-excitation */
131 str = track(strtrack);
132 }
133
134 f0v = xdvalloc(t->num_frames());
135 mcep = xdmalloc(t->num_frames(),t->num_channels()-1);
136
137 for (i=0; i<t->num_frames(); i++)
138 {
139 f0v->data[i] = t->a(i,0);
140 for (j=1; j<t->num_channels(); j++)
141 mcep->data[i][j-1] = t->a(i,j);
142 }
143
144 if (t->num_frames() > 1)
145 shift = 1000.0*(t->t(1)-t->t(0));
146 else
147 shift = 5.0;
148
149 ALPHA = FLONM(siod_get_lval("mlsa_alpha_param",((*siod_get_lval("mlsa_alpha_param", "mlsa: mlsa_alpha_param not set"
)).storage_as.flonum.data)
150 "mlsa: mlsa_alpha_param not set"))((*siod_get_lval("mlsa_alpha_param", "mlsa: mlsa_alpha_param not set"
)).storage_as.flonum.data);
151 BETA = FLONM(siod_get_lval("mlsa_beta_param",((*siod_get_lval("mlsa_beta_param", "mlsa: mlsa_beta_param not set"
)).storage_as.flonum.data)
152 "mlsa: mlsa_beta_param not set"))((*siod_get_lval("mlsa_beta_param", "mlsa: mlsa_beta_param not set"
)).storage_as.flonum.data);
153
154 w = synthesis_body(mcep,f0v,str,sr,shift,ALPHA,BETA);
155
156 wave = new EST_Wave(w->length,1,sr);
157
158 for (i=0; i<w->length; i++)
159 wave->a(i) = (short)w->data[i];
160
161 xdmfree(mcep);
162 xdvfree(f0v);
163 xdvfree(w);
164
165 return siod(wave);
166}
167
168
169DVECTOR synthesis_body(DMATRIX mcep, // input mel-cep sequence
170 DVECTOR f0v, // input F0 sequence
171 EST_Track *str, // str for mixed excitation
172 double fs, // sampling frequency (Hz)
173 double framem, // FFT length
174 double alpha,
175 double beta)
176{
177 long t, pos;
178 int framel;
179 double f0;
180 VocoderSetup vs;
181 DVECTOR xd = NODATA__null;
182 DVECTOR syn = NODATA__null;
183 int i,j;
184
185 framel = (int)(framem * fs / 1000.0);
186 init_vocoder(fs, framel, mcep->col - 1, &vs);
187
188 if (str != NULL__null)
189 {
190 /* Mixed excitation filters */
191 LISP filters = siod_get_lval("me_mix_filters",
192 "mlsa: me_mix_filters not set");
193 LISP f;
194 int fl;
195 for (fl=0,f=filters; f; fl++)
196 f=cdr(f);
197 for (fl=0,f=filters; f; fl++)
198 f=cdr(f);
199 vs.ME_num = 5;
200 vs.ME_order = fl/vs.ME_num;
201
202 for (i=0; i < vs.ME_num; i++)
203 {
204 for (j=0; j<vs.ME_order; j++)
205 {
206 vs.h[i][j] = FLONM(car(filters))((*car(filters)).storage_as.flonum.data);
207 filters = cdr(filters);
208 }
209 }
210 vs.gauss = MFALSE;
211 }
212
213 // synthesize waveforms by MLSA filter
214 xd = xdvalloc(mcep->row * (framel + 2));
215 for (t = 0, pos = 0; t < mcep->row; t++) {
216 if (t >= f0v->length) f0 = 0.0;
217 else f0 = f0v->data[t];
218
219 vocoder(f0, mcep->data[t],
220 str, t,
221 mcep->col - 1,
222 alpha, beta, &vs,
223 xd->data, &pos);
224
225 }
226 syn = xdvcut(xd, 0, pos);
227
228 // normalized amplitude
229 /* waveampcheck(syn, XFALSE); */
230 wavecompressor(syn);
231
232 // memory free
233 xdvfree(xd);
234 free_vocoder(&vs);
235
236 return syn;
237}
238
239static void wavecompressor(DVECTOR wav)
240{
241 /* a somewhat over specific compressor */
242 int i;
243 double maxvalue, absv, d;
244 int sign;
245
246 maxvalue = MAX(FABS(dvmax(wav, NULL)), FABS(dvmin(wav, NULL)))((((dvmax(wav, __null)) >= 0.0 ? (dvmax(wav, __null)) : -(
dvmax(wav, __null)))) > (((dvmin(wav, __null)) >= 0.0 ?
(dvmin(wav, __null)) : -(dvmin(wav, __null)))) ? (((dvmax(wav
, __null)) >= 0.0 ? (dvmax(wav, __null)) : -(dvmax(wav, __null
)))) : (((dvmin(wav, __null)) >= 0.0 ? (dvmin(wav, __null)
) : -(dvmin(wav, __null)))));
247 for (i=0; i < wav->length; i++)
248 {
249 sign = ( wav->data[i] < 0 ) ? -1 : 1;
250 absv = FABS(wav->data[i])((wav->data[i]) >= 0.0 ? (wav->data[i]) : -(wav->
data[i]));
251 if (absv > 15000)
252 {
253 d = absv - 15000;
254 d /= maxvalue-15000;
255 wav->data[i] = sign*(12500+(5000*d));
256 }
257 else if (absv > 10000)
258 wav->data[i] = sign*(10000+((absv-10000)/2.0));
259 }
260
261}
262
263static void waveampcheck(DVECTOR wav, XBOOLint msg_flag)
264{
265 double value;
266 int k;
267
268 value = MAX(FABS(dvmax(wav, NULL)), FABS(dvmin(wav, NULL)))((((dvmax(wav, __null)) >= 0.0 ? (dvmax(wav, __null)) : -(
dvmax(wav, __null)))) > (((dvmin(wav, __null)) >= 0.0 ?
(dvmin(wav, __null)) : -(dvmin(wav, __null)))) ? (((dvmax(wav
, __null)) >= 0.0 ? (dvmax(wav, __null)) : -(dvmax(wav, __null
)))) : (((dvmin(wav, __null)) >= 0.0 ? (dvmin(wav, __null)
) : -(dvmin(wav, __null)))));
269 if (value >= 32000.0) {
270 if (msg_flag == XTRUE1) {
271 fprintf(stderrstderr, "amplitude is too big: %f\n", value);
272 fprintf(stderrstderr, "execute normalization\n");
273 }
274 /* was dvscoper(wav, "*", 32000.0 / value); */
275 for (k = 0; k < wav->length; k++) {
276 wav->data[k] = wav->data[k] * (32000.0/value);
277 if (wav->imag != NULL__null) {
278 wav->imag[k] = wav->imag[k] * (32000.0/value);
279 }
280 }
281 }
282
283 return;
284}
285
286static void init_vocoder(double fs, int framel, int m, VocoderSetup *vs)
287{
288 // initialize global parameter
289 int i;
290
291 vs->fprd = framel;
292 vs->iprd = 1;
293 vs->seed = 1;
294 vs->pd = 5;
295
296 vs->next =1;
297 vs->gauss = MTRUE;
298
299 vs->pade[ 0]=1.0;
300 vs->pade[ 1]=1.0; vs->pade[ 2]=0.0;
301 vs->pade[ 3]=1.0; vs->pade[ 4]=0.0; vs->pade[ 5]=0.0;
302 vs->pade[ 6]=1.0; vs->pade[ 7]=0.0; vs->pade[ 8]=0.0; vs->pade[ 9]=0.0;
303 vs->pade[10]=1.0; vs->pade[11]=0.4999273; vs->pade[12]=0.1067005; vs->pade[13]=0.01170221; vs->pade[14]=0.0005656279;
304 vs->pade[15]=1.0; vs->pade[16]=0.4999391; vs->pade[17]=0.1107098; vs->pade[18]=0.01369984; vs->pade[19]=0.0009564853;
305 vs->pade[20]=0.00003041721;
306
307 vs->rate = fs;
308 vs->c = wcalloc(double,3 * (m + 1) + 3 * (vs->pd + 1) + vs->pd * (m + 2))((double *)safe_wcalloc(sizeof(double)*(3 * (m + 1) + 3 * (vs
->pd + 1) + vs->pd * (m + 2))));
309
310 vs->p1 = -1;
311 vs->sw = 0;
312 vs->x = 0x55555555;
313
314 // for postfiltering
315 vs->mc = NULL__null;
316 vs->o = 0;
317 vs->d = NULL__null;
318 vs->irleng= 64;
319
320 // for MIXED EXCITATION
321 vs->ME_order = 48;
322 vs->ME_num = 5;
323 vs->hpulse = walloc(double,vs->ME_order)((double *)safe_walloc(sizeof(double)*(vs->ME_order)));
324 vs->hnoise = walloc(double,vs->ME_order)((double *)safe_walloc(sizeof(double)*(vs->ME_order)));
325 vs->xpulsesig = walloc(double,vs->ME_order)((double *)safe_walloc(sizeof(double)*(vs->ME_order)));
326 vs->xnoisesig = walloc(double,vs->ME_order)((double *)safe_walloc(sizeof(double)*(vs->ME_order)));
327 vs->h = walloc(double *,vs->ME_num)((double * *)safe_walloc(sizeof(double *)*(vs->ME_num)));
328 for (i=0; i< vs->ME_num; i++)
329 vs->h[i] = walloc(double,vs->ME_order)((double *)safe_walloc(sizeof(double)*(vs->ME_order)));
330
331 return;
332}
333
334static double plus_or_minus_one()
335{
336 /* Randomly return 1 or -1 */
337 if (rand() > RAND_MAX2147483647/2.0)
338 return 1.0;
339 else
340 return -1.0;
341}
342
343static void vocoder(double p, double *mc,
344 EST_Track *str, int t,
345 int m, double a, double beta,
346 VocoderSetup *vs, double *wav, long *pos)
347{
348 double inc, x, e1, e2;
349 int i, j, k;
350 double xpulse, xnoise;
351 double fxpulse, fxnoise;
352
353 if (str != NULL__null) /* MIXED-EXCITATION */
354 {
355 /* Copy in str's and build hpulse and hnoise for this frame */
356 for (i=0; i<vs->ME_order; i++)
357 {
358 vs->hpulse[i] = vs->hnoise[i] = 0.0;
359 for (j=0; j<vs->ME_num; j++)
360 {
361 vs->hpulse[i] += str->a(t,j) * vs->h[j][i];
362 vs->hnoise[i] += (1 - str->a(t,j)) * vs->h[j][i];
363 }
364 }
365 printf("awb_debug str %f %f %f %f %f\n",
366 str->a(t,0),
367 str->a(t,1),
368 str->a(t,2),
369 str->a(t,3),
370 str->a(t,4));
371 }
372
373 if (p != 0.0)
374 p = vs->rate / p; // f0 -> pitch
375
376 if (vs->p1 < 0) {
377 if (vs->gauss & (vs->seed != 1))
378 vs->next = srnd((unsigned)vs->seed);
379
380 vs->p1 = p;
381 vs->pc = vs->p1;
382 vs->cc = vs->c + m + 1;
383 vs->cinc = vs->cc + m + 1;
384 vs->d1 = vs->cinc + m + 1;
385
386 mc2b(mc, vs->c, m, a);
387
388 if (beta > 0.0 && m > 1) {
389 e1 = b2en(vs->c, m, a, vs);
390 vs->c[1] -= beta * a * mc[2];
391 for (k=2;k<=m;k++)
392 vs->c[k] *= (1.0 + beta);
393 e2 = b2en(vs->c, m, a, vs);
394 vs->c[0] += log(e1/e2)/2;
395 }
396
397 return;
398 }
399
400 mc2b(mc, vs->cc, m, a);
401 if (beta>0.0 && m > 1) {
402 e1 = b2en(vs->cc, m, a, vs);
403 vs->cc[1] -= beta * a * mc[2];
404 for (k = 2; k <= m; k++)
405 vs->cc[k] *= (1.0 + beta);
406 e2 = b2en(vs->cc, m, a, vs);
407 vs->cc[0] += log(e1 / e2) / 2.0;
408 }
409
410 for (k=0; k<=m; k++)
411 vs->cinc[k] = (vs->cc[k] - vs->c[k]) *
412 (double)vs->iprd / (double)vs->fprd;
413
414 if (vs->p1!=0.0 && p!=0.0) {
415 inc = (p - vs->p1) * (double)vs->iprd / (double)vs->fprd;
416 } else {
417 inc = 0.0;
418 vs->pc = p;
419 vs->p1 = 0.0;
420 }
421
422 for (j = vs->fprd, i = (vs->iprd + 1) / 2; j--;) {
423 if (vs->p1 == 0.0) {
424 if (vs->gauss)
425 x = (double) nrandom(vs);
426 else
427 x = plus_or_minus_one();
428
429 if (str != NULL__null) /* MIXED EXCITATION */
430 {
431 xnoise = x;
432 xpulse = 0.0;
433 }
434 } else {
435 if ((vs->pc += 1.0) >= vs->p1)
436 {
437 x = sqrt (vs->p1);
438 vs->pc = vs->pc - vs->p1;
439 }
440 else
441 x = 0.0;
442
443 if (str != NULL__null) /* MIXED EXCITATION */
444 {
445 xpulse = x;
446 xnoise = plus_or_minus_one();
447 }
448 }
449
450 /* MIXED EXCITATION */
451 /* The real work -- apply shaping filters to pulse and noise */
452 if (str != NULL__null)
453 {
454 fxpulse = fxnoise = 0.0;
455 for (k=vs->ME_order-1; k>0; k--)
456 {
457 fxpulse += vs->hpulse[k] * vs->xpulsesig[k];
458 fxnoise += vs->hnoise[k] * vs->xnoisesig[k];
459
460 vs->xpulsesig[k] = vs->xpulsesig[k-1];
461 vs->xnoisesig[k] = vs->xnoisesig[k-1];
462 }
463 fxpulse += vs->hpulse[0] * xpulse;
464 fxnoise += vs->hnoise[0] * xnoise;
465 vs->xpulsesig[0] = xpulse;
466 vs->xnoisesig[0] = xnoise;
467
468 x = fxpulse + fxnoise; /* excitation is pulse plus noise */
469 printf("awb_debug %f\n",(float)x);
470 }
471
472 x *= exp(vs->c[0]);
473
474 x = mlsadf(x, vs->c, m, a, vs->pd, vs->d1, vs);
475
476 wav[*pos] = x;
477 *pos += 1;
478
479 if (!--i) {
480 vs->p1 += inc;
481 for (k = 0; k <= m; k++) vs->c[k] += vs->cinc[k];
482 i = vs->iprd;
483 }
484 }
485
486 vs->p1 = p;
487 memmove(vs->c,vs->cc,sizeof(double)*(m+1));
488
489 return;
490}
491
492static double mlsadf(double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs)
493{
494
495 vs->ppade = &(vs->pade[pd*(pd+1)/2]);
496
497 x = mlsadf1 (x, b, m, a, pd, d, vs);
498 x = mlsadf2 (x, b, m, a, pd, &d[2*(pd+1)], vs);
499
500 return(x);
501}
502
503static double mlsadf1(double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs)
504{
505 double v, out = 0.0, *pt, aa;
506 register int i;
507
508 aa = 1 - a*a;
509 pt = &d[pd+1];
510
511 for (i=pd; i>=1; i--) {
512 d[i] = aa*pt[i-1] + a*d[i];
513 pt[i] = d[i] * b[1];
514 v = pt[i] * vs->ppade[i];
515 x += (1 & i) ? v : -v;
516 out += v;
517 }
518
519 pt[0] = x;
520 out += x;
521
522 return(out);
523}
524
525static double mlsadf2 (double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs)
526{
527 double v, out = 0.0, *pt, aa;
528 register int i;
529
530 aa = 1 - a*a;
Value stored to 'aa' is never read
531 pt = &d[pd * (m+2)];
532
533 for (i=pd; i>=1; i--) {
534 pt[i] = mlsafir (pt[i-1], b, m, a, &d[(i-1)*(m+2)]);
535 v = pt[i] * vs->ppade[i];
536
537 x += (1&i) ? v : -v;
538 out += v;
539 }
540
541 pt[0] = x;
542 out += x;
543
544 return(out);
545}
546
547static double mlsafir (double x, double *b, int m, double a, double *d)
548{
549 double y = 0.0;
550 double aa;
551 register int i;
552
553 aa = 1 - a*a;
554
555 d[0] = x;
556 d[1] = aa*d[0] + a*d[1];
557
558 for (i=2; i<=m; i++) {
559 d[i] = d[i] + a*(d[i+1]-d[i-1]);
560 y += d[i]*b[i];
561 }
562
563 for (i=m+1; i>1; i--)
564 d[i] = d[i-1];
565
566 return(y);
567}
568
569static double nrandom (VocoderSetup *vs)
570{
571 if (vs->sw == 0) {
572 vs->sw = 1;
573 do {
574 vs->r1 = 2.0 * rnd(&vs->next) - 1.0;
575 vs->r2 = 2.0 * rnd(&vs->next) - 1.0;
576 vs->s = vs->r1 * vs->r1 + vs->r2 * vs->r2;
577 } while (vs->s > 1 || vs->s == 0);
578
579 vs->s = sqrt (-2 * log(vs->s) / vs->s);
580
581 return(vs->r1*vs->s);
582 }
583 else {
584 vs->sw = 0;
585
586 return (vs->r2*vs->s);
587 }
588}
589
590static double rnd (unsigned long *next)
591{
592 double r;
593
594 *next = *next * 1103515245L + 12345;
595 r = (*next / 65536L) % 32768L;
596
597 return(r/RANDMAX32767);
598}
599
600static unsigned long srnd ( unsigned long seed )
601{
602 return(seed);
603}
604
605static int mseq (VocoderSetup *vs)
606{
607 register int x0, x28;
608
609 vs->x >>= 1;
610
611 if (vs->x & B00x00000001)
612 x0 = 1;
613 else
614 x0 = -1;
615
616 if (vs->x & B280x10000000)
617 x28 = 1;
618 else
619 x28 = -1;
620
621 if (x0 + x28)
622 vs->x &= B31_0x7fffffff;
623 else
624 vs->x |= B310x80000000;
625
626 return(x0);
627}
628
629// mc2b : transform mel-cepstrum to MLSA digital fillter coefficients
630static void mc2b (double *mc, double *b, int m, double a)
631{
632 b[m] = mc[m];
633
634 for (m--; m>=0; m--)
635 b[m] = mc[m] - a * b[m+1];
636
637 return;
638}
639
640
641static double b2en (double *b, int m, double a, VocoderSetup *vs)
642{
643 double en;
644 int k;
645
646 if (vs->o<m) {
647 if (vs->mc != NULL__null)
648 wfree(vs->mc);
649
650 vs->mc = wcalloc(double,(m + 1) + 2 * vs->irleng)((double *)safe_wcalloc(sizeof(double)*((m + 1) + 2 * vs->
irleng)));
651 vs->cep = vs->mc + m+1;
652 vs->ir = vs->cep + vs->irleng;
653 }
654
655 b2mc(b, vs->mc, m, a);
656 freqt(vs->mc, m, vs->cep, vs->irleng-1, -a, vs);
657 c2ir(vs->cep, vs->irleng, vs->ir, vs->irleng);
658 en = 0.0;
659
660 for (k=0;k<vs->irleng;k++)
661 en += vs->ir[k] * vs->ir[k];
662
663 return(en);
664}
665
666
667// b2bc : transform MLSA digital filter coefficients to mel-cepstrum
668static void b2mc (double *b, double *mc, int m, double a)
669{
670 double d, o;
671
672 d = mc[m] = b[m];
673 for (m--; m>=0; m--) {
674 o = b[m] + a * d;
675 d = b[m];
676 mc[m] = o;
677 }
678
679 return;
680}
681
682// freqt : frequency transformation
683static void freqt (double *c1, int m1, double *c2, int m2, double a, VocoderSetup *vs)
684{
685 register int i, j;
686 double b;
687
688 if (vs->d==NULL__null) {
689 vs->size = m2;
690 vs->d = wcalloc(double,vs->size + vs->size + 2)((double *)safe_wcalloc(sizeof(double)*(vs->size + vs->
size + 2)));
691 vs->g = vs->d+vs->size+1;
692 }
693
694 if (m2>vs->size) {
695 wfree(vs->d);
696 vs->size = m2;
697 vs->d = wcalloc(double,vs->size + vs->size + 2)((double *)safe_wcalloc(sizeof(double)*(vs->size + vs->
size + 2)));
698 vs->g = vs->d+vs->size+1;
699 }
700
701 b = 1-a*a;
702 for (i=0; i<m2+1; i++)
703 vs->g[i] = 0.0;
704
705 for (i=-m1; i<=0; i++) {
706 if (0 <= m2)
707 vs->g[0] = c1[-i]+a*(vs->d[0]=vs->g[0]);
708 if (1 <= m2)
709 vs->g[1] = b*vs->d[0]+a*(vs->d[1]=vs->g[1]);
710 for (j=2; j<=m2; j++)
711 vs->g[j] = vs->d[j-1]+a*((vs->d[j]=vs->g[j])-vs->g[j-1]);
712 }
713
714 memmove(c2,vs->g,sizeof(double)*(m2+1));
715
716 return;
717}
718
719// c2ir : The minimum phase impulse response is evaluated from the minimum phase cepstrum
720static void c2ir (double *c, int nc, double *h, int leng)
721{
722 register int n, k, upl;
723 double d;
724
725 h[0] = exp(c[0]);
726 for (n=1; n<leng; n++) {
727 d = 0;
728 upl = (n>=nc) ? nc-1 : n;
729 for (k=1; k<=upl; k++)
730 d += k*c[k]*h[n-k];
731 h[n] = d/n;
732 }
733
734 return;
735}
736
737#if 0
738static double get_dpow(double *rmcep, double *cmcep, int m, double a,
739 VocoderSetup *vs)
740{
741 double e1, e2, dpow;
742
743 if (vs->p1 < 0) {
744 vs->p1 = 1;
745 vs->cc = vs->c + m + 1;
746 vs->cinc = vs->cc + m + 1;
747 vs->d1 = vs->cinc + m + 1;
748 }
749
750 mc2b(rmcep, vs->c, m, a);
751 e1 = b2en(vs->c, m, a, vs);
752
753 mc2b(cmcep, vs->cc, m, a);
754 e2 = b2en(vs->cc, m, a, vs);
755
756 dpow = log(e1 / e2) / 2.0;
757
758 return dpow;
759}
760#endif
761
762static void free_vocoder(VocoderSetup *vs)
763{
764 int i;
765
766 wfree(vs->c);
767 wfree(vs->mc);
768 wfree(vs->d);
769
770 vs->c = NULL__null;
771 vs->mc = NULL__null;
772 vs->d = NULL__null;
773 vs->ppade = NULL__null;
774 vs->cc = NULL__null;
775 vs->cinc = NULL__null;
776 vs->d1 = NULL__null;
777 vs->g = NULL__null;
778 vs->cep = NULL__null;
779 vs->ir = NULL__null;
780
781 wfree(vs->hpulse);
782 wfree(vs->hnoise);
783 wfree(vs->xpulsesig);
784 wfree(vs->xnoisesig);
785 for (i=0; i<vs->ME_num; i++)
786 wfree(vs->h[i]);
787 wfree(vs->h);
788
789 return;
790}
791
792/* from vector.cc */
793
794static DVECTOR xdvalloc(long length)
795{
796 DVECTOR x;
797
798 length = MAX(length, 0)((length) > (0) ? (length) : (0));
799 x = wcalloc(struct DVECTOR_STRUCT,1)((struct DVECTOR_STRUCT *)safe_wcalloc(sizeof(struct DVECTOR_STRUCT
)*(1)));
800 x->data = wcalloc(double,MAX(length, 1))((double *)safe_wcalloc(sizeof(double)*(((length) > (1) ? (
length) : (1)))));
801 x->imag = NULL__null;
802 x->length = length;
803
804 return x;
805}
806
807static void xdvfree(DVECTOR x)
808{
809 if (x != NULL__null) {
810 if (x->data != NULL__null) {
811 wfree(x->data);
812 }
813 if (x->imag != NULL__null) {
814 wfree(x->imag);
815 }
816 wfree(x);
817 }
818
819 return;
820}
821
822static void dvialloc(DVECTOR x)
823{
824 if (x->imag != NULL__null) {
825 wfree(x->imag);
826 }
827 x->imag = wcalloc(double,x->length)((double *)safe_wcalloc(sizeof(double)*(x->length)));
828
829 return;
830}
831
832static DVECTOR xdvcut(DVECTOR x, long offset, long length)
833{
834 long k;
835 long pos;
836 DVECTOR y;
837
838 y = xdvalloc(length);
839 if (x->imag != NULL__null) {
840 dvialloc(y);
841 }
842
843 for (k = 0; k < y->length; k++) {
844 pos = k + offset;
845 if (pos >= 0 && pos < x->length) {
846 y->data[k] = x->data[pos];
847 if (y->imag != NULL__null) {
848 y->imag[k] = x->imag[pos];
849 }
850 } else {
851 y->data[k] = 0.0;
852 if (y->imag != NULL__null) {
853 y->imag[k] = 0.0;
854 }
855 }
856 }
857
858 return y;
859}
860
861static DMATRIX xdmalloc(long row, long col)
862{
863 DMATRIX matrix;
864 int i;
865
866 matrix = wcalloc(struct DMATRIX_STRUCT,1)((struct DMATRIX_STRUCT *)safe_wcalloc(sizeof(struct DMATRIX_STRUCT
)*(1)));
867 matrix->data = wcalloc(double *,row)((double * *)safe_wcalloc(sizeof(double *)*(row)));
868 for (i=0; i<row; i++)
869 matrix->data[i] = wcalloc(double,col)((double *)safe_wcalloc(sizeof(double)*(col)));
870 matrix->imag = NULL__null;
871 matrix->row = row;
872 matrix->col = col;
873
874 return matrix;
875}
876
877void xdmfree(DMATRIX matrix)
878{
879 int i;
880
881 if (matrix != NULL__null) {
882 if (matrix->data != NULL__null) {
883 for (i=0; i<matrix->row; i++)
884 wfree(matrix->data[i]);
885 wfree(matrix->data);
886 }
887 if (matrix->imag != NULL__null) {
888 for (i=0; i<matrix->row; i++)
889 wfree(matrix->imag[i]);
890 wfree(matrix->imag);
891 }
892 wfree(matrix);
893 }
894
895 return;
896}
897
898
899/* from voperate.cc */
900static double dvmax(DVECTOR x, long *index)
901{
902 long k;
903 long ind;
904 double max;
905
906 ind = 0;
907 max = x->data[ind];
908 for (k = 1; k < x->length; k++) {
909 if (max < x->data[k]) {
910 ind = k;
911 max = x->data[k];
912 }
913 }
914
915 if (index != NULL__null) {
916 *index = ind;
917 }
918
919 return max;
920}
921
922static double dvmin(DVECTOR x, long *index)
923{
924 long k;
925 long ind;
926 double min;
927
928 ind = 0;
929 min = x->data[ind];
930 for (k = 1; k < x->length; k++) {
931 if (min > x->data[k]) {
932 ind = k;
933 min = x->data[k];
934 }
935 }
936
937 if (index != NULL__null) {
938 *index = ind;
939 }
940
941 return min;
942}