#include "config.h"

/* $Id: nsynth.c,v 1.13 1994/11/08 13:30:50 a904209 Exp a904209 $
 */

/* 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 "useconfig.h"
#include <stdio.h>
#include <math.h>
#include "proto.h"
#include "nsynth.h"
#ifndef PI
#ifndef M_PI                      /* <math.h> */
#define PI               3.1415927
#else /* M_PI */
#define PI               M_PI
#endif /* M_PI */
#endif

#ifdef __STDC__
#define ONE 1.0F
#else
#define ONE 1.0
#endif

typedef struct
{
	const char *name;
	float a;
	float b;
	float c;
	float p1;
	float p2;
}
resonator_t, *resonator_ptr;

/* Various global variables */

int time_count = 0;
static int warnsw;                    /* JPI added for f0 flutter */
float global_gain = 1;            /* global gain parameter for voice */

/* COUNTERS */

static int32_t nper;              /* Current loc in voicing period   40000 samp/s */

/* COUNTER LIMITS */

static int32_t T0;                /* Fundamental period in output samples times 4 */
static int32_t nopen;             /* Number of samples in open phase of period  */
static int32_t nmod;              /* Position in period to begin noise amp. modul */

/* Variables that have to stick around for awhile, and thus locals
   are not appropriate
 */

/* Incoming parameter Variables which need to be updated synchronously  */

static int32_t F0hz10;            /* Voicing fund freq in Hz  */
static int32_t AVdb;              /* Amp of voicing in dB,    0 to   70  */
static int32_t Kskew;             /* Skewness of alternate periods,0 to   40  */

/* Various amplitude variables used in main loop */

static float amp_voice;           /* AVdb converted to linear gain  */
static float amp_bypas;           /* AB converted to linear gain  */
static float par_amp_voice;       /* AVpdb converted to linear gain  */
static float amp_aspir;           /* AP converted to linear gain  */
static float amp_frica;           /* AF converted to linear gain  */
static float amp_breth;           /* ATURB converted to linear gain  */

/* State variables of sound sources */

static int32_t skew;              /* Alternating jitter, in half-period units  */

static float natglot_a;           /* Makes waveshape of glottal pulse when open  */
static float natglot_b;           /* Makes waveshape of glottal pulse when open  */
static float vwave;               /* Ditto, but before multiplication by AVdb  */
static float vlast;               /* Previous output of voice  */
static float nlast;               /* Previous output of random number generator  */
static float glotlast;            /* Previous value of glotout  */
static float decay;               /* TLTdb converted to exponential time const  */
static float onemd;               /* in voicing one-pole low-pass filter  */
static float minus_pi_t;          /* func. of sample rate */
static float two_pi_t;            /* func. of sample rate */


/* INTERNAL MEMORY FOR DIGITAL RESONATORS AND ANTIRESONATOR  */

static resonator_t rnpp =
{"parallel nasal pole"};
static resonator_t r1p =
{"parallel 1st formant"};
static resonator_t r2p =
{"parallel 2nd formant"};
static resonator_t r3p =
{"parallel 3rd formant"};
static resonator_t r4p =
{"parallel 4th formant"};
static resonator_t r5p =
{"parallel 5th formant"};
static resonator_t r6p =
{"parallel 6th formant"};
static resonator_t r1c =
{"cascade 1st formant"};
static resonator_t r2c =
{"cascade 2nd formant"};
static resonator_t r3c =
{"cascade 3rd formant"};
static resonator_t r4c =
{"cascade 4th formant"};
static resonator_t r5c =
{"cascade 5th formant"};
static resonator_t r6c =
{"cascade 6th formant"};
static resonator_t r7c =
{"cascade 7th formant"};
static resonator_t r8c =
{"cascade 8th formant"};
static resonator_t rnpc =
{"cascade nasal pole"};
static resonator_t rnz =
{"cascade nasal zero"};
static resonator_t rgl =
{"crit-damped glot low-pass filter"};
static resonator_t rlp =
{"downsamp low-pass filter"};
static resonator_t rout =
{"output low-pass"};

/*
 * Constant natglot[] controls shape of glottal pulse as a function
 * of desired duration of open phase N0
 * (Note that N0 is specified in terms of 40,000 samples/sec of speech)
 *
 *    Assume voicing waveform V(t) has form: k1 t**2 - k2 t**3
 *
 *    If the radiation characterivative, a temporal derivative
 *      is folded in, and we go from continuous time to discrete
 *      integers n:  dV/dt = vwave[n]
 *                         = sum over i=1,2,...,n of { a - (i * b) }
 *                         = a n  -  b/2 n**2
 *
 *      where the  constants a and b control the detailed shape
 *      and amplitude of the voicing waveform over the open
 *      potion of the voicing cycle "nopen".
 *
 *    Let integral of dV/dt have no net dc flow --> a = (b * nopen) / 3
 *
 *    Let maximum of dUg(n)/dn be constant --> b = gain / (nopen * nopen)
 *      meaning as nopen gets bigger, V has bigger peak proportional to n
 *
 *    Thus, to generate the table below for 40 <= nopen <= 263:
 *
 *      natglot[nopen - 40] = 1920000 / (nopen * nopen)
 */
static const short natglot[224] =
{
	1200, 1142, 1088, 1038, 991, 948, 907, 869, 833, 799,
	768, 738, 710, 683, 658, 634, 612, 590, 570, 551,
	533, 515, 499, 483, 468, 454, 440, 427, 415, 403,
	391, 380, 370, 360, 350, 341, 332, 323, 315, 307,
	300, 292, 285, 278, 272, 265, 259, 253, 247, 242,
	237, 231, 226, 221, 217, 212, 208, 204, 199, 195,
	192, 188, 184, 180, 177, 174, 170, 167, 164, 161,
	158, 155, 153, 150, 147, 145, 142, 140, 137, 135,
	133, 131, 128, 126, 124, 122, 120, 119, 117, 115,
	113, 111, 110, 108, 106, 105, 103, 102, 100, 99,
	97, 96, 95, 93, 92, 91, 90, 88, 87, 86,
	85, 84, 83, 82, 80, 79, 78, 77, 76, 75,
	75, 74, 73, 72, 71, 70, 69, 68, 68, 67,
	66, 65, 64, 64, 63, 62, 61, 61, 60, 59,
	59, 58, 57, 57, 56, 56, 55, 55, 54, 54,
	53, 53, 52, 52, 51, 51, 50, 50, 49, 49,
	48, 48, 47, 47, 46, 46, 45, 45, 44, 44,
	43, 43, 42, 42, 41, 41, 41, 41, 40, 40,
	39, 39, 38, 38, 38, 38, 37, 37, 36, 36,
	36, 36, 35, 35, 35, 35, 34, 34, 33, 33,
	33, 33, 32, 32, 32, 32, 31, 31, 31, 31,
	30, 30, 30, 30, 29, 29, 29, 29, 28, 28,
	28, 28, 27, 27
};

/*
 * Convertion table, db to linear, 87 dB --> 32767
 *                                 86 dB --> 29491 (1 dB down = 0.5**1/6)
 *                                 ...
 *                                 81 dB --> 16384 (6 dB down = 0.5)
 *                                 ...
 *                                  0 dB -->     0
 *
 * The just noticeable difference for a change in intensity of a vowel
 *   is approximately 1 dB.  Thus all amplitudes are quantized to 1 dB
 *   steps.
 */

static const float amptable[88] =
{
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 6.0, 7.0,
	8.0, 9.0, 10.0, 11.0, 13.0,
	14.0, 16.0, 18.0, 20.0, 22.0,
	25.0, 28.0, 32.0, 35.0, 40.0,
	45.0, 51.0, 57.0, 64.0, 71.0,
	80.0, 90.0, 101.0, 114.0, 128.0,
	142.0, 159.0, 179.0, 202.0, 227.0,
	256.0, 284.0, 318.0, 359.0, 405.0,
	455.0, 512.0, 568.0, 638.0, 719.0,
	811.0, 911.0, 1024.0, 1137.0, 1276.0,
	1438.0, 1622.0, 1823.0, 2048.0, 2273.0,
	2552.0, 2875.0, 3244.0, 3645.0, 4096.0,
	4547.0, 5104.0, 5751.0, 6488.0, 7291.0,
	8192.0, 9093.0, 10207.0, 11502.0, 12976.0,
	14582.0, 16384.0, 18350.0, 20644.0, 23429.0,
	26214.0, 29491.0, 32767.0
};

#if 0
static const char *par_name[] =
{
	"F0hz10",
	"AVdb",
	"F1hz", "B1hz",
	"F2hz", "B2hz",
	"F3hz", "B3hz",
	"F4hz", "B4hz",
	"F5hz", "B5hz",
	"F6hz", "B6hz",
	"FNZhz", "BNZhz",
	"FNPhz", "BNPhz",
	"AP",
	"Kopen",
	"Aturb",
	"TLTdb",
	"AF",
	"Kskew",
	"A1", "B1phz",
	"A2", "B2phz",
	"A3", "B3phz",
	"A4", "B4phz",
	"A5", "B5phz",
	"A6", "B6phz",
	"ANP",
	"AB",
	"AVpdb",
	"Gain0"
};
#endif

static void flutter PROTO((klatt_global_ptr globals, klatt_frame_ptr pars));
static float resonator PROTO((resonator_ptr r, Float input));
static float antiresonator PROTO((resonator_ptr r, Float input));
static float impulsive_source PROTO((int32_t nper));
static float natural_source PROTO((int32_t nper));
static void setabc PROTO((int32_t f, int32_t bw, resonator_ptr rp));
static void setabcg PROTO((int32_t f, int32_t bw, resonator_ptr rp, Float gain));
static void setzeroabc PROTO((int32_t f, int32_t bw, resonator_ptr rp));
static float DBtoLIN PROTO((klatt_global_ptr globals, int32_t dB));
static void overload_warning PROTO((klatt_global_ptr globals, int32_t arg));
static short clip PROTO((klatt_global_ptr globals, Float input));
static void pitch_synch_par_reset PROTO((klatt_global_ptr globals,
                                         klatt_frame_ptr frame, int32_t ns));
static void frame_init PROTO((klatt_global_ptr globals, klatt_frame_ptr frame));

#ifdef _MSC_VER /* disable some silent conversion warnings */
#pragma warning(disable:4244)
#endif

#if 0
void show_parms PROTO((klatt_global_ptr globals, int *pars));
void show_parms( klatt_global_ptr globals,int *pars )
{
	int i;
	static int names;
	if ((names++ % 64) == 0)
	{
		for (i = 0; i < NPAR; i++)
			printf( "%s ", par_name[i] );
		printf( "\n" );
	}
	for (i = 0; i < NPAR; i++)
	{
		printf( "%*d ", (int) strlen( par_name[i] ), pars[i] );
	}
	printf( "\n" );
}
#endif

/*
   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.
 */
static void flutter( klatt_global_ptr globals,klatt_frame_ptr pars )
{
	int32_t original_f0 = pars->F0hz10 / 10;
	double fla = (double) globals->f0_flutter / 50;
	double flb = (double) original_f0 / 100;
	double flc = sin( 2 * PI * 12.7 * time_count );
	double fld = sin( 2 * PI * 7.1 * time_count );
	double fle = sin( 2 * PI * 4.7 * time_count );
	double delta_f0 = fla * flb * (flc + fld + fle) * 10;
	F0hz10 += (int32_t) delta_f0;
}

static float impulsive_source( int32_t nper )
{
	static float doublet[] =
	{0., 13000000., -13000000.};
	if (nper < 3)
	{
		vwave = doublet[nper];
	}
	else
	{
		vwave = 0.0;
	}
	/* Low-pass filter the differenciated impulse with a critically-damped
	   second-order filter, time constant proportional to Kopen */
	return resonator( &rgl, vwave );
}


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

static float natural_source( int32_t nper )
{
	float lgtemp;
	/* See if glottis open */
	if (nper < nopen)
	{
		natglot_a -= natglot_b;
		vwave += natglot_a;
		lgtemp = vwave * 0.028f;   /* function of samp_rate ? */
		return (lgtemp);
	}
	else
	{
		/* Glottis closed */
		vwave = 0.0;
		return (0.0);
	}
}

/*----------------------------------------------------------------------------*/
/* Convert formant freqencies and bandwidth into
   resonator difference equation coefficents
 */
static void setabc( int32_t f,int32_t bw,resonator_ptr rp )
//int32_t f;                        /* Frequency of resonator in Hz  */
//int32_t bw;                       /* Bandwidth of resonator in Hz  */
//resonator_ptr rp;                 /* Are output coefficients  */
{
	double arg = minus_pi_t * bw;
	float r = exp( arg );         /* Let r  =  exp(-pi bw t) */
	rp->c = -(r * r);             /* Let c  =  -r**2 */
	arg = two_pi_t * f;
	rp->b = r * cos( arg ) * 2.0f; /* Let b = r * 2*cos(2 pi f t) */
	rp->a = 1.0f - rp->b - rp->c;  /* Let a = 1.0 - b - c */
}

/* Convienience function for setting parallel resonators with gain */
static void setabcg( int32_t f, int32_t bw, resonator_ptr rp, Float gain )
/*int32_t f;                        Frequency of resonator in Hz  */
/*int32_t bw;                       Bandwidth of resonator in Hz  */
/*resonator_ptr rp;                 Are output coefficients  */
/*Float gain;*/
{
	setabc( f, bw, rp );
	rp->a *= gain;
}

/* Convert formant freqencies and bandwidth into
 *      anti-resonator difference equation constants
 */
static void setzeroabc( int32_t f,int32_t bw,resonator_ptr rp )
//int32_t f;                        /* Frequency of resonator in Hz  */
//int32_t bw;                       /* Bandwidth of resonator in Hz  */
//resonator_ptr rp;                 /* Are output coefficients  */
{
	setabc( f, bw, rp );          /* First compute ordinary resonator coefficients */
	/* Now convert to antiresonator coefficients */
	rp->a = 1.0f / rp->a;          /* a'=  1/a */
	rp->b *= -rp->a;              /* b'= -b/a */
	rp->c *= -rp->a;              /* c'= -c/a */
}

/*----------------------------------------------------------------------------*/


/* Convert from decibels to a linear scale factor */
static float DBtoLIN( klatt_global_ptr globals,int32_t dB )
{
	/* Check limits or argument (can be removed in final product) */
	if (dB < 0)
		dB = 0;
	else if (dB >= 88)
	{
		#ifdef NSYNTH_VERBOSE
		if (!globals->quiet_flag)
			printf( "Try to compute amptable[%d]\n", dB );
		#endif
		dB = 87;
	}
	return amptable[dB] * 0.001f;
}

/* WHAT WERE THESE FOR ? */
#define ACOEF           0.005
#define BCOEF           (1.0 - ACOEF)   /* Slight decay to remove dc */

#ifdef NSYNTH_VERBOSE
static float dBconvert PROTO((int32_t arg));
static float dBconvert( int32_t arg )
{
	return 20.0 * log10((double) arg / 32767.0 );
}

static void overload_warning( klatt_global_ptr globals,int32_t arg )
{
	if (warnsw == 0)
	{
		warnsw++;
		if (!globals->quiet_flag)
		{
			printf( "\n* * * WARNING: " );
			printf( " Signal at output of synthesizer (+%3.1f dB) exceeds 0 dB\n",
			        dBconvert( arg ));
		}
	}
}
#else
static void overload_warning(klatt_global_ptr globals,int32_t arg) { }
#endif

/* Reset selected parameters pitch-synchronously */

static void pitch_synch_par_reset( klatt_global_ptr globals,klatt_frame_ptr frame,int32_t ns )
{
	int32_t temp;
	float temp1;
	if (F0hz10 > 0)
	{
		T0 = (40 * globals->samrate) / F0hz10;
		/* Period in samp*4 */
		amp_voice = DBtoLIN( globals, AVdb );

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

		/* Breathiness of voicing waveform */

		amp_breth = DBtoLIN( globals, frame->Aturb ) * 0.1f;

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

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

		if (nopen >= (T0 - 1))
		{
			nopen = T0 - 2;
			#ifdef NSYNTH_VERBOSE
			if (!globals->quiet_flag)
				printf( "Warning: glottal open period cannot exceed T0, truncated\n" );
			#endif
		}

		if (nopen < 40)
		{
			nopen = 40;           /* F0 max = 1000 Hz */
			#ifdef NSYNTH_VERBOSE
			if (!globals->quiet_flag)
			{
				printf( "Warning: minimum glottal open period is 10 samples.\n" );
				printf( "truncated, nopen = %d\n", nopen );
			}
			#endif
		}

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

		natglot_b = natglot[nopen - 40];
		natglot_a = (natglot_b * nopen) * .333f;

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

		temp = globals->samrate / nopen;
		setabc( 0L, temp, &rgl );

		/* Make gain at F1 about constant */

		temp1 = nopen * .00833f;
		rgl.a *= (temp1 * temp1);

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

		temp = T0 - nopen;
		if (Kskew > temp)
		{
			#ifdef NSYNTH_VERBOSE
			if (!globals->quiet_flag)
				printf( "Kskew duration=%d > glottal closed period=%d, truncate\n",
				        Kskew, T0 - nopen );
			#endif
			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.0;
		nmod = T0;
		amp_breth = 0.0;
		natglot_a = 0.0;
		natglot_b = 0.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 = (0.033f * frame->TLTdb); /* Function of samp_rate ? */
		if (decay > 0.0)
		{
			onemd = 1.0f - decay;
		}
		else
		{
			onemd = 1.0;
		}
	}
}

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

static void frame_init( klatt_global_ptr globals,klatt_frame_ptr frame )
{
	int32_t Gain0;                /* Overall gain, 60 dB is unity  0 to   60  */
	float amp_parF1;              /* A1 converted to linear gain  */
	float amp_parFN;              /* ANP converted to linear gain  */
	float amp_parF2;              /* A2 converted to linear gain  */
	float amp_parF3;              /* A3 converted to linear gain  */
	float amp_parF4;              /* A4 converted to linear gain  */
	float amp_parF5;              /* A5 converted to linear gain  */
	float amp_parF6;              /* A6 converted to linear gain  */

#if 0
	show_parms((int *) frame );
#endif

	/*
	   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 = frame->F0hz10;
	AVdb = frame->AVdb - 7;
	if (AVdb < 0)
		AVdb = 0;

	amp_aspir = DBtoLIN( globals, frame->ASP ) * .05f;
	amp_frica = DBtoLIN( globals, frame->AF ) * 0.25f;

	Kskew = frame->Kskew;
	par_amp_voice = DBtoLIN( globals, frame->AVpdb );

	/* Fudge factors (which comprehend affects of formants on each other?)
	   with these in place ALL_PARALLEL should sound as close as
	   possible to CASCADE_PARALLEL.
	   Possible problem feeding in Holmes's amplitudes given this.
	 */
	amp_parF1 = DBtoLIN( globals, frame->A1 ) * 0.4f; /* -7.96 dB */
	amp_parF2 = DBtoLIN( globals, frame->A2 ) * 0.15f; /* -16.5 dB */
	amp_parF3 = DBtoLIN( globals, frame->A3 ) * 0.06f; /* -24.4 dB */
	amp_parF4 = DBtoLIN( globals, frame->A4 ) * 0.04f; /* -28.0 dB */
	amp_parF5 = DBtoLIN( globals, frame->A5 ) * 0.022f; /* -33.2 dB */
	amp_parF6 = DBtoLIN( globals, frame->A6 ) * 0.03f; /* -30.5 dB */
	amp_parFN = DBtoLIN( globals, frame->ANP ) * 0.6f; /* -4.44 dB */
	amp_bypas = DBtoLIN( globals, frame->AB ) * 0.05f; /* -26.0 db */

	if (globals->nfcascade >= 8)
	{
		/* Inside Nyquist rate ? */
		if (globals->samrate >= 16000)
			setabc( 7500, 600, &r8c );
		else
			globals->nfcascade = 6;
	}

	if (globals->nfcascade >= 7)
	{
		/* Inside Nyquist rate ? */
		if (globals->samrate >= 16000)
			setabc( 6500, 500, &r7c );
		else
			globals->nfcascade = 6;
	}

	/* Set coefficients of variable cascade resonators */

	if (globals->nfcascade >= 6)
		setabc( frame->F6hz, frame->B6hz, &r6c );

	if (globals->nfcascade >= 5)
		setabc( frame->F5hz, frame->B5hz, &r5c );

	setabc( frame->F4hz, frame->B4hz, &r4c );
	setabc( frame->F3hz, frame->B3hz, &r3c );
	setabc( frame->F2hz, frame->B2hz, &r2c );
	setabc( frame->F1hz, frame->B1hz, &r1c );

	/* Set coeficients of nasal resonator and zero antiresonator */
	setabc( frame->FNPhz, frame->BNPhz, &rnpc );
	setzeroabc( frame->FNZhz, frame->BNZhz, &rnz );

	/* Set coefficients of parallel resonators, and amplitude of outputs */
	setabcg( frame->F1hz, frame->B1phz, &r1p, amp_parF1 );
	setabcg( frame->FNPhz, frame->BNPhz, &rnpp, amp_parFN );
	setabcg( frame->F2hz, frame->B2phz, &r2p, amp_parF2 );
	setabcg( frame->F3hz, frame->B3phz, &r3p, amp_parF3 );
	setabcg( frame->F4hz, frame->B4phz, &r4p, amp_parF4 );
	setabcg( frame->F5hz, frame->B5phz, &r5p, amp_parF5 );
	setabcg( frame->F6hz, frame->B6phz, &r6p, amp_parF6 );


	/* fold overall gain into output resonator */
	Gain0 = frame->Gain0 - 3;
	if (Gain0 <= 0)
		Gain0 = 57;
	/* output low-pass filter - resonator with freq 0 and BW = globals->samrate
	   Thus 3db point is globals->samrate/2 i.e. Nyquist limit.
	   Only 3db down seems rather mild...
	 */
	setabcg( 0L, (int32_t) globals->samrate, &rout, DBtoLIN( globals, Gain0 ));
}

static short clip( klatt_global_ptr globals, Float input )
{
	int32_t temp = input;
	/* clip on boundaries of 16-bit word */
	if (temp < -32767)
	{
		overload_warning( globals, -temp );
		temp = -32767;
	}
	else if (temp > 32767)
	{
		overload_warning( globals, temp );
		temp = 32767;
	}
	return (temp);
}

/* Generic resonator function */
static float resonator( resonator_ptr r, Float input )
{
	register float x = r->a * input + r->b * r->p1 + r->c * r->p2;
	r->p2 = r->p1;
	r->p1 = x;
	return x;
}

/* Generic anti-resonator function
   Same as resonator except that a,b,c need to be set with setzeroabc()
   and we save inputs in p1/p2 rather than outputs.
   There is currently only one of these - "rnz"
 */
/*  Output = (rnz.a * input) + (rnz.b * oldin1) + (rnz.c * oldin2) */

static float antiresonator( resonator_ptr r, Float input )
{
	register float x = r->a * input + r->b * r->p1 + r->c * r->p2;
	r->p2 = r->p1;
	r->p1 = input;
	return x;
}

/*
   function PARWAV

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

void parwave( klatt_global_ptr globals,klatt_frame_ptr frame,short int * jwave )
{
	int32_t ns;
	float out = 0.0;
	/* Output of cascade branch, also final output  */

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

	frame_init( globals, frame );

	if (globals->f0_flutter != 0)
	{
		time_count++;             /* used for f0 flutter */
		flutter( globals, frame ); /* add f0 flutter */
	}

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

	for (ns = 0; ns < globals->nspfr; ns++)
	{
		static uint32_t seed = 5; /* fixed staring value */
		float noise;
		int n4;
		float sourc;              /* Sound source if all-parallel config used  */
		float glotout;            /* Output of glottal sound source  */
		float par_glotout;        /* Output of parallelglottal sound sourc  */
		float voice;              /* Current sample of voicing waveform  */
		float frics;              /* Frication sound source  */
		float aspiration;         /* Aspiration sound source  */
		int32_t nrand;            /* Varible used by random number generator  */

		/* Our own code like rand(), but portable
		   whole upper 31 bits of seed random
		   assumes 32-bit unsigned arithmetic
		 */
		seed = seed * 1664525 + 1;

		/* Shift top bits of seed up to top of long then back down to LS 14 bits */
		/* Assumes 8 bits per sizeof unit i.e. a "byte" */
		nrand = ((int32_t) seed) >> (8 * sizeof(int32_t) - 14);

		/* 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.75f * nlast); /* Function of samp_rate ? */
		nlast = 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 */
		sourc = 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++)
		{
			if (globals->glsource == IMPULSIVE)
			{
				/* Use impulsive glottal source */
				voice = impulsive_source( nper );
			}
			else
			{
				/* Or use a more-natural-shaped source waveform with excitation
				   occurring both upon opening and upon closure, stronest at closure */
				voice = natural_source( nper );
			}

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

			/* Low-pass filter voicing waveform before downsampling from 4*globals->samrate */
			/* to globals->samrate samples/sec.  Resonator f=.09*globals->samrate, bw=.06*globals->samrate  */

			voice = resonator( &rlp, voice ); /* in=voice, out=voice */

			/* Increment counter that keeps track of 4*globals->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;

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

		par_glotout = glotout;

		if (globals->synthesis_model != ALL_PARALLEL)
		{
			/* Cascade vocal tract, excited by laryngeal sources.
			   Nasal antiresonator, then formants FNP, F5, F4, F3, F2, F1
			 */
			float rnzout = antiresonator( &rnz, glotout ); /* Output of cascade nazal zero resonator  */
			float casc_next_in = resonator( &rnpc, rnzout ); /* in=rnzout, out=rnpc.p1 */

			/* Recoded from sequence of if's to use C's fall through switch
			   semantics. May allow compiler to optimize
			 */
			switch (globals->nfcascade)
			{
			case 8:
				casc_next_in = resonator( &r8c, casc_next_in ); /* Do not use unless samrat = 16000 */
			case 7:
				casc_next_in = resonator( &r7c, casc_next_in ); /* Do not use unless samrat = 16000 */
			case 6:
				casc_next_in = resonator( &r6c, casc_next_in ); /* Do not use unless long vocal tract or samrat increased */
			case 5:
				casc_next_in = resonator( &r5c, casc_next_in );
			case 4:
				casc_next_in = resonator( &r4c, casc_next_in );
			case 3:
				casc_next_in = resonator( &r3c, casc_next_in );
			case 2:
				casc_next_in = resonator( &r2c, casc_next_in );
			case 1:
				out = resonator( &r1c, casc_next_in );
				break;
			default:
				out = 0.0;
			}
#if 0
			/* Excite parallel F1 and FNP by voicing waveform */
			/* Source is voicing plus aspiration */
			/* Add in phase, boost lows for nasalized */
			out += (resonator( &rnpp, par_glotout ) + resonator( &r1p, par_glotout ));
#endif
		}
		else
		{
			/* Is ALL_PARALLEL */
			/* NIS - rsynth "hack"
			   As Holmes' scheme is weak at nasals and (physically) nasal cavity
			   is "back near glottis" feed glottal source through nasal resonators
			   Don't think this is quite right, but improves things a bit
			 */
			par_glotout = antiresonator( &rnz, par_glotout );
			par_glotout = resonator( &rnpc, par_glotout );
			/* And just use r1p NOT rnpp */
			out = resonator( &r1p, par_glotout );
			/* Sound sourc for other parallel resonators is frication
			   plus first difference of voicing waveform.
			 */
			sourc += (par_glotout - glotlast);
			glotlast = par_glotout;
		}

		/* Standard parallel vocal tract
		   Formants F6,F5,F4,F3,F2, outputs added with alternating sign
		 */
		out = resonator( &r6p, sourc ) - out;
		out = resonator( &r5p, sourc ) - out;
		out = resonator( &r4p, sourc ) - out;
		out = resonator( &r3p, sourc ) - out;
		out = resonator( &r2p, sourc ) - out;

		out = amp_bypas * sourc - out;

		out = resonator( &rout, out );
		out *= global_gain;
		*jwave++ = clip( globals, out ); /* Convert back to integer */
	}
}

void parwave_init( klatt_global_ptr globals )
{
	int32_t FLPhz = (950 * globals->samrate) / 10000;
	int32_t BLPhz = (630 * globals->samrate) / 10000;

	minus_pi_t = -PI / globals->samrate;
	two_pi_t = -2.0f * minus_pi_t;

	setabc( FLPhz, BLPhz, &rlp );
	nper = 0;                     /* LG */
	T0 = 0;                       /* LG */

	rnpp.p1 = 0;                  /* parallel nasal pole  */
	rnpp.p2 = 0;

	r1p.p1 = 0;                   /* parallel 1st formant */
	r1p.p2 = 0;

	r2p.p1 = 0;                   /* parallel 2nd formant */
	r2p.p2 = 0;

	r3p.p1 = 0;                   /* parallel 3rd formant */
	r3p.p2 = 0;

	r4p.p1 = 0;                   /* parallel 4th formant */
	r4p.p2 = 0;

	r5p.p1 = 0;                   /* parallel 5th formant */
	r5p.p2 = 0;

	r6p.p1 = 0;                   /* parallel 6th formant */
	r6p.p2 = 0;

	r1c.p1 = 0;                   /* cascade 1st formant  */
	r1c.p2 = 0;

	r2c.p1 = 0;                   /* cascade 2nd formant  */
	r2c.p2 = 0;

	r3c.p1 = 0;                   /* cascade 3rd formant  */
	r3c.p2 = 0;

	r4c.p1 = 0;                   /* cascade 4th formant  */
	r4c.p2 = 0;

	r5c.p1 = 0;                   /* cascade 5th formant  */
	r5c.p2 = 0;

	r6c.p1 = 0;                   /* cascade 6th formant  */
	r6c.p2 = 0;

	r7c.p1 = 0;
	r7c.p2 = 0;

	r8c.p1 = 0;
	r8c.p2 = 0;

	rnpc.p1 = 0;                  /* cascade nasal pole  */
	rnpc.p2 = 0;

	rnz.p1 = 0;                   /* cascade nasal zero  */
	rnz.p2 = 0;

	rgl.p1 = 0;                   /* crit-damped glot low-pass filter */
	rgl.p2 = 0;

	rlp.p1 = 0;                   /* downsamp low-pass filter  */
	rlp.p2 = 0;

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