| File: | modules/hts21_engine/hts21_mlsa_resynthesis.cc | 
| Location: | line 461, column 4 | 
| Description: | Value stored to 'aa' is never read | 
| 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 fetvox 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 "hts21_mlsa_resynthesis.h" | 
| 64 | |
| 65 | LISP hts21_mlsa_resynthesis(LISP ltrack) | 
| 66 | { | 
| 67 | /* Resynthesizes a wave from given track */ | 
| 68 | EST_Track *t; | 
| 69 | EST_Wave *wave = 0; | 
| 70 | DVECTOR w; | 
| 71 | DMATRIX mcep; | 
| 72 | DVECTOR f0v; | 
| 73 | int sr = 16000; | 
| 74 | int i,j; | 
| 75 | double shift; | 
| 76 | |
| 77 | if ((ltrack == NULL__null) || | 
| 78 |         (TYPEP(ltrack,tc_string)( (ltrack != __null) && ((((ltrack) == ((struct obj * ) 0)) ? 0 : ((*(ltrack)).type)) == (13)) ) &&  | 
| 79 | (streq(get_c_string(ltrack),"nil")(strcmp(get_c_string(ltrack),"nil")==0)))) | 
| 80 | return siod(new EST_Wave(0,1,sr)); | 
| 81 | |
| 82 | t = track(ltrack); | 
| 83 | |
| 84 | f0v = xdvalloc(t->num_frames()); | 
| 85 | mcep = xdmalloc(t->num_frames(),t->num_channels()-1); | 
| 86 | |
| 87 | for (i=0; i<t->num_frames(); i++) | 
| 88 | { | 
| 89 | f0v->data[i] = t->a(i,0); | 
| 90 | for (j=1; j<t->num_channels(); j++) | 
| 91 | mcep->data[i][j-1] = t->a(i,j); | 
| 92 | } | 
| 93 | |
| 94 | if (t->num_frames() > 1) | 
| 95 | shift = 1000.0*(t->t(1)-t->t(0)); | 
| 96 | else | 
| 97 | shift = 5.0; | 
| 98 | |
| 99 | w = synthesis_body(mcep,f0v,NULL__null,sr,shift); | 
| 100 | |
| 101 | wave = new EST_Wave(w->length,1,sr); | 
| 102 | |
| 103 | for (i=0; i<w->length; i++) | 
| 104 | wave->a(i) = (short)w->data[i]; | 
| 105 | |
| 106 | xdmfree(mcep); | 
| 107 | xdvfree(f0v); | 
| 108 | xdvfree(w); | 
| 109 | |
| 110 | return siod(wave); | 
| 111 | } | 
| 112 | |
| 113 | |
| 114 | DVECTOR synthesis_body(DMATRIX mcep, // input mel-cep sequence | 
| 115 | DVECTOR f0v, // input F0 sequence | 
| 116 | DVECTOR dpow, // input diff-pow sequence | 
| 117 | double fs, // sampling frequency (Hz) | 
| 118 | double framem) // FFT length | 
| 119 | { | 
| 120 | long t, pos; | 
| 121 | int framel; | 
| 122 | double f0; | 
| 123 | VocoderSetup vs; | 
| 124 | DVECTOR xd = NODATA__null; | 
| 125 | DVECTOR syn = NODATA__null; | 
| 126 | |
| 127 | framel = (int)(framem * fs / 1000.0); | 
| 128 | init_vocoder(fs, framel, mcep->col - 1, &vs); | 
| 129 | |
| 130 | // synthesize waveforms by MLSA filter | 
| 131 | xd = xdvalloc(mcep->row * (framel + 2)); | 
| 132 | for (t = 0, pos = 0; t < mcep->row; t++) { | 
| 133 | if (t >= f0v->length) f0 = 0.0; | 
| 134 | else f0 = f0v->data[t]; | 
| 135 | if (dpow == NODATA__null) | 
| 136 | vocoder(f0, mcep->data[t], mcep->col - 1, ALPHA0.42, 0.0, &vs, | 
| 137 | xd->data, &pos); | 
| 138 | else | 
| 139 | vocoder(f0, mcep->data[t], dpow->data[t], mcep->col - 1, ALPHA0.42, | 
| 140 | 0.0, &vs, xd->data, &pos); | 
| 141 | } | 
| 142 | syn = xdvcut(xd, 0, pos); | 
| 143 | |
| 144 | // normalized amplitude | 
| 145 | waveampcheck(syn, XFALSE0); | 
| 146 | |
| 147 | // memory free | 
| 148 | xdvfree(xd); | 
| 149 | free_vocoder(&vs); | 
| 150 | |
| 151 | return syn; | 
| 152 | } | 
| 153 | |
| 154 | #if 0 | 
| 155 | static DVECTOR get_dpowvec(DMATRIX rmcep, DMATRIX cmcep) | 
| 156 | { | 
| 157 | long t; | 
| 158 | DVECTOR dpow = NODATA__null; | 
| 159 | VocoderSetup pvs; | 
| 160 | |
| 161 | // error check | 
| 162 | if (rmcep->col != cmcep->col) { | 
| 163 | fprintf(stderrstderr, "Error: Different number of dimensions\n"); | 
| 164 | exit(1); | 
| 165 | } | 
| 166 | if (rmcep->row != cmcep->row) { | 
| 167 | fprintf(stderrstderr, "Error: Different number of frames\n"); | 
| 168 | exit(1); | 
| 169 | } | 
| 170 | |
| 171 | // memory allocation | 
| 172 | dpow = xdvalloc(rmcep->row); | 
| 173 | init_vocoder(16000.0, 80, rmcep->col - 1, &pvs); | 
| 174 | |
| 175 | // calculate differential power | 
| 176 | for (t = 0; t < rmcep->row; t++) | 
| 177 | dpow->data[t] = get_dpow(rmcep->data[t], cmcep->data[t], | 
| 178 | rmcep->col - 1, ALPHA0.42, &pvs); | 
| 179 | |
| 180 | // memory free | 
| 181 | free_vocoder(&pvs); | 
| 182 | |
| 183 | return dpow; | 
| 184 | } | 
| 185 | #endif | 
| 186 | |
| 187 | static void waveampcheck(DVECTOR wav, XBOOLint msg_flag) | 
| 188 | { | 
| 189 | double value; | 
| 190 | int k; | 
| 191 | |
| 192 |     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)))));  | 
| 193 | if (value >= 32000.0) { | 
| 194 | if (msg_flag == XTRUE1) { | 
| 195 | fprintf(stderrstderr, "amplitude is too big: %f\n", value); | 
| 196 | fprintf(stderrstderr, "execute normalization\n"); | 
| 197 | } | 
| 198 | /* was dvscoper(wav, "*", 32000.0 / value); */ | 
| 199 | for (k = 0; k < wav->length; k++) { | 
| 200 | wav->data[k] = wav->data[k] * (32000.0/value); | 
| 201 | if (wav->imag != NULL__null) { | 
| 202 | wav->imag[k] = wav->imag[k] * (32000.0/value); | 
| 203 | } | 
| 204 | } | 
| 205 | } | 
| 206 | |
| 207 | return; | 
| 208 | } | 
| 209 | |
| 210 | static void init_vocoder(double fs, int framel, int m, VocoderSetup *vs) | 
| 211 | { | 
| 212 | // initialize global parameter | 
| 213 | vs->fprd = framel; | 
| 214 | vs->iprd = 1; | 
| 215 | vs->seed = 1; | 
| 216 | vs->pd = 5; | 
| 217 | |
| 218 | vs->next =1; | 
| 219 | vs->gauss = MTRUE; | 
| 220 | |
| 221 | vs->pade[ 0]=1.0; | 
| 222 | vs->pade[ 1]=1.0; vs->pade[ 2]=0.0; | 
| 223 | vs->pade[ 3]=1.0; vs->pade[ 4]=0.0; vs->pade[ 5]=0.0; | 
| 224 | vs->pade[ 6]=1.0; vs->pade[ 7]=0.0; vs->pade[ 8]=0.0; vs->pade[ 9]=0.0; | 
| 225 | 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; | 
| 226 | 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; | 
| 227 | vs->pade[20]=0.00003041721; | 
| 228 | |
| 229 | vs->rate = fs; | 
| 230 |     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))));  | 
| 231 | |
| 232 | vs->p1 = -1; | 
| 233 | vs->sw = 0; | 
| 234 | vs->x = 0x55555555; | 
| 235 | |
| 236 | // for postfiltering | 
| 237 | vs->mc = NULL__null; | 
| 238 | vs->o = 0; | 
| 239 | vs->d = NULL__null; | 
| 240 | vs->irleng= 64; | 
| 241 | |
| 242 | return; | 
| 243 | } | 
| 244 | |
| 245 | static void vocoder(double p, double *mc, int m, double a, double beta, | 
| 246 | VocoderSetup *vs, double *wav, long *pos) | 
| 247 | { | 
| 248 | double inc, x, e1, e2; | 
| 249 | int i, j, k; | 
| 250 | |
| 251 | if (p != 0.0) | 
| 252 | p = vs->rate / p; // f0 -> pitch | 
| 253 | |
| 254 | if (vs->p1 < 0) { | 
| 255 | if (vs->gauss & (vs->seed != 1)) | 
| 256 | vs->next = srnd((unsigned)vs->seed); | 
| 257 | |
| 258 | vs->p1 = p; | 
| 259 | vs->pc = vs->p1; | 
| 260 | vs->cc = vs->c + m + 1; | 
| 261 | vs->cinc = vs->cc + m + 1; | 
| 262 | vs->d1 = vs->cinc + m + 1; | 
| 263 | |
| 264 | mc2b(mc, vs->c, m, a); | 
| 265 | |
| 266 | if (beta > 0.0 && m > 1) { | 
| 267 | e1 = b2en(vs->c, m, a, vs); | 
| 268 | vs->c[1] -= beta * a * mc[2]; | 
| 269 | for (k=2;k<=m;k++) | 
| 270 | vs->c[k] *= (1.0 + beta); | 
| 271 | e2 = b2en(vs->c, m, a, vs); | 
| 272 | vs->c[0] += log(e1/e2)/2; | 
| 273 | } | 
| 274 | |
| 275 | return; | 
| 276 | } | 
| 277 | |
| 278 | mc2b(mc, vs->cc, m, a); | 
| 279 | if (beta>0.0 && m > 1) { | 
| 280 | e1 = b2en(vs->cc, m, a, vs); | 
| 281 | vs->cc[1] -= beta * a * mc[2]; | 
| 282 | for (k = 2; k <= m; k++) | 
| 283 | vs->cc[k] *= (1.0 + beta); | 
| 284 | e2 = b2en(vs->cc, m, a, vs); | 
| 285 | vs->cc[0] += log(e1 / e2) / 2.0; | 
| 286 | } | 
| 287 | |
| 288 | for (k=0; k<=m; k++) | 
| 289 | vs->cinc[k] = (vs->cc[k] - vs->c[k]) * | 
| 290 | (double)vs->iprd / (double)vs->fprd; | 
| 291 | |
| 292 | if (vs->p1!=0.0 && p!=0.0) { | 
| 293 | inc = (p - vs->p1) * (double)vs->iprd / (double)vs->fprd; | 
| 294 | } else { | 
| 295 | inc = 0.0; | 
| 296 | vs->pc = p; | 
| 297 | vs->p1 = 0.0; | 
| 298 | } | 
| 299 | |
| 300 | for (j = vs->fprd, i = (vs->iprd + 1) / 2; j--;) { | 
| 301 | if (vs->p1 == 0.0) { | 
| 302 | if (vs->gauss) | 
| 303 | x = (double) nrandom(vs); | 
| 304 | else | 
| 305 | x = mseq(vs); | 
| 306 | } else { | 
| 307 | if ((vs->pc += 1.0) >= vs->p1) { | 
| 308 | x = sqrt (vs->p1); | 
| 309 | vs->pc = vs->pc - vs->p1; | 
| 310 | } else x = 0.0; | 
| 311 | } | 
| 312 | |
| 313 | x *= exp(vs->c[0]); | 
| 314 | |
| 315 | x = mlsadf(x, vs->c, m, a, vs->pd, vs->d1, vs); | 
| 316 | |
| 317 | wav[*pos] = x; | 
| 318 | *pos += 1; | 
| 319 | |
| 320 | if (!--i) { | 
| 321 | vs->p1 += inc; | 
| 322 | for (k = 0; k <= m; k++) vs->c[k] += vs->cinc[k]; | 
| 323 | i = vs->iprd; | 
| 324 | } | 
| 325 | } | 
| 326 | |
| 327 | vs->p1 = p; | 
| 328 | memmove(vs->c,vs->cc,sizeof(double)*(m+1)); | 
| 329 | |
| 330 | return; | 
| 331 | } | 
| 332 | |
| 333 | static void vocoder(double p, double *mc, double dpow, int m, double a, double beta, | 
| 334 | VocoderSetup *vs, double *wav, long *pos) | 
| 335 | { | 
| 336 | double inc, x, e1, e2; | 
| 337 | int i, j, k; | 
| 338 | |
| 339 | if (p != 0.0) | 
| 340 | p = vs->rate / p; // f0 -> pitch | 
| 341 | |
| 342 | if (vs->p1 < 0) { | 
| 343 | if (vs->gauss & (vs->seed != 1)) | 
| 344 | vs->next = srnd((unsigned)vs->seed); | 
| 345 | |
| 346 | vs->p1 = p; | 
| 347 | vs->pc = vs->p1; | 
| 348 | vs->cc = vs->c + m + 1; | 
| 349 | vs->cinc = vs->cc + m + 1; | 
| 350 | vs->d1 = vs->cinc + m + 1; | 
| 351 | |
| 352 | mc2b(mc, vs->c, m, a); | 
| 353 | vs->c[0] += dpow; | 
| 354 | |
| 355 | if (beta > 0.0 && m > 1) { | 
| 356 | e1 = b2en(vs->c, m, a, vs); | 
| 357 | vs->c[1] -= beta * a * mc[2]; | 
| 358 | for (k=2;k<=m;k++) | 
| 359 | vs->c[k] *= (1.0 + beta); | 
| 360 | e2 = b2en(vs->c, m, a, vs); | 
| 361 | vs->c[0] += log(e1/e2)/2; | 
| 362 | } | 
| 363 | |
| 364 | return; | 
| 365 | } | 
| 366 | |
| 367 | mc2b(mc, vs->cc, m, a); | 
| 368 | vs->cc[0] += dpow; | 
| 369 | if (beta>0.0 && m > 1) { | 
| 370 | e1 = b2en(vs->cc, m, a, vs); | 
| 371 | vs->cc[1] -= beta * a * mc[2]; | 
| 372 | for (k = 2; k <= m; k++) | 
| 373 | vs->cc[k] *= (1.0 + beta); | 
| 374 | e2 = b2en(vs->cc, m, a, vs); | 
| 375 | vs->cc[0] += log(e1 / e2) / 2.0; | 
| 376 | } | 
| 377 | |
| 378 | for (k=0; k<=m; k++) | 
| 379 | vs->cinc[k] = (vs->cc[k] - vs->c[k]) * | 
| 380 | (double)vs->iprd / (double)vs->fprd; | 
| 381 | |
| 382 | if (vs->p1!=0.0 && p!=0.0) { | 
| 383 | inc = (p - vs->p1) * (double)vs->iprd / (double)vs->fprd; | 
| 384 | } else { | 
| 385 | inc = 0.0; | 
| 386 | vs->pc = p; | 
| 387 | vs->p1 = 0.0; | 
| 388 | } | 
| 389 | |
| 390 | for (j = vs->fprd, i = (vs->iprd + 1) / 2; j--;) { | 
| 391 | if (vs->p1 == 0.0) { | 
| 392 | if (vs->gauss) | 
| 393 | x = (double) nrandom(vs); | 
| 394 | else | 
| 395 | x = mseq(vs); | 
| 396 | } else { | 
| 397 | if ((vs->pc += 1.0) >= vs->p1) { | 
| 398 | x = sqrt (vs->p1); | 
| 399 | vs->pc = vs->pc - vs->p1; | 
| 400 | } else x = 0.0; | 
| 401 | } | 
| 402 | |
| 403 | x *= exp(vs->c[0]); | 
| 404 | |
| 405 | x = mlsadf(x, vs->c, m, a, vs->pd, vs->d1, vs); | 
| 406 | |
| 407 | wav[*pos] = x; | 
| 408 | *pos += 1; | 
| 409 | |
| 410 | if (!--i) { | 
| 411 | vs->p1 += inc; | 
| 412 | for (k = 0; k <= m; k++) vs->c[k] += vs->cinc[k]; | 
| 413 | i = vs->iprd; | 
| 414 | } | 
| 415 | } | 
| 416 | |
| 417 | vs->p1 = p; | 
| 418 | memmove(vs->c,vs->cc,sizeof(double)*(m+1)); | 
| 419 | |
| 420 | return; | 
| 421 | } | 
| 422 | |
| 423 | static double mlsadf(double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs) | 
| 424 | { | 
| 425 | |
| 426 | vs->ppade = &(vs->pade[pd*(pd+1)/2]); | 
| 427 | |
| 428 | x = mlsadf1 (x, b, m, a, pd, d, vs); | 
| 429 | x = mlsadf2 (x, b, m, a, pd, &d[2*(pd+1)], vs); | 
| 430 | |
| 431 | return(x); | 
| 432 | } | 
| 433 | |
| 434 | static double mlsadf1(double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs) | 
| 435 | { | 
| 436 | double v, out = 0.0, *pt, aa; | 
| 437 | register int i; | 
| 438 | |
| 439 | aa = 1 - a*a; | 
| 440 | pt = &d[pd+1]; | 
| 441 | |
| 442 | for (i=pd; i>=1; i--) { | 
| 443 | d[i] = aa*pt[i-1] + a*d[i]; | 
| 444 | pt[i] = d[i] * b[1]; | 
| 445 | v = pt[i] * vs->ppade[i]; | 
| 446 | x += (1 & i) ? v : -v; | 
| 447 | out += v; | 
| 448 | } | 
| 449 | |
| 450 | pt[0] = x; | 
| 451 | out += x; | 
| 452 | |
| 453 | return(out); | 
| 454 | } | 
| 455 | |
| 456 | static double mlsadf2 (double x, double *b, int m, double a, int pd, double *d, VocoderSetup *vs) | 
| 457 | { | 
| 458 | double v, out = 0.0, *pt, aa; | 
| 459 | register int i; | 
| 460 | |
| 461 | aa = 1 - a*a; | 
Value stored to 'aa' is never read  | |
| 462 | pt = &d[pd * (m+2)]; | 
| 463 | |
| 464 | for (i=pd; i>=1; i--) { | 
| 465 | pt[i] = mlsafir (pt[i-1], b, m, a, &d[(i-1)*(m+2)]); | 
| 466 | v = pt[i] * vs->ppade[i]; | 
| 467 | |
| 468 | x += (1&i) ? v : -v; | 
| 469 | out += v; | 
| 470 | } | 
| 471 | |
| 472 | pt[0] = x; | 
| 473 | out += x; | 
| 474 | |
| 475 | return(out); | 
| 476 | } | 
| 477 | |
| 478 | static double mlsafir (double x, double *b, int m, double a, double *d) | 
| 479 | { | 
| 480 | double y = 0.0; | 
| 481 | double aa; | 
| 482 | register int i; | 
| 483 | |
| 484 | aa = 1 - a*a; | 
| 485 | |
| 486 | d[0] = x; | 
| 487 | d[1] = aa*d[0] + a*d[1]; | 
| 488 | |
| 489 | for (i=2; i<=m; i++) { | 
| 490 | d[i] = d[i] + a*(d[i+1]-d[i-1]); | 
| 491 | y += d[i]*b[i]; | 
| 492 | } | 
| 493 | |
| 494 | for (i=m+1; i>1; i--) | 
| 495 | d[i] = d[i-1]; | 
| 496 | |
| 497 | return(y); | 
| 498 | } | 
| 499 | |
| 500 | static double nrandom (VocoderSetup *vs) | 
| 501 | { | 
| 502 | if (vs->sw == 0) { | 
| 503 | vs->sw = 1; | 
| 504 | do { | 
| 505 | vs->r1 = 2.0 * rnd(&vs->next) - 1.0; | 
| 506 | vs->r2 = 2.0 * rnd(&vs->next) - 1.0; | 
| 507 | vs->s = vs->r1 * vs->r1 + vs->r2 * vs->r2; | 
| 508 | } while (vs->s > 1 || vs->s == 0); | 
| 509 | |
| 510 | vs->s = sqrt (-2 * log(vs->s) / vs->s); | 
| 511 | |
| 512 | return(vs->r1*vs->s); | 
| 513 | } | 
| 514 | else { | 
| 515 | vs->sw = 0; | 
| 516 | |
| 517 | return (vs->r2*vs->s); | 
| 518 | } | 
| 519 | } | 
| 520 | |
| 521 | static double rnd (unsigned long *next) | 
| 522 | { | 
| 523 | double r; | 
| 524 | |
| 525 | *next = *next * 1103515245L + 12345; | 
| 526 | r = (*next / 65536L) % 32768L; | 
| 527 | |
| 528 | return(r/RANDMAX32767); | 
| 529 | } | 
| 530 | |
| 531 | static unsigned long srnd ( unsigned long seed ) | 
| 532 | { | 
| 533 | return(seed); | 
| 534 | } | 
| 535 | |
| 536 | static int mseq (VocoderSetup *vs) | 
| 537 | { | 
| 538 | register int x0, x28; | 
| 539 | |
| 540 | vs->x >>= 1; | 
| 541 | |
| 542 | if (vs->x & B00x00000001) | 
| 543 | x0 = 1; | 
| 544 | else | 
| 545 | x0 = -1; | 
| 546 | |
| 547 | if (vs->x & B280x10000000) | 
| 548 | x28 = 1; | 
| 549 | else | 
| 550 | x28 = -1; | 
| 551 | |
| 552 | if (x0 + x28) | 
| 553 | vs->x &= B31_0x7fffffff; | 
| 554 | else | 
| 555 | vs->x |= B310x80000000; | 
| 556 | |
| 557 | return(x0); | 
| 558 | } | 
| 559 | |
| 560 | // mc2b : transform mel-cepstrum to MLSA digital fillter coefficients | 
| 561 | static void mc2b (double *mc, double *b, int m, double a) | 
| 562 | { | 
| 563 | b[m] = mc[m]; | 
| 564 | |
| 565 | for (m--; m>=0; m--) | 
| 566 | b[m] = mc[m] - a * b[m+1]; | 
| 567 | |
| 568 | return; | 
| 569 | } | 
| 570 | |
| 571 | |
| 572 | static double b2en (double *b, int m, double a, VocoderSetup *vs) | 
| 573 | { | 
| 574 | double en; | 
| 575 | int k; | 
| 576 | |
| 577 | if (vs->o<m) { | 
| 578 | if (vs->mc != NULL__null) | 
| 579 | wfree(vs->mc); | 
| 580 | |
| 581 |       vs->mc = wcalloc(double,(m + 1) + 2 * vs->irleng)((double *)safe_wcalloc(sizeof(double)*((m + 1) + 2 * vs-> irleng)));  | 
| 582 | vs->cep = vs->mc + m+1; | 
| 583 | vs->ir = vs->cep + vs->irleng; | 
| 584 | } | 
| 585 | |
| 586 | b2mc(b, vs->mc, m, a); | 
| 587 | freqt(vs->mc, m, vs->cep, vs->irleng-1, -a, vs); | 
| 588 | c2ir(vs->cep, vs->irleng, vs->ir, vs->irleng); | 
| 589 | en = 0.0; | 
| 590 | |
| 591 | for (k=0;k<vs->irleng;k++) | 
| 592 | en += vs->ir[k] * vs->ir[k]; | 
| 593 | |
| 594 | return(en); | 
| 595 | } | 
| 596 | |
| 597 | |
| 598 | // b2bc : transform MLSA digital filter coefficients to mel-cepstrum | 
| 599 | static void b2mc (double *b, double *mc, int m, double a) | 
| 600 | { | 
| 601 | double d, o; | 
| 602 | |
| 603 | d = mc[m] = b[m]; | 
| 604 | for (m--; m>=0; m--) { | 
| 605 | o = b[m] + a * d; | 
| 606 | d = b[m]; | 
| 607 | mc[m] = o; | 
| 608 | } | 
| 609 | |
| 610 | return; | 
| 611 | } | 
| 612 | |
| 613 | // freqt : frequency transformation | 
| 614 | static void freqt (double *c1, int m1, double *c2, int m2, double a, VocoderSetup *vs) | 
| 615 | { | 
| 616 | register int i, j; | 
| 617 | double b; | 
| 618 | |
| 619 | if (vs->d==NULL__null) { | 
| 620 | vs->size = m2; | 
| 621 |       vs->d    = wcalloc(double,vs->size + vs->size + 2)((double *)safe_wcalloc(sizeof(double)*(vs->size + vs-> size + 2)));  | 
| 622 | vs->g = vs->d+vs->size+1; | 
| 623 | } | 
| 624 | |
| 625 | if (m2>vs->size) { | 
| 626 | wfree(vs->d); | 
| 627 | vs->size = m2; | 
| 628 |       vs->d    = wcalloc(double,vs->size + vs->size + 2)((double *)safe_wcalloc(sizeof(double)*(vs->size + vs-> size + 2)));  | 
| 629 | vs->g = vs->d+vs->size+1; | 
| 630 | } | 
| 631 | |
| 632 | b = 1-a*a; | 
| 633 | for (i=0; i<m2+1; i++) | 
| 634 | vs->g[i] = 0.0; | 
| 635 | |
| 636 | for (i=-m1; i<=0; i++) { | 
| 637 | if (0 <= m2) | 
| 638 | vs->g[0] = c1[-i]+a*(vs->d[0]=vs->g[0]); | 
| 639 | if (1 <= m2) | 
| 640 | vs->g[1] = b*vs->d[0]+a*(vs->d[1]=vs->g[1]); | 
| 641 | for (j=2; j<=m2; j++) | 
| 642 | vs->g[j] = vs->d[j-1]+a*((vs->d[j]=vs->g[j])-vs->g[j-1]); | 
| 643 | } | 
| 644 | |
| 645 | memmove(c2,vs->g,sizeof(double)*(m2+1)); | 
| 646 | |
| 647 | return; | 
| 648 | } | 
| 649 | |
| 650 | // c2ir : The minimum phase impulse response is evaluated from the minimum phase cepstrum | 
| 651 | static void c2ir (double *c, int nc, double *h, int leng) | 
| 652 | { | 
| 653 | register int n, k, upl; | 
| 654 | double d; | 
| 655 | |
| 656 | h[0] = exp(c[0]); | 
| 657 | for (n=1; n<leng; n++) { | 
| 658 | d = 0; | 
| 659 | upl = (n>=nc) ? nc-1 : n; | 
| 660 | for (k=1; k<=upl; k++) | 
| 661 | d += k*c[k]*h[n-k]; | 
| 662 | h[n] = d/n; | 
| 663 | } | 
| 664 | |
| 665 | return; | 
| 666 | } | 
| 667 | |
| 668 | #if 0 | 
| 669 | static double get_dpow(double *rmcep, double *cmcep, int m, double a, | 
| 670 | VocoderSetup *vs) | 
| 671 | { | 
| 672 | double e1, e2, dpow; | 
| 673 | |
| 674 | if (vs->p1 < 0) { | 
| 675 | vs->p1 = 1; | 
| 676 | vs->cc = vs->c + m + 1; | 
| 677 | vs->cinc = vs->cc + m + 1; | 
| 678 | vs->d1 = vs->cinc + m + 1; | 
| 679 | } | 
| 680 | |
| 681 | mc2b(rmcep, vs->c, m, a); | 
| 682 | e1 = b2en(vs->c, m, a, vs); | 
| 683 | |
| 684 | mc2b(cmcep, vs->cc, m, a); | 
| 685 | e2 = b2en(vs->cc, m, a, vs); | 
| 686 | |
| 687 | dpow = log(e1 / e2) / 2.0; | 
| 688 | |
| 689 | return dpow; | 
| 690 | } | 
| 691 | #endif | 
| 692 | |
| 693 | static void free_vocoder(VocoderSetup *vs) | 
| 694 | { | 
| 695 | wfree(vs->c); | 
| 696 | wfree(vs->mc); | 
| 697 | wfree(vs->d); | 
| 698 | |
| 699 | vs->c = NULL__null; | 
| 700 | vs->mc = NULL__null; | 
| 701 | vs->d = NULL__null; | 
| 702 | vs->ppade = NULL__null; | 
| 703 | vs->cc = NULL__null; | 
| 704 | vs->cinc = NULL__null; | 
| 705 | vs->d1 = NULL__null; | 
| 706 | vs->g = NULL__null; | 
| 707 | vs->cep = NULL__null; | 
| 708 | vs->ir = NULL__null; | 
| 709 | |
| 710 | return; | 
| 711 | } | 
| 712 | |
| 713 | /* from vector.cc */ | 
| 714 | |
| 715 | static DVECTOR xdvalloc(long length) | 
| 716 | { | 
| 717 | DVECTOR x; | 
| 718 | |
| 719 | length = MAX(length, 0)((length) > (0) ? (length) : (0)); | 
| 720 |     x = wcalloc(struct DVECTOR_STRUCT,1)((struct DVECTOR_STRUCT *)safe_wcalloc(sizeof(struct DVECTOR_STRUCT )*(1)));  | 
| 721 |     x->data = wcalloc(double,MAX(length, 1))((double *)safe_wcalloc(sizeof(double)*(((length) > (1) ? ( length) : (1)))));  | 
| 722 | x->imag = NULL__null; | 
| 723 | x->length = length; | 
| 724 | |
| 725 | return x; | 
| 726 | } | 
| 727 | |
| 728 | static void xdvfree(DVECTOR x) | 
| 729 | { | 
| 730 | if (x != NULL__null) { | 
| 731 | if (x->data != NULL__null) { | 
| 732 | wfree(x->data); | 
| 733 | } | 
| 734 | if (x->imag != NULL__null) { | 
| 735 | wfree(x->imag); | 
| 736 | } | 
| 737 | wfree(x); | 
| 738 | } | 
| 739 | |
| 740 | return; | 
| 741 | } | 
| 742 | |
| 743 | static void dvialloc(DVECTOR x) | 
| 744 | { | 
| 745 | if (x->imag != NULL__null) { | 
| 746 | wfree(x->imag); | 
| 747 | } | 
| 748 | x->imag = wcalloc(double,x->length)((double *)safe_wcalloc(sizeof(double)*(x->length))); | 
| 749 | |
| 750 | return; | 
| 751 | } | 
| 752 | |
| 753 | static DVECTOR xdvcut(DVECTOR x, long offset, long length) | 
| 754 | { | 
| 755 | long k; | 
| 756 | long pos; | 
| 757 | DVECTOR y; | 
| 758 | |
| 759 | y = xdvalloc(length); | 
| 760 | if (x->imag != NULL__null) { | 
| 761 | dvialloc(y); | 
| 762 | } | 
| 763 | |
| 764 | for (k = 0; k < y->length; k++) { | 
| 765 | pos = k + offset; | 
| 766 | if (pos >= 0 && pos < x->length) { | 
| 767 | y->data[k] = x->data[pos]; | 
| 768 | if (y->imag != NULL__null) { | 
| 769 | y->imag[k] = x->imag[pos]; | 
| 770 | } | 
| 771 | } else { | 
| 772 | y->data[k] = 0.0; | 
| 773 | if (y->imag != NULL__null) { | 
| 774 | y->imag[k] = 0.0; | 
| 775 | } | 
| 776 | } | 
| 777 | } | 
| 778 | |
| 779 | return y; | 
| 780 | } | 
| 781 | |
| 782 | static DMATRIX xdmalloc(long row, long col) | 
| 783 | { | 
| 784 | DMATRIX matrix; | 
| 785 | int i; | 
| 786 | |
| 787 |     matrix = wcalloc(struct DMATRIX_STRUCT,1)((struct DMATRIX_STRUCT *)safe_wcalloc(sizeof(struct DMATRIX_STRUCT )*(1)));  | 
| 788 | matrix->data = wcalloc(double *,row)((double * *)safe_wcalloc(sizeof(double *)*(row))); | 
| 789 | for (i=0; i<row; i++) | 
| 790 | matrix->data[i] = wcalloc(double,col)((double *)safe_wcalloc(sizeof(double)*(col))); | 
| 791 | matrix->imag = NULL__null; | 
| 792 | matrix->row = row; | 
| 793 | matrix->col = col; | 
| 794 | |
| 795 | return matrix; | 
| 796 | } | 
| 797 | |
| 798 | void xdmfree(DMATRIX matrix) | 
| 799 | { | 
| 800 | int i; | 
| 801 | |
| 802 | if (matrix != NULL__null) { | 
| 803 | if (matrix->data != NULL__null) { | 
| 804 | for (i=0; i<matrix->row; i++) | 
| 805 | wfree(matrix->data[i]); | 
| 806 | wfree(matrix->data); | 
| 807 | } | 
| 808 | if (matrix->imag != NULL__null) { | 
| 809 | for (i=0; i<matrix->row; i++) | 
| 810 | wfree(matrix->imag[i]); | 
| 811 | wfree(matrix->imag); | 
| 812 | } | 
| 813 | wfree(matrix); | 
| 814 | } | 
| 815 | |
| 816 | return; | 
| 817 | } | 
| 818 | |
| 819 | |
| 820 | /* from voperate.cc */ | 
| 821 | static double dvmax(DVECTOR x, long *index) | 
| 822 | { | 
| 823 | long k; | 
| 824 | long ind; | 
| 825 | double max; | 
| 826 | |
| 827 | ind = 0; | 
| 828 | max = x->data[ind]; | 
| 829 | for (k = 1; k < x->length; k++) { | 
| 830 | if (max < x->data[k]) { | 
| 831 | ind = k; | 
| 832 | max = x->data[k]; | 
| 833 | } | 
| 834 | } | 
| 835 | |
| 836 | if (index != NULL__null) { | 
| 837 | *index = ind; | 
| 838 | } | 
| 839 | |
| 840 | return max; | 
| 841 | } | 
| 842 | |
| 843 | static double dvmin(DVECTOR x, long *index) | 
| 844 | { | 
| 845 | long k; | 
| 846 | long ind; | 
| 847 | double min; | 
| 848 | |
| 849 | ind = 0; | 
| 850 | min = x->data[ind]; | 
| 851 | for (k = 1; k < x->length; k++) { | 
| 852 | if (min > x->data[k]) { | 
| 853 | ind = k; | 
| 854 | min = x->data[k]; | 
| 855 | } | 
| 856 | } | 
| 857 | |
| 858 | if (index != NULL__null) { | 
| 859 | *index = ind; | 
| 860 | } | 
| 861 | |
| 862 | return min; | 
| 863 | } |