
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                        P A R W A V E . C                            */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*   Copyright            1982                    by Dennis H. Klatt
 *
 *      Klatt synthesizer
 *         Modified version of synthesizer described in
 *         J. Acoust. Soc. Am., Mar. 1980. -- new voicing
 *         source.
 *
 * Edit history
 * 000001 10-Mar-83 DK  Initial creation.
 * 000002  5-May-83 DK  Fix bug in operation of parallel F1
 * 000003  7-Jul-83 DK  Allow parallel B1 to vary, and if ALL_PARALLEL,
 *                      also allow B2 and B3 to vary
 * 000004 26-Jul-83 DK  Get rid of mulsh, use short for VAX
 * 000005 24-Oct-83 DK  Split off parwavtab.c, change short to int
 * 000006 16-Nov-83 DK  Make samrate a variable, use exp(), cos() rand()
 * 000007 17-Nov-83 DK  Convert to float, remove  cpsw, add set outsl
 * 000008 28-Nov-83 DK  Add simple impulsive glottal source option
 * 000009  7-Dec-83 DK  Use spkrdef[7] to select impulse or natural voicing
 *                       and update cascade F1,..,F6 at update times
 * 000010 19-Dec-83 DK  Add subroutine no_rad_char() to get rid of rad char
 * 000011 28-Jan-84 DK  Allow up to 8 formants in cascade branch F7 fixed
 *                       at 6.5 kHz, F8 fixed at 7.5 kHz
 * 000012 14-Feb-84 DK  Fix bug in 'os' options so os>12 works
 * 000013 17-May-84 DK  Add G0 code
 * 000014 12-Mar-85 DHW modify for Haskins environment
 * 000015 11-Jul-87 LG  modificiations for PC
 * 000016 20-Apr-91 ATS Modified for SPARCSTATION
 */

#include <stdio.h>
#include <math.h>
#include "proto.h"
#include "parwavta.def"

extern int rand();
extern int srand();

long truncate();
float DBtoLIN();
float impulsive_source();
float natural_source();


/*
function FLUTTER

This function adds F0 flutter, as specified in:

"Analysis, synthesis and perception of voice quality variations among
female and male talkers" D.H. Klatt and L.C. Klatt JASA 87(2) February 1990.
Flutter is added by applying a quasi-random element constructed from three
slowly varying sine waves.
*/

void flutter()

{
  long original_f0;
  double delta_f0;
  double fla,flb,flc,fld,fle;

  original_f0 = pars[0] / 10;

  fla = (double) f0_flutter / 50;
  flb = (double) original_f0 / 100;
  flc = sin(2*PI*12.7*time_count);
  fld = sin(2*PI*7.1*time_count);
  fle = sin(2*PI*4.7*time_count);
  delta_f0 =  fla * flb * (flc + fld + fle) * 10;
  F0hz10 = F0hz10 + (long) delta_f0;
}




/*
function PARWAV

CONVERT FRAME OF PARAMETER DATA TO A WAVEFORM CHUNK
Synthesize nspfr samples of waveform and store in jwave[].
*/


void
parwav(jwave)
short *jwave;

{


/* Initialize synthesizer and get specification for current speech
   frame from host microcomputer */

  gethost();

  time_count++; /* used for f0 flutter */
  flutter();    /* add f0 flutter */

/* MAIN LOOP, for each output sample of current frame: */

  for (ns=0; ns<nspfr; ns++)
  {

/*    Get low-passed random number for aspiration and frication noise */
    gen_noise();                    /* output variable is 'noise' */

/*    Amplitude modulate noise (reduce noise amplitude during
 *    second half of glottal period) if voicing simultaneously present */
    if (nper > nmod)
    {
      noise *= 0.5;
    }

/*    Compute frication noise */
    frics = amp_frica * noise;

/*  Compute voicing waveform */
/*    (run glottal source simulation at 4 times normal sample rate to minimize
 *    quantization noise in period of female voice) */

    for (n4=0; n4<4; n4++)
    {

/*        Use impulsive glottal source */
      if (glsource == IMPULSIVE)
      {
	voice = impulsive_source();
      }
/*        Or use a more-natural-shaped source waveform with excitation
          occurring both upon opening and upon closure, stronest at closure */
      else
      {
	if (glsource == NATURAL)
	{
	  voice = natural_source();
	}
      }

/*        Reset period when counter 'nper' reaches T0 */
      if (nper >= T0)
      {
	nper = 0;
	pitch_synch_par_reset();
      }


/*        Low-pass filter voicing waveform before downsampling from 4*samrate
 *        to samrate samples/sec.  Resonator f=.09*samrate, bw=.06*samrate */
      resonlp();          /* in=voice, out=voice */

/*        Increment counter that keeps track of 4*samrate samples/sec */
      nper++;
    }

/*    Tilt spectrum of voicing source down by soft low-pass filtering, amount
 *    of tilt determined by TLTdb */
    voice = (voice * onemd) + (vlast * decay);
    vlast = voice;

/*    Add breathiness during glottal open phase */
    if (nper < nopen)
    {
/*        Amount of breathiness determined by parameter Aturb */
/*        Use nrand rather than noise because noise is low-passed */
      voice += amp_breth * nrand;
    }


/*    Set amplitude of voicing */
    glotout = amp_voice * voice;
    par_glotout = par_amp_voice * voice;


/*    Compute aspiration amplitude and add to voicing source */
    aspiration = amp_aspir * noise;
    glotout += aspiration;
    par_glotout += aspiration;


/*  Cascade vocal tract, excited by laryngeal sources.
 *    Nasal antiresonator, then formants FNP, F5, F4, F3, F2, F1 */

    if(synthesis_model != ALL_PARALLEL)
    {
      resoncnz();             /* in=glotout, out=rnzout   */
      resoncnp();             /* in=rnzout, out=rnpc_1 */

      casc_next_in = rnpc_1;

      if (nfcascade >= 8)
      {
	resonc8();          /* Do not use unless samrat = 16000 */
	casc_next_in = r8c_1;
      }

      if (nfcascade >= 7)
      {
	resonc7();          /* Do not use unless samrat = 16000 */
	casc_next_in = r7c_1;
      }

      if (nfcascade >= 6)
      {
	resonc6();          /* Do not use unless long vocal tract or */
	casc_next_in = r6c_1;       /* samrat increased */
      }

      if (nfcascade >= 5)
      {
	resonc5();
	casc_next_in = r5c_1;
      }

      if (nfcascade >= 4)
      {
	resonc4();
	casc_next_in = r4c_1;
      }

      if (nfcascade >= 3)
      {
	resonc3();
	casc_next_in = r3c_1;
      }

      if (nfcascade >= 2)
      {
	resonc2();
	casc_next_in = r2c_1;
      }

      if (nfcascade >= 1)
      {
	resonc1();
      }

      if (outsl > 0)
      {
	if (outsl < 13)
	{
	  if (outsl ==  1)        out = voice;
	  if (outsl ==  2)        out = aspiration;
	  if (outsl ==  3)        out = frics;
	  if (outsl ==  4)        out = glotout * 2.0;
	  if (outsl ==  5)        out = par_glotout;
	  if (outsl ==  6)        out = rnzout;
	  if (outsl ==  7)        out = rnpc_1;
	  if (outsl ==  8)        out = r5c_1;
	  if (outsl ==  9)        out = r4c_1;
	  if (outsl == 10)        out = r3c_1;
	  if (outsl == 11)        out = r2c_1;
	  if (outsl == 12)        out = r1c_1;
	  if (outsl <= 4)
	  {
	    no_rad_char(out);   /* Take out effects of radiation char */
	  }
	}
      }
      else
      {
	out = r1c_1;
      }
    }
    else
    {
      out = 0;
    }

/*    Excite parallel F1 and FNP by voicing waveform */

    sourc = par_glotout;        /* Source is voicing plus aspiration */
    reson1p();                 /* in=sourc, out=rnpp_1 */
    resonnpp();                /* in=sourc, out=r1c_1 */
    out += rnpp_1 + r1p_1;     /* Add in phase, boost lows for nasalized */


/*  Standard parallel vocal tract
 *  Formants F6,F5,F4,F3,F2, outputs added with alternating sign */

/*    Sound sourc for other parallel resonators is frication plus */
/*    first difference of voicing waveform */

    sourc = frics + par_glotout - glotlast;
    glotlast = par_glotout;


    reson6p();              /* in=sourc, out=r6p_1 */
    out = r6p_1 - out;

    reson5p();              /* in=sourc, out=r5p_1 */
    out = r5p_1 - out;

    reson4p();              /* in=sourc, out=r4p_1 */
    out = r4p_1 - out;

    reson3p();              /* in=sourc, out=r3p_1 */
    out = r3p_1 - out;

    reson2p();              /* in=sourc, out=r2p_1 */
    out = r2p_1 - out;

    outbypas = amp_bypas * sourc;
    out = outbypas - out;

    if (outsl > 12)
    {
      if (outsl == 13)        out = r6p_1;
      if (outsl == 14)        out = r5p_1;
      if (outsl == 15)        out = r4p_1;
      if (outsl == 16)        out = r3p_1;
      if (outsl == 17)        out = r2p_1;
      if (outsl == 18)        out = r1p_1;
      if (outsl == 19)        out = rnpp_1;
      if (outsl == 20)        out = outbypas;
    }
    casc_next_in = out;
    resoncOut();
    out = rout_1;
    temp = out * amp_gain0;         /* Convert back to integer */
    getmax(temp,&sigmx);
    *jwave++ = truncate(temp);   /* Truncate if exceeds 16 bits */
  }
}




/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                   S U B R O U T I N E   G E T H O S T               */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


/*
Get variable parameters from host computer,
initially also get definition of fixed pars */



gethost()

{
  unsigned useed;         /* for Lattice random number gen LG */

/*  Initialize speaker definition */

  if (initsw == 0)
  {
    initsw++;

    FLPhz = (950 * samrate) / 10000;
    BLPhz = (630 * samrate) / 10000;

    ranseed = 1;

    minus_pi_t = -PI / samrate;
    two_pi_t = -2. * minus_pi_t;
    setabc(FLPhz,BLPhz,&rlpa,&rlpb,&rlpc);
    useed = ranseed;
    srand(useed);     /* Init ran # generator */
    resonclr ();      /* LG adds 7/11/87 initializes resonantors for
			 new utterance  */
    nper = 0;         /* LG */
    T0 = 0;           /* LG */
  }


/*
  Read  speech frame definition into temp store
  and move some parameters into active use immediately
  (voice-excited ones are updated pitch synchronously
  to avoid waveform glitches).
*/

  F0hz10= pars[ 0];
  AVdb  = pars[ 1] - 7;
  if (AVdb < 0)
  {
    AVdb = 0;
  }

  F1hz  = pars[ 2];
  B1hz  = pars[ 3];


  F2hz  = pars[ 4];
  B2hz  = pars[ 5];

  F3hz  = pars[ 6];
  B3hz  = pars[ 7];

  F4hz  = pars[ 8];
  B4hz  = pars[ 9];

  F5hz  = pars[10];
  B5hz  = pars[11];


  F6hz  = pars[12];        /* f  of parallel 6th formant */
  B6hz  = pars[13];        /* bw of cascade 6th formant (doesn't exist) */

  FNZhz = pars[14];
  BNZhz = pars[15];

  FNPhz = pars[16];
  BNPhz = pars[17];

  AP    = pars[18];
  amp_aspir = DBtoLIN(AP) * .05;
  Kopen = pars[19];

  Aturb = pars[20];
  TLTdb = pars[21];

  AF    = pars[22];
  amp_frica = DBtoLIN(AF) * 0.25;
  Kskew = pars[23];
  AVpdb = pars[38];
  par_amp_voice = DBtoLIN(AVpdb);

  A1    = pars[24];
  amp_parF1 = DBtoLIN(A1) * 0.4;
  B1phz = pars[25];

  A2    = pars[26];
  amp_parF2 = DBtoLIN(A2) * 0.15;
  B2phz = pars[27];

  A3    = pars[28];
  amp_parF3 = DBtoLIN(A3) * 0.06;
  B3phz = pars[29];

  A4    = pars[30];
  amp_parF4 = DBtoLIN(A4) * 0.04;
  B4phz = pars[31];

  A5    = pars[32];
  amp_parF5 = DBtoLIN(A5) * 0.022;
  B5phz = pars[33];

  A6    = pars[34];
  amp_parF6 = DBtoLIN(A6) * 0.03;
  B6phz = pars[35];

  ANP   = pars[36];
  amp_parFNP = DBtoLIN(ANP) * 0.6;

  AB    = pars[37];
  amp_bypas = DBtoLIN(AB) * 0.05;

  Gain0 = pars[39] - 3;
  if (Gain0 <= 0)
  {
    Gain0 = 57;
  }
  amp_gain0 = DBtoLIN(Gain0);



/*    Set coefficients of variable cascade resonators */

  if (nfcascade >= 8)
  {
    setabc(7500,600,&r8ca,&r8cb,&r8cc);
  }
  if (nfcascade >= 7)
  {
    setabc(6500,500,&r7ca,&r7cb,&r7cc);
  }
  if (nfcascade >= 6)
  {
    setabc(F6hz,B6hz,&r6ca,&r6cb,&r6cc);
  }
  if (nfcascade >= 5)
  {
    setabc(F5hz,B5hz,&r5ca,&r5cb,&r5cc);
  }
  setabc(F4hz,B4hz,&r4ca,&r4cb,&r4cc);
  setabc(F3hz,B3hz,&r3ca,&r3cb,&r3cc);
  setabc(F2hz,B2hz,&r2ca,&r2cb,&r2cc);
  setabc(F1hz,B1hz,&r1ca,&r1cb,&r1cc);

/*    Set coeficients of nasal resonator and zero antiresonator */

  setabc(FNPhz,BNPhz,&rnpca,&rnpcb,&rnpcc);
  setzeroabc(FNZhz,BNZhz,&rnza,&rnzb,&rnzc);

/*    Set coefficients of parallel resonators, and amplitude of outputs */

  setabc(F1hz,B1phz,&r1pa,&r1pb,&r1pc);
  r1pa *= amp_parF1;
  setabc(FNPhz,BNPhz,&rnppa,&rnppb,&rnppc);
  rnppa *= amp_parFNP;
  setabc(F2hz,B2phz,&r2pa,&r2pb,&r2pc);
  r2pa *= amp_parF2;
  setabc(F3hz,B3phz,&r3pa,&r3pb,&r3pc);
  r3pa *= amp_parF3;
  setabc(F4hz,B4phz,&r4pa,&r4pb,&r4pc);
  r4pa *= amp_parF4;
  setabc(F5hz,B5phz,&r5pa,&r5pb,&r5pc);
  r5pa *= amp_parF5;
  setabc(F6hz,B6phz,&r6pa,&r6pb,&r6pc);
  r6pa *= amp_parF6;

/* output low-pass filter */

  setabc((long)0.0,(long)8000.0,&routa,&routb,&routc);

  dispt += nspfr;
  disptcum += nspfr;
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*       S U B R O U T I N E   I M P U L S I V E - S O U R C E         */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

static float doublet[] = { 0., 13000000., -13000000. };

float impulsive_source() {

    if (nper < 3)
    {
      vwave = doublet[nper];
    }
    else
    {
      vwave = 0.;
    }

/*    Low-pass filter the differenciated impulse with a critically-damped
      second-order filter, time constant proportional to Kopen */

    resonglot();
    return(rgl_1);
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*         S U B R O U T I N E   N A T U R A L - S O U R C E           */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*  Vwave is the differentiated glottal flow waveform, there is a weak
    spectral zero around 800 Hz, magic constants a,b reset pitch-synch */

float natural_source()

{

/*    See if glottis open */

  float lgtemp;

  if (nper < nopen)
  {
    a -= b;
    vwave += a;
    lgtemp=vwave * 0.028;
    return(lgtemp);
  }

/*    Glottis closed */
  else
  {
    vwave = 0.;
    return(0.0);
  }
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*  S U B R O U T I N E   P I T C H _ S Y N C _ P A R _ U P D A T E    */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Reset selected parameters pitch-synchronously */

pitch_synch_par_reset()

{
  long zero;

  if (F0hz10 > 0)
  {
    T0 = (40 * samrate) / F0hz10;        /* Period in samp*4 */
    amp_voice = DBtoLIN(AVdb);

/*        Duration of period before amplitude modulation */
    nmod = T0;
    if (AVdb > 0)
    {
      nmod >>= 1;
    }

/*        Breathiness of voicing waveform */

    amp_breth = DBtoLIN(Aturb) * 0.1;

/*        Set open phase of glottal period */
/*        where  40 <= open phase <= 263 */

    nopen = 4 * Kopen;
    if ((glsource == IMPULSIVE) && (nopen > 263))
    {
      nopen = 263;
    }

    if (nopen >= (T0-1))
    {
      nopen = T0 - 2;
      if(quiet_flag == FALSE)
      {
	printf("Warning: glottal open period cannot exceed T0, truncated\n");
      }
    }

    if (nopen < 40)
    {
      nopen = 40;    /* F0 max = 1000 Hz */
      if(quiet_flag == FALSE)
      {
	printf("Warning: minimum glottal open period is 10 samples.\n");
	printf("truncated, nopen = %d\n",nopen);
      }
    }

/*        Reset a & b, which determine shape of "natural" glottal waveform */

    b = B0[nopen-40];
    a = (b * nopen) * .333;

/*        Reset width of "impulsive" glottal pulse */

    temp = samrate / nopen;
    zero = 0;
    setabc(zero,temp,&rgla,&rglb,&rglc);

/*        Make gain at F1 about constant */

    temp1 = nopen *.00833;
    rgla *= temp1 * temp1;

/*        Truncate skewness so as not to exceed duration of closed phase
          of glottal period */

    temp = T0 - nopen;
    if (Kskew > temp)
    {
      if(quiet_flag == FALSE)
      {
	printf("Kskew duration=%d > glottal closed period=%d, truncate\n",
	       Kskew, T0-nopen);
      }
      Kskew = temp;
    }
    if (skew >= 0)
    {
      skew = Kskew;   /* Reset skew to requested Kskew */
    }
    else
    {
      skew = - Kskew;
    }

/*        Add skewness to closed portion of voicing period */

    T0 = T0 + skew;
    skew = - skew;
  }

  else
  {
    T0 = 4;                     /* Default for f0 undefined */
    amp_voice = 0.;
    nmod = T0;
    amp_breth = 0.;
    a = 0.;
    b = 0.;
  }

/*    Reset these pars pitch synchronously or at update rate if f0=0 */

  if ((T0 != 4) || (ns == 0))
  {

/*        Set one-pole low-pass filter that tilts glottal source */

    decay = (.033 * TLTdb);
    if (decay > 0.)
    {
      onemd = 1. - decay;
    }
    else
    {
      onemd = 1.;
    }
  }
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                   S U B R O U T I N E   S E T A B C                 */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*      Convert formant freqencies and bandwidth into
 *      resonator difference equation constants */


setabc(f,bw,acoef,bcoef,ccoef)

long    f;                      /* Frequency of resonator in Hz         */
long    bw;                     /* Bandwidth of resonator in Hz         */
float   *acoef, *bcoef, *ccoef; /* Are output coefficients              */

{
  float r;
  double arg,exp(),cos();

/*    Let r  =  exp(-pi bw t) */

  arg = minus_pi_t * bw;
  r = exp(arg);

/*    Let c  =  -r**2 */

  *ccoef = -(r * r);

/*    Let b = r * 2*cos(2 pi f t) */

  arg = two_pi_t * f;
  *bcoef = r * cos(arg) * 2.;

/*    Let a = 1.0 - b - c */

  *acoef = 1.0 - *bcoef - *ccoef;
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*              S U B R O U T I N E   S E T Z E R O A B C              */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*      Convert formant freqencies and bandwidth into
 *      anti-resonator difference equation constants */

setzeroabc(f,bw,acoef,bcoef,ccoef)

long    f;                      /* Frequency of resonator in Hz         */
long    bw;                     /* Bandwidth of resonator in Hz         */
float   *acoef, *bcoef, *ccoef; /* Are output coefficients              */

{
  float r;
  double arg,exp(),cos();

/*    First compute ordinary resonator coefficients */
/*    Let r  =  exp(-pi bw t) */

  arg = minus_pi_t * bw;
  r = exp(arg);

/*    Let c  =  -r**2 */

  *ccoef = -(r * r);

/*    Let b = r * 2*cos(2 pi f t) */

  arg = two_pi_t * f;
  *bcoef = r * cos(arg) * 2.;

/*    Let a = 1.0 - b - c */

  *acoef = 1. - *bcoef - *ccoef;

/*    Now convert to antiresonator coefficients (a'=1/a, b'=b/a, c'=c/a) */

  *acoef = 1.0 / *acoef;
  *ccoef *= -*acoef;
  *bcoef *= -*acoef;
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*            S U B R O U T I N E   G E N - N O I S E                  */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Random number generator (return a number between -8191 and +8191) */

gen_noise()

{
  long temp;

  temp = rand();
  nrand = (temp >> 17) - 8191;

/*    Tilt down noise spectrum by soft low-pass filter having
 *    a pole near the origin in the z-plane, i.e.
 *    output = input + (0.75 * lastoutput) */

  noise = nrand + (0.75 * nlast);
  nlast = noise;
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                   S U B R O U T I N E   D B t o L I N               */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*      Convert from decibels to a linear scale factor
 */

float DBtoLIN(dB)

long dB;

{

/*    Check limits or argument (can be removed in final product) */

    float lgtemp;
    if (dB < 0)
    {
      if(quiet_flag == FALSE)
      {
	printf("Try to compute amptable[%d]\n", dB);
	return(0);
      }
    }
    lgtemp=amptable[dB] * .001;
    return(lgtemp);
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                   S U B R O U T I N E   R E S O N L P               */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Critically-damped Low-Pass Resonator of Impulsive Glottal Source */

resonglot()

{
  register long temp3,temp4;

  temp4 = rglc * rgl_2;          /*   (ccoef * old2)     */
  rgl_2 = rgl_1;

  temp3 = rglb * rgl_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = rgla * vwave;          /* + (acoef * input)    */
  rgl_1 = temp4 + temp3;
}


/* Low-Pass Downsampling Resonator of Glottal Source */

resonlp()

{
  register long temp3,temp4;

  temp4 = rlpc * rlp_2;          /*   (ccoef * old2)     */
  rlp_2 = rlp_1;

  temp3 = rlpb * rlp_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = rlpa * voice;          /* + (acoef * input)    */
  rlp_1 = temp4 + temp3;
  voice = rlp_1;
}


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                S U B R O U T I N E   R E S O N C                    */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Nasal Antiresonator of Cascade Vocal Tract:
 *  Output = (rnza * input) + (rnzb * oldin1) + (rnzc * oldin2) */

resoncnz()

{
  register long temp3,temp4;

  temp4 = rnzc * rnz_2;          /*   (ccoef * oldin2)   */
  rnz_2 = rnz_1;

  temp3 = rnzb * rnz_1;          /* + (bcoef * oldin1)   */
  temp4 += temp3;

  temp3 = rnza * glotout;            /* + (acoef * input)    */
  rnz_1 = glotout;
  rnzout = temp4 + temp3;
}

/* Nasal Resonator of Cascade Vocal Tract */

resoncnp()

{
  register long temp3,temp4;

  temp4 = rnpcc * rnpc_2;        /*   (ccoef * old2)     */
  rnpc_2 = rnpc_1;

  temp3 = rnpcb * rnpc_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = rnpca * rnzout;          /* + (acoef * input)    */
  rnpc_1 = temp4 + temp3;
}

/* Eighth cascaded Formant */

resonc8()

{
  register long temp3,temp4;

  temp4 = r8cc * r8c_2;          /*   (ccoef * old2)     */
  r8c_2 = r8c_1;

  temp3 = r8cb * r8c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r8ca * casc_next_in;   /* + (acoef * input)    */
  r8c_1 = temp4 + temp3;
}

/* Seventh cascaded Formant */

resonc7()

{
  register long temp3,temp4;

  temp4 = r7cc * r7c_2;          /*   (ccoef * old2)     */
  r7c_2 = r7c_1;

  temp3 = r7cb * r7c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r7ca * casc_next_in;   /* + (acoef * input)    */
  r7c_1 = temp4 + temp3;
}

/* Sixth cascaded Formant */

resonc6()

{
  register long temp3,temp4;

  temp4 = r6cc * r6c_2;          /*   (ccoef * old2)     */
  r6c_2 = r6c_1;

  temp3 = r6cb * r6c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r6ca * casc_next_in;   /* + (acoef * input)    */
  r6c_1 = temp4 + temp3;
}

/* Fifth Formant */

resonc5()

{
  register long temp3,temp4;

  temp4 = r5cc * r5c_2;          /*   (ccoef * old2)     */
  r5c_2 = r5c_1;

  temp3 = r5cb * r5c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r5ca * casc_next_in;   /* + (acoef * input)    */
  r5c_1 = temp4 + temp3;
}

/* Fourth Formant */

resonc4()

{
  register long temp3,temp4;

  temp4 = r4cc * r4c_2;          /*   (ccoef * old2)     */
  r4c_2 = r4c_1;

  temp3 = r4cb * r4c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r4ca * casc_next_in;   /* + (acoef * input)    */
  r4c_1 = temp4 + temp3;
}

/* Third Formant */

resonc3()

{
  register long temp3,temp4;

  temp4 = r3cc * r3c_2;          /*   (ccoef * old2)     */
  r3c_2 = r3c_1;

  temp3 = r3cb * r3c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r3ca * casc_next_in;   /* + (acoef * input)    */
  r3c_1 = temp4 + temp3;
}

/* Second Formant */

resonc2()

{
  register long temp3,temp4;

  temp4 = r2cc * r2c_2;          /*   (ccoef * old2)     */
  r2c_2 = r2c_1;

  temp3 = r2cb * r2c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r2ca * casc_next_in;   /* + (acoef * input)    */
  r2c_1 = temp4 + temp3;
}

/* First Formant of Cascade Vocal Tract */

resonc1()

{
  register long temp3,temp4;

  temp4 = r1cc * r1c_2;          /*   (ccoef * old2)     */
  r1c_2 = r1c_1;

  temp3 = r1cb * r1c_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = r1ca * casc_next_in;   /* + (acoef * input)    */
  r1c_1 = temp4 + temp3;
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                  S U B R O U T I N E   R E S O N P                  */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/*   Output = (acoef * input) + (bcoef * old1) + (ccoef * old2); */

/* Sixth Formant of Parallel Vocal Tract */

reson6p()

{
  register long temp3,temp4;

  temp4 = r6pc * r6p_2;
  r6p_2 = r6p_1;

  temp3 = r6pb * r6p_1;
  temp4 += temp3;

  temp3 = r6pa * sourc;
  r6p_1 = temp4 + temp3;
}

/* Fifth Formant of Parallel Vocal Tract */

reson5p()

{
  register long temp3,temp4;

  temp4 = r5pc * r5p_2;
  r5p_2 = r5p_1;

  temp3 = r5pb * r5p_1;
  temp4 += temp3;

  temp3 = r5pa * sourc;
  r5p_1 = temp4 + temp3;
}

/* Fourth Formant of Parallel Vocal Tract */

reson4p()

{
  register long temp3,temp4;

  temp4 = r4pc * r4p_2;
  r4p_2 = r4p_1;

  temp3 = r4pb * r4p_1;
  temp4 += temp3;

  temp3 = r4pa * sourc;
  r4p_1 = temp4 + temp3;
}

/* Third Formant of Parallel Vocal Tract */

reson3p()

{
  register long temp3,temp4;

  temp4 = r3pc * r3p_2;
  r3p_2 = r3p_1;

  temp3 = r3pb * r3p_1;
  temp4 += temp3;

  temp3 = r3pa * sourc;
  r3p_1 = temp4 + temp3;
}

/* Second Formant of Parallel Vocal Tract */

reson2p()

{
  register long temp3,temp4;

  temp4 = r2pc * r2p_2;
  r2p_2 = r2p_1;

  temp3 = r2pb * r2p_1;
  temp4 += temp3;

  temp3 = r2pa * sourc;
  r2p_1 = temp4 + temp3;
}


/* First Formant of Parallel Vocal Tract */

reson1p()

{
  register long temp3,temp4;

  temp4 = r1pc * r1p_2;
  r1p_2 = r1p_1;

  temp3 = r1pb * r1p_1;
  temp4 += temp3;

  temp3 = r1pa * sourc;
  r1p_1 = temp4 + temp3;
}

/* Nasal Formant of Parallel Vocal Tract */

resonnpp()

{
  register long temp3,temp4;

  temp4 = rnppc * rnpp_2;
  rnpp_2 = rnpp_1;

  temp3 = rnppb * rnpp_1;
  temp4 += temp3;

  temp3 = rnppa * sourc;
  rnpp_1 = temp4 + temp3;
}




/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*               S U B R O U T I N E   N O - R A D - C H A R           */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#define ACOEF           0.005
#define BCOEF           (1.0 - ACOEF)   /* Slight decay to remove dc */

no_rad_char(in)

float in;

{
  static float lastin;

  out = (ACOEF * in) + (BCOEF * lastin);
  lastin = out;
  out = -100. * out;      /* Scale up to make visible */
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                  S U B R O U T I N E   G E T M A X                  */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Find absolute maximum of arg1 & arg2, save in arg2 */

getmax(arg1,arg2)

long arg1;
long *arg2;

{
  if (arg1 < 0)
  {
    arg1 = - arg1;
  }

  if (arg1 > *arg2)
  {
    *arg2 = arg1;
  }
}



/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                     */
/*                 S U B R O U T I N E   T R U N C A T E               */
/*                                                                     */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* Truncate arg to fit into 16-bit word */

long truncate(arg)

long arg;

{
  if (arg < -32768)
  {
    overload_warning(-arg);
    arg = -32768;
  }
  if (arg >  32767)
  {
    overload_warning(arg);
    arg =  32767;
  }
  return(arg);
}

float dBconvert(arg)

long arg;

{
  double x,log10();
  float db;

  x = arg / 32767.;
  x = log10(x);
  db = 20. * x;
  return(db);
}


overload_warning(arg) long arg;

{
  if (warnsw == 0)
  {
    warnsw++;
    if(quiet_flag==FALSE)
    {
      printf("\n* * * WARNING: ");
      printf(" Signal at output of synthesizer (+%3.1f dB) exceeds 0 dB\n",
	     dBconvert(arg));
      printf("    at synthesis parameter frame = %d\n",disptcum/nspfr);
    }
  }
}

resonclr ()

{
  rnpp_1 = 0;  /* Last output sample from parallel nasal pole  */
  rnpp_2 = 0;  /* Second-previous output sample                */

  r1p_1  = 0;  /* Last output sample from parallel 1st formant */
  r1p_2  = 0;  /* Second-previous output sample                */

  r2p_1  = 0;  /* Last output sample from parallel 2nd formant */
  r2p_2  = 0;  /* Second-previous output sample                */

  r3p_1  = 0;  /* Last output sample from parallel 3rd formant */
  r3p_2  = 0;  /* Second-previous output sample                */

  r4p_1  = 0;  /* Last output sample from parallel 4th formant */
  r4p_2  = 0;  /* Second-previous output sample                */

  r5p_1  = 0;  /* Last output sample from parallel 5th formant */
  r5p_2  = 0;  /* Second-previous output sample                */

  r6p_1  = 0;  /* Last output sample from parallel 6th formant */
  r6p_2  = 0;  /* Second-previous output sample                */

  r1c_1  = 0;  /* Last output sample from cascade 1st formant  */
  r1c_2  = 0;  /* Second-previous output sample                */

  r2c_1  = 0;  /* Last output sample from cascade 2nd formant  */
  r2c_2  = 0;  /* Second-previous output sample                */

  r3c_1  = 0;  /* Last output sample from cascade 3rd formant  */
  r3c_2  = 0;  /* Second-previous output sample                */

  r4c_1  = 0;  /* Last output sample from cascade 4th formant  */
  r4c_2  = 0;  /* Second-previous output sample                */

  r5c_1  = 0;  /* Last output sample from cascade 5th formant  */
  r5c_2  = 0;  /* Second-previous output sample                */

  r6c_1  = 0;  /* Last output sample from cascade 6th formant  */
  r6c_2  = 0;  /* Second-previous output sample                */

  r7c_1 = 0;
  r7c_2 = 0;

  r8c_1 = 0;
  r8c_2 = 0;

  rnpc_1 = 0;  /* Last output sample from cascade nasal pole   */
  rnpc_2 = 0;  /* Second-previous output sample                */

  rnz_1  = 0;  /* Last output sample from cascade nasal zero   */
  rnz_2  = 0;  /* Second-previous output sample                */

  rgl_1  = 0;  /* Last output crit-damped glot low-pass filter */
  rgl_2  = 0;  /* Second-previous output sample                */

  rlp_1  = 0;  /* Last output from downsamp low-pass filter    */
  rlp_2  = 0;  /* Second-previous output sample                */

  vlast  = 0;  /* Previous output of voice                     */
  nlast  = 0;  /* Previous output of random number generator   */
  glotlast = 0; /* Previous value of glotout                   */
}

resoncOut()

{
  register long temp3,temp4;

  temp4 = routc * rout_2;          /*   (ccoef * old2)     */
  rout_2 = rout_1;

  temp3 = routb * rout_1;          /* + (bcoef * old1)     */
  temp4 += temp3;

  temp3 = routa * casc_next_in;   /* + (acoef * input)    */
  rout_1 = temp4 + temp3;
}




