home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / rsynth / src / nsynth.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  32.2 KB  |  1,278 lines
  1.  
  2. /* Copyright            1982                    by Dennis H. Klatt
  3.  *
  4.  *      Klatt synthesizer
  5.  *         Modified version of synthesizer described in
  6.  *         J. Acoust. Soc. Am., Mar. 1980. -- new voicing
  7.  *         source.
  8.  *
  9.  * Edit history
  10.  * 000001 10-Mar-83 DK  Initial creation.
  11.  * 000002  5-May-83 DK  Fix bug in operation of parallel F1
  12.  * 000003  7-Jul-83 DK  Allow parallel B1 to vary, and if ALL_PARALLEL,
  13.  *                      also allow B2 and B3 to vary
  14.  * 000004 26-Jul-83 DK  Get rid of mulsh, use short for VAX
  15.  * 000005 24-Oct-83 DK  Split off parwavtab.c, change short to int
  16.  * 000006 16-Nov-83 DK  Make samrate a variable, use exp(), cos() rand()
  17.  * 000007 17-Nov-83 DK  Convert to float, remove  cpsw, add set outsl
  18.  * 000008 28-Nov-83 DK  Add simple impulsive glottal source option
  19.  * 000009  7-Dec-83 DK  Use spkrdef[7] to select impulse or natural voicing
  20.  *                       and update cascade F1,..,F6 at update times
  21.  * 000010 19-Dec-83 DK  Add subroutine no_rad_char() to get rid of rad char
  22.  * 000011 28-Jan-84 DK  Allow up to 8 formants in cascade branch F7 fixed
  23.  *                       at 6.5 kHz, F8 fixed at 7.5 kHz
  24.  * 000012 14-Feb-84 DK  Fix bug in 'os' options so os>12 works
  25.  * 000013 17-May-84 DK  Add G0 code
  26.  * 000014 12-Mar-85 DHW modify for Haskins environment
  27.  * 000015 11-Jul-87 LG  modificiations for PC
  28.  * 000016 20-Apr-91 ATS Modified for SPARCSTATION
  29.  */
  30.  
  31. #define SUSPECT
  32. #define KLATT_USE_OUTS
  33.  
  34. #include <stdio.h>
  35. #include <stdlib.h>
  36. #include <math.h>
  37. #include "proto.h"
  38. #include "parwave.h"
  39.  
  40. extern int rand PROTO((void));
  41. /* BENOIST
  42. extern int srand PROTO((unsigned int));
  43. */
  44.  
  45. typedef struct
  46.  {
  47.   char *name;
  48.   float a;
  49.   float b;
  50.   float c;
  51.   float p1;
  52.   float p2;
  53.  }
  54. resonator_t, *resonator_ptr;
  55.  
  56. /* VARIABLES TO HOLD SPEAKER DEFINITION FROM HOST:                    */
  57.  
  58. /* N.I-S 23rd June 1993
  59.    Moved definitions of some globals here, and hence to parawav.o as
  60.    they are needed by all users.
  61.    We can set sensible defaults ...
  62. */
  63.  
  64. int quiet_flag = TRUE;
  65. int f0_flutter = 0;
  66. long initsw = 0;
  67. long warnsw = 0;
  68. long dispt = 0;
  69. long disptcum = 0;
  70. long outsl = 0;  /* Output waveform selector - normal  */
  71. long glsource = NATURAL;
  72.                                 /* 1->impulsive, 2->natural voicing source      */
  73. int time_count = 0;
  74.                                 /* JPI added for f0 flutter */
  75. klatt_t pars;      /* Incoming parameters */
  76. long sigmx = 0;  /* Maximum level reached */
  77.  
  78. /* COUNTERS */
  79.  
  80. static long nper;
  81.                                 /* Current loc in voicing period   40000 samp/s */
  82.  
  83. /* COUNTER LIMITS */
  84.  
  85. static long T0;  /* Fundamental period in output samples times 4 */
  86. static long nopen;
  87.                                 /* Number of samples in open phase of period    */
  88. static long nmod;
  89.                                 /* Position in period to begin noise amp. modul */
  90.  
  91. /* VARIABLES THAT HAVE TO STICK AROUND FOR AWHILE, AND THUS LOCALS */
  92. /* ARE NOT APPROPRIATE */
  93.  
  94. /* Incoming parameter Variables which need to be updated synchronously  */
  95.  
  96. static long F0hz10;
  97.                                 /* Voicing fund freq in Hz                       */
  98. static long AVdb;
  99.                                 /* Amp of voicing in dB,    0 to   70           */
  100. static long Aturb;
  101.                                 /* Breathiness in voicing,      0 to   80       */
  102. static long TLTdb;
  103.                                 /* Voicing spectral tilt in dB,  0 to   24      */
  104. static long Kopen;
  105.                                 /* # of samples in open period, 10 to   65      */
  106. static long Kskew;
  107.                                 /* Skewness of alternate periods,0 to   40      */
  108.  
  109. /* Various amplitude variables used in main loop */
  110.  
  111. static float amp_voice;         /* AVdb converted to linear gain            */
  112. static float amp_bypas;         /* AB converted to linear gain              */
  113. static float par_amp_voice;     /* AVpdb converted to linear gain           */
  114. static float amp_aspir;         /* AP converted to linear gain              */
  115. static float amp_frica;         /* AF converted to linear gain              */
  116. static float amp_breth;         /* ATURB converted to linear gain           */
  117.  
  118. /* State variables of sound sources */
  119.  
  120. static long nrand;
  121.                                 /* Varible used by random number generator      */
  122. static long skew;
  123.                                 /* Alternating jitter, in half-period units     */
  124.  
  125. static float a;  /* Makes waveshape of glottal pulse when open   */
  126. static float b;  /* Makes waveshape of glottal pulse when open   */
  127. static float vwave;
  128.                                 /* Ditto, but before multiplication by AVdb     */
  129. static float vlast;
  130.                                 /* Previous output of voice                     */
  131. static float nlast;
  132.                                 /* Previous output of random number generator   */
  133. static float glotlast;
  134.                                 /* Previous value of glotout                    */
  135. static float decay;
  136.                                 /* TLTdb converted to exponential time const    */
  137. static float onemd;
  138.                                 /* in voicing one-pole low-pass filter          */
  139. static float minus_pi_t;
  140.                                 /* func. of sample rate */
  141. static float two_pi_t;
  142.                                 /* func. of sample rate */
  143.  
  144.  
  145. /* INTERNAL MEMORY FOR DIGITAL RESONATORS AND ANTIRESONATOR  */
  146.  
  147. static resonator_t rnpp =
  148. {"parallel nasal pole"};
  149. static resonator_t r1p =
  150. {"parallel 1st formant"};
  151. static resonator_t r2p =
  152. {"parallel 2nd formant"};
  153. static resonator_t r3p =
  154. {"parallel 3rd formant"};
  155. static resonator_t r4p =
  156. {"parallel 4th formant"};
  157. static resonator_t r5p =
  158. {"parallel 5th formant"};
  159. static resonator_t r6p =
  160. {"parallel 6th formant"};
  161. static resonator_t r1c =
  162. {"cascade 1st formant"};
  163. static resonator_t r2c =
  164. {"cascade 2nd formant"};
  165. static resonator_t r3c =
  166. {"cascade 3rd formant"};
  167. static resonator_t r4c =
  168. {"cascade 4th formant"};
  169. static resonator_t r5c =
  170. {"cascade 5th formant"};
  171. static resonator_t r6c =
  172. {"cascade 6th formant"};
  173. static resonator_t r7c =
  174. {"cascade 7th formant"};
  175. static resonator_t r8c =
  176. {"cascade 8th formant"};
  177. static resonator_t rnpc =
  178. {"cascade nasal pole"};
  179. static resonator_t rnz =
  180. {"cascade nasal zero"};
  181. static resonator_t rgl =
  182. {"crit-damped glot low-pass filter"};
  183. static resonator_t rlp =
  184. {"downsamp low-pass filter"};
  185. static resonator_t rout =
  186. {"output low-pass"};
  187.  
  188. /*
  189.  * Constant B0 controls shape of glottal pulse as a function
  190.  * of desired duration of open phase N0
  191.  * (Note that N0 is specified in terms of 40,000 samples/sec of speech)
  192.  *
  193.  *    Assume voicing waveform V(t) has form: k1 t**2 - k2 t**3
  194.  *
  195.  *    If the radiation characterivative, a temporal derivative
  196.  *      is folded in, and we go from continuous time to discrete
  197.  *      integers n:  dV/dt = vwave[n]
  198.  *                         = sum over i=1,2,...,n of { a - (i * b) }
  199.  *                         = a n  -  b/2 n**2
  200.  *
  201.  *      where the  constants a and b control the detailed shape
  202.  *      and amplitude of the voicing waveform over the open
  203.  *      potion of the voicing cycle "nopen".
  204.  *
  205.  *    Let integral of dV/dt have no net dc flow --> a = (b * nopen) / 3
  206.  *
  207.  *    Let maximum of dUg(n)/dn be constant --> b = gain / (nopen * nopen)
  208.  *      meaning as nopen gets bigger, V has bigger peak proportional to n
  209.  *
  210.  *    Thus, to generate the table below for 40 <= nopen <= 263:
  211.  *
  212.  *      B0[nopen - 40] = 1920000 / (nopen * nopen)
  213.  */
  214. static const short B0[224] =
  215. {
  216.  1200, 1142, 1088, 1038, 991,
  217.  948, 907, 869, 833, 799,
  218.  768, 738, 710, 683, 658,
  219.  634, 612, 590, 570, 551,
  220.  533, 515, 499, 483, 468,
  221.  454, 440, 427, 415, 403,
  222.  391, 380, 370, 360, 350,
  223.  341, 332, 323, 315, 307,
  224.  300, 292, 285, 278, 272,
  225.  265, 259, 253, 247, 242,
  226.  237, 231, 226, 221, 217,
  227.  212, 208, 204, 199, 195,
  228.  192, 188, 184, 180, 177,
  229.  174, 170, 167, 164, 161,
  230.  158, 155, 153, 150, 147,
  231.  145, 142, 140, 137, 135,
  232.  133, 131, 128, 126, 124,
  233.  122, 120, 119, 117, 115,
  234.  113, 111, 110, 108, 106,
  235.  105, 103, 102, 100, 99,
  236.  97, 96, 95, 93, 92,
  237.  91, 90, 88, 87, 86,
  238.  85, 84, 83, 82, 80,
  239.  79, 78, 77, 76, 75,
  240.  75, 74, 73, 72, 71,
  241.  70, 69, 68, 68, 67,
  242.  66, 65, 64, 64, 63,
  243.  62, 61, 61, 60, 59,
  244.  59, 58, 57, 57, 56,
  245.  56, 55, 55, 54, 54,
  246.  53, 53, 52, 52, 51,
  247.  51, 50, 50, 49, 49,
  248.  48, 48, 47, 47, 46,
  249.  46, 45, 45, 44, 44,
  250.  43, 43, 42, 42, 41,
  251.  41, 41, 41, 40, 40,
  252.  39, 39, 38, 38, 38,
  253.  38, 37, 37, 36, 36,
  254.  36, 36, 35, 35, 35,
  255.  35, 34, 34, 33, 33,
  256.  33, 33, 32, 32, 32,
  257.  32, 31, 31, 31, 31,
  258.  30, 30, 30, 30, 29,
  259.  29, 29, 29, 28, 28,
  260.  28, 28, 27, 27
  261. };
  262.  
  263. /*
  264.  * Convertion table, db to linear, 87 dB --> 32767
  265.  *                                 86 dB --> 29491 (1 dB down = 0.5**1/6)
  266.  *                                 ...
  267.  *                                 81 dB --> 16384 (6 dB down = 0.5)
  268.  *                                 ...
  269.  *                                  0 dB -->     0
  270.  *
  271.  * The just noticeable difference for a change in intensity of a vowel
  272.  *   is approximately 1 dB.  Thus all amplitudes are quantized to 1 dB
  273.  *   steps.
  274.  */
  275.  
  276. static const float amptable[88] =
  277. {
  278.  0., 0., 0., 0., 0.,
  279.  0., 0., 0., 0., 0.,
  280.  0., 0., 0., 6., 7.,
  281.  8., 9., 10., 11., 13.,
  282.  14., 16., 18., 20., 22.,
  283.  25., 28., 32., 35., 40.,
  284.  45., 51., 57., 64., 71.,
  285.  80., 90., 101., 114., 128.,
  286.  142., 159., 179., 202., 227.,
  287.  256., 284., 318., 359., 405.,
  288.  455., 512., 568., 638., 719.,
  289.  811., 911., 1024., 1137., 1276.,
  290.  1438., 1622., 1823., 2048., 2273.,
  291.  2552., 2875., 3244., 3645., 4096.,
  292.  4547., 5104., 5751., 6488., 7291.,
  293.  8192., 9093., 10207., 11502., 12976.,
  294.  14582., 16384., 18350., 20644., 23429.,
  295.  26214., 29491., 32767.
  296. };
  297.  
  298. const char *par_name[] =
  299. {
  300.  "F0hz10",
  301.  "AVdb",
  302.  "F1hz",
  303.  "B1hz",
  304.  "F2hz",
  305.  "B2hz",
  306.  "F3hz",
  307.  "B3hz",
  308.  "F4hz",
  309.  "B4hz",
  310.  "F5hz",
  311.  "B5hz",
  312.  "F6hz",
  313.  "B6hz",
  314.  "FNZhz",
  315.  "BNZhz",
  316.  "FNPhz",
  317.  "BNPhz",
  318.  "AP",
  319.  "Kopen",
  320.  "Aturb",
  321.  "TLTdb",
  322.  "AF",
  323.  "Kskew",
  324.  "A1",
  325.  "B1phz",
  326.  "A2",
  327.  "B2phz",
  328.  "A3",
  329.  "B3phz",
  330.  "A4",
  331.  "B4phz",
  332.  "A5",
  333.  "B5phz",
  334.  "A6",
  335.  "B6phz",
  336.  "ANP",
  337.  "AB",
  338.  "AVpdb",
  339.  "Gain0"
  340. };
  341.  
  342. /*
  343. function FLUTTER
  344.  
  345. This function adds F0 flutter, as specified in:
  346.  
  347. "Analysis, synthesis and perception of voice quality variations among
  348. female and male talkers" D.H. Klatt and L.C. Klatt JASA 87(2) February 1990.
  349. Flutter is added by applying a quasi-random element constructed from three
  350. slowly varying sine waves.
  351.  
  352. */
  353.  
  354. static void flutter PROTO((klatt_ptr pars));
  355. static void
  356. flutter(pars)
  357. klatt_ptr pars;
  358. {
  359.  long original_f0;
  360.  double delta_f0;
  361.  double fla, flb, flc, fld, fle;
  362.  
  363.  original_f0 = pars->F0hz10 / 10;
  364.  
  365.  fla = (double) f0_flutter / 50;
  366.  flb = (double) original_f0 / 100;
  367.  flc = sin(2 * PI * 12.7 * time_count);
  368.  fld = sin(2 * PI * 7.1 * time_count);
  369.  fle = sin(2 * PI * 4.7 * time_count);
  370.  delta_f0 = fla * flb * (flc + fld + fle) * 10;
  371.  F0hz10 = F0hz10 + (long) delta_f0;
  372. }
  373.  
  374. static float resonator PROTO((resonator_ptr r, float input));
  375.  
  376. /* Generic resonator function */
  377. static float
  378. resonator(r, input)
  379. resonator_ptr r;
  380. float input;
  381. {
  382.  register float x = r->a * input + r->b * r->p1 + r->c * r->p2;
  383.  r->p2 = r->p1;
  384.  r->p1 = x;
  385.  return x;
  386. }
  387.  
  388. /* Generic anti-resonator function
  389.    Same as resonator except that a,b,c need to be set with setzeroabc()
  390.    and we save inputs in p1/p2 rather than outputs.
  391.    There is currently only one of these - "rnz"
  392. */
  393. /*  Output = (rnz.a * input) + (rnz.b * oldin1) + (rnz.c * oldin2) */
  394. static float antiresonator PROTO((resonator_ptr r, float input));
  395.  
  396. static float
  397. antiresonator(r, input)
  398. resonator_ptr r;
  399. float input;
  400. {
  401.  register float x = r->a * input + r->b * r->p1 + r->c * r->p2;
  402.  r->p2 = r->p1;
  403.  r->p1 = input;
  404.  return x;
  405. }
  406.  
  407. static float impulsive_source PROTO((long nper));
  408.  
  409. static float
  410. impulsive_source(nper)
  411. long nper;
  412. {
  413.  static float doublet[] =
  414.  {0., 13000000., -13000000.};
  415.  if (nper < 3)
  416.   {
  417.    vwave = doublet[nper];
  418.   }
  419.  else
  420.   {
  421.    vwave = 0.0;
  422.   }
  423.  /* Low-pass filter the differenciated impulse with a critically-damped
  424.     second-order filter, time constant proportional to Kopen */
  425.  return resonator(&rgl, vwave);
  426. }
  427.  
  428.  
  429. /* Vwave is the differentiated glottal flow waveform, there is a weak
  430.     spectral zero around 800 Hz, magic constants a,b reset pitch-synch */
  431.  
  432. static float natural_source PROTO((long nper));
  433.  
  434. static float
  435. natural_source(nper)
  436. long nper;
  437. {
  438.  float lgtemp;
  439.  /* See if glottis open */
  440.  if (nper < nopen)
  441.   {
  442.    a -= b;
  443.    vwave += a;
  444.    lgtemp = vwave * 0.028;
  445.    return (lgtemp);
  446.   }
  447.  else
  448.   {
  449.    /* Glottis closed */
  450.    vwave = 0.;
  451.    return (0.0);
  452.   }
  453. }
  454.  
  455. /* Convert formant freqencies and bandwidth into
  456.    resonator difference equation coefficents
  457. */
  458.  
  459. static void setabc PROTO((long int f, long int bw, resonator_ptr rp));
  460.  
  461. static void
  462. setabc(f, bw, rp)
  463. long int f;          /* Frequency of resonator in Hz         */
  464. long int bw;        /* Bandwidth of resonator in Hz         */
  465. resonator_ptr rp;
  466.                                 /* Are output coefficients              */
  467. {
  468.  float r;
  469.  double arg;
  470.  
  471. /* Let r  =  exp(-pi bw t) */
  472.  
  473.  arg = minus_pi_t * bw;
  474.  r = exp(arg);
  475.  
  476. /* Let c  =  -r**2 */
  477.  
  478.  rp->c = -(r * r);
  479.  
  480. /* Let b = r * 2*cos(2 pi f t) */
  481.  
  482.  arg = two_pi_t * f;
  483.  rp->b = r * cos(arg) * 2.;
  484.  
  485. /* Let a = 1.0 - b - c */
  486.  
  487.  rp->a = 1.0 - rp->b - rp->c;
  488. #if 0
  489.  printf("%20s %6.4g %6.4g %6.4g\n", rp->name, rp->a, rp->b, rp->c);
  490. #endif
  491. }
  492.  
  493. /* Convert formant freqencies and bandwidth into
  494.  *      anti-resonator difference equation constants */
  495.  
  496. static void setzeroabc PROTO((long int f, long int bw, resonator_ptr rp));
  497.  
  498. static void
  499. setzeroabc(f, bw, rp)
  500. long int f;          /* Frequency of resonator in Hz         */
  501. long int bw;        /* Bandwidth of resonator in Hz         */
  502. resonator_ptr rp;
  503.                                 /* Are output coefficients              */
  504. {
  505.  float r;
  506.  double arg;
  507.  
  508. /* First compute ordinary resonator coefficients */
  509. /* Let r  =  exp(-pi bw t) */
  510.  
  511.  arg = minus_pi_t * bw;
  512.  r = exp(arg);
  513.  
  514. /* Let c  =  -r**2 */
  515.  
  516.  rp->c = -(r * r);
  517.  
  518. /* Let b = r * 2*cos(2 pi f t) */
  519.  
  520.  arg = two_pi_t * f;
  521.  rp->b = r * cos(arg) * 2.;
  522.  
  523. /* Let a = 1.0 - b - c */
  524.  
  525.  rp->a = 1.0 - rp->b - rp->c;
  526.  
  527. /* Now convert to antiresonator coefficients (a'=1/a, b'=b/a, c'=c/a) */
  528.  
  529.  rp->a = 1.0 / rp->a;
  530.  rp->c *= -rp->a;
  531.  rp->b *= -rp->a;
  532. }
  533.  
  534.  
  535. /* Random number generator (return a number between -8191 and +8191) */
  536. static float gen_noise PROTO((void));
  537.  
  538. static float
  539. gen_noise()
  540. {
  541.  float noise;
  542.  long temp = rand();
  543.  nrand = (temp >> 17) - 8191;
  544.  
  545. /* Tilt down noise spectrum by soft low-pass filter having
  546.  *    a pole near the origin in the z-plane, i.e.
  547.  *    output = input + (0.75 * lastoutput) */
  548.  
  549.  noise = nrand + (0.75 * nlast);
  550.  nlast = noise;
  551.  return noise;
  552. }
  553.  
  554. /* Convert from decibels to a linear scale factor */
  555.  
  556. static float DBtoLIN PROTO((long int dB));
  557.  
  558. static float
  559. DBtoLIN(dB)
  560. long int dB;
  561. {
  562.  /* Check limits or argument (can be removed in final product) */
  563.  if (dB < 0)
  564.   dB = 0;
  565.  else if (dB >= 88)
  566.   {
  567.    if (!quiet_flag)
  568.     printf("Try to compute amptable[%d]\n", dB);
  569.    dB = 87;
  570.   }
  571.  return amptable[dB] * 0.001;
  572. }
  573.  
  574. #define ACOEF           0.005
  575. #define BCOEF           (1.0 - ACOEF) /* Slight decay to remove dc */
  576.  
  577. #if 0
  578. static float no_rad_char PROTO((float in));
  579.  
  580. static float
  581. no_rad_char(in)
  582. float in;
  583. {
  584.  static float lastin;
  585.  float out = (ACOEF * in) + (BCOEF * lastin);
  586.  lastin = out;
  587.  out = -100. * out;
  588.  /* Scale up to make visible */
  589.  return out;
  590. }
  591. #endif
  592.  
  593. static float dBconvert PROTO((long int arg));
  594.  
  595. static float
  596. dBconvert(arg)
  597. long int arg;
  598. {
  599.  return 20.0 * log10((double) arg / 32767.0);
  600. }
  601.  
  602. static void overload_warning PROTO((long int arg));
  603.  
  604. static void
  605. overload_warning(arg)
  606. long int arg;
  607. {
  608.  if (warnsw == 0)
  609.   {
  610.    warnsw++;
  611.    if (!quiet_flag)
  612.     {
  613.      printf("\n* * * WARNING: ");
  614.      printf(" Signal at output of synthesizer (+%3.1f dB) exceeds 0 dB\n",
  615.             dBconvert(arg));
  616.      printf("    at synthesis parameter frame = %d\n", disptcum / nspfr);
  617.     }
  618.   }
  619. }
  620.  
  621. static short clip PROTO((float input));
  622.  
  623. static short
  624. clip(input)
  625. float input;
  626. {
  627.  long temp = input;
  628.  long x = (temp < 0) ? -temp : temp;
  629.  if (x > sigmx)
  630.   sigmx = x;
  631.  
  632.  /* clip on boundaries of 16-bit word */
  633.  if (temp < -32767)
  634.   {
  635.    overload_warning(x);
  636.    temp = -32767;
  637.   }
  638.  else if (temp > 32767)
  639.   {
  640.    overload_warning(temp);
  641.    temp = 32767;
  642.   }
  643.  return (temp);
  644. }
  645.  
  646. static void resonclr PROTO((void));
  647.  
  648. static void
  649. resonclr()
  650. {
  651.  rnpp.p1 = 0;       /* parallel nasal pole  */
  652.  rnpp.p2 = 0;
  653.  
  654.  r1p.p1 = 0;         /* parallel 1st formant */
  655.  r1p.p2 = 0;
  656.  
  657.  r2p.p1 = 0;         /* parallel 2nd formant */
  658.  r2p.p2 = 0;
  659.  
  660.  r3p.p1 = 0;         /* parallel 3rd formant */
  661.  r3p.p2 = 0;
  662.  
  663.  r4p.p1 = 0;         /* parallel 4th formant */
  664.  r4p.p2 = 0;
  665.  
  666.  r5p.p1 = 0;         /* parallel 5th formant */
  667.  r5p.p2 = 0;
  668.  
  669.  r6p.p1 = 0;         /* parallel 6th formant */
  670.  r6p.p2 = 0;
  671.  
  672.  r1c.p1 = 0;         /* cascade 1st formant  */
  673.  r1c.p2 = 0;
  674.  
  675.  r2c.p1 = 0;         /* cascade 2nd formant  */
  676.  r2c.p2 = 0;
  677.  
  678.  r3c.p1 = 0;         /* cascade 3rd formant  */
  679.  r3c.p2 = 0;
  680.  
  681.  r4c.p1 = 0;         /* cascade 4th formant  */
  682.  r4c.p2 = 0;
  683.  
  684.  r5c.p1 = 0;         /* cascade 5th formant  */
  685.  r5c.p2 = 0;
  686.  
  687.  r6c.p1 = 0;         /* cascade 6th formant  */
  688.  r6c.p2 = 0;
  689.  
  690.  r7c.p1 = 0;
  691.  r7c.p2 = 0;
  692.  
  693.  r8c.p1 = 0;
  694.  r8c.p2 = 0;
  695.  
  696.  rnpc.p1 = 0;       /* cascade nasal pole   */
  697.  rnpc.p2 = 0;
  698.  
  699.  rnz.p1 = 0;         /* cascade nasal zero   */
  700.  rnz.p2 = 0;
  701.  
  702.  rgl.p1 = 0;         /* crit-damped glot low-pass filter */
  703.  rgl.p2 = 0;
  704.  
  705.  rlp.p1 = 0;         /* downsamp low-pass filter    */
  706.  rlp.p2 = 0;
  707.  
  708.  vlast = 0;           /* Previous output of voice                     */
  709.  nlast = 0;           /* Previous output of random number generator   */
  710.  glotlast = 0;     /* Previous value of glotout                   */
  711. }
  712.  
  713.  
  714. /* Reset selected parameters pitch-synchronously */
  715.  
  716. static void pitch_synch_par_reset PROTO((long ns));
  717.  
  718. static void
  719. pitch_synch_par_reset(ns)
  720. long ns;
  721. {
  722.  long temp;
  723.  float temp1;
  724.  if (F0hz10 > 0)
  725.   {
  726.    T0 = (40 * samp_rate) / F0hz10;
  727.    /* Period in samp*4 */
  728.    amp_voice = DBtoLIN(AVdb);
  729.  
  730.    /* Duration of period before amplitude modulation */
  731.    nmod = T0;
  732.    if (AVdb > 0)
  733.     {
  734.      nmod >>= 1;
  735.     }
  736.  
  737.    /* Breathiness of voicing waveform */
  738.  
  739.    amp_breth = DBtoLIN(Aturb) * 0.1;
  740.  
  741.    /* Set open phase of glottal period */
  742.    /* where  40 <= open phase <= 263 */
  743.  
  744.    nopen = 4 * Kopen;
  745.    if ((glsource == IMPULSIVE) && (nopen > 263))
  746.     {
  747.      nopen = 263;
  748.     }
  749.  
  750.    if (nopen >= (T0 - 1))
  751.     {
  752.      nopen = T0 - 2;
  753.      if (!quiet_flag)
  754.       {
  755.        printf("Warning: glottal open period cannot exceed T0, truncated\n");
  756.       }
  757.     }
  758.  
  759.    if (nopen < 40)
  760.     {
  761.      nopen = 40;     /* F0 max = 1000 Hz */
  762.      if (!quiet_flag)
  763.       {
  764.        printf("Warning: minimum glottal open period is 10 samples.\n");
  765.        printf("truncated, nopen = %d\n", nopen);
  766.       }
  767.     }
  768.  
  769.    /* Reset a & b, which determine shape of "natural" glottal waveform */
  770.  
  771.    b = B0[nopen - 40];
  772.    a = (b * nopen) * .333;
  773.  
  774.    /* Reset width of "impulsive" glottal pulse */
  775.  
  776.    temp = samp_rate / nopen;
  777.    setabc(0L, temp, &rgl);
  778.  
  779.    /* Make gain at F1 about constant */
  780.  
  781.    temp1 = nopen * .00833;
  782.    rgl.a *= (temp1 * temp1);
  783.  
  784.    /* Truncate skewness so as not to exceed duration of closed phase
  785.       of glottal period */
  786.  
  787.    temp = T0 - nopen;
  788.    if (Kskew > temp)
  789.     {
  790.      if (!quiet_flag)
  791.       {
  792.        printf("Kskew duration=%d > glottal closed period=%d, truncate\n",
  793.               Kskew, T0 - nopen);
  794.       }
  795.      Kskew = temp;
  796.     }
  797.    if (skew >= 0)
  798.     {
  799.      skew = Kskew; /* Reset skew to requested Kskew */
  800.     }
  801.    else
  802.     {
  803.      skew = -Kskew;
  804.     }
  805.  
  806.    /* Add skewness to closed portion of voicing period */
  807.  
  808.    T0 = T0 + skew;
  809.    skew = -skew;
  810.   }
  811.  else
  812.   {
  813.    T0 = 4;               /* Default for f0 undefined */
  814.    amp_voice = 0.;
  815.    nmod = T0;
  816.    amp_breth = 0.;
  817.    a = 0.0;
  818.    b = 0.0;
  819.   }
  820.  
  821. /* Reset these pars pitch synchronously or at update rate if f0=0 */
  822.  
  823.  if ((T0 != 4) || (ns == 0))
  824.   {
  825.    /* Set one-pole low-pass filter that tilts glottal source */
  826.    decay = (0.033 * TLTdb);
  827.    if (decay > 0.0)
  828.     {
  829.      onemd = 1.0 - decay;
  830.     }
  831.    else
  832.     {
  833.      onemd = 1.0;
  834.     }
  835.   }
  836. }
  837.  
  838. void show_parms PROTO((int *pars));
  839.  
  840. void
  841. show_parms(pars)
  842. int *pars;
  843. {
  844.  int i;
  845.  if (!initsw)
  846.   {
  847.    for (i = 0; i < NPAR; i++)
  848.     printf("%s ", par_name[i]);
  849.    printf("\n");
  850.   }
  851.  for (i = 0; i < NPAR; i++)
  852.   {
  853.    printf("%*d ", (int) strlen(par_name[i]), pars[i]);
  854.   }
  855.  printf("\n");
  856. }
  857.  
  858. /* Get variable parameters from host computer,
  859.    initially also get definition of fixed pars
  860. */
  861.  
  862. static void gethost PROTO((klatt_ptr pars));
  863. static void
  864. gethost(pars)
  865. klatt_ptr pars;
  866. {
  867.  /* VARIABLES TO HOLD INPUT PARAMETERS FROM HOST:                      */
  868.  long F1hz;           /* First formant freq in Hz,  200 to 1300      */
  869.  long F2hz;           /* Second formant freq in Hz,  550 to 3000     */
  870.  long F3hz;           /* Third formant freq in Hz, 1200 to 4999      */
  871.  long F4hz;           /* Fourth formant freq in Hz, 1200 to 4999     */
  872.  long F5hz;           /* Fifth formant freq in Hz, 1200 to 4999      */
  873.  long F6hz;           /* Sixth formant freq in Hz, 1200 to 4999      */
  874.  long FNZhz;         /* Nasal zero freq in Hz,  248 to  528      */
  875.  long FNPhz;         /* Nasal pole freq in Hz,  248 to  528      */
  876.  
  877.  long B1hz;           /* First formant bw in Hz,   40 to 1000      */
  878.  long B2hz;           /* Second formant bw in Hz,   40 to 1000      */
  879.  long B3hz;           /* Third formant bw in Hz,   40 to 1000      */
  880.  long B4hz;           /* Fourth formant bw in Hz,  40 to 1000      */
  881.  long B5hz;           /* Fifth formant bw in Hz,   40 to 1000      */
  882.  long B6hz;           /* Sixth formant bw in Hz,   40 to 2000      */
  883.  long B1phz;         /* Par. 1st formant bw in Hz,   40 to 1000     */
  884.  long B2phz;         /* Par. 2nd formant bw in Hz,   40 to 1000     */
  885.  long B3phz;         /* Par. 3rd formant bw in Hz,   40 to 1000     */
  886.  long B4phz;         /* Par. 4th formant bw in Hz,  40 to 1000      */
  887.  long B5phz;         /* Par. 5th formant bw in Hz,   40 to 1000     */
  888.  long B6phz;         /* Par. 6th formant bw in Hz,   40 to 2000     */
  889.  long BNZhz;         /* Nasal zero bw in Hz,   40 to 1000      */
  890.  long BNPhz;         /* Nasal pole bw in Hz,   40 to 1000      */
  891.  
  892.  long AVpdb;         /* Amp of voicing, par in dB,    0 to   70      */
  893.  long AF;               /* Amp of frication in dB,    0 to   80         */
  894.  long A1;               /* Amp of par 1st formant in dB, 0 to   80      */
  895.  long ANP;             /* Amp of par nasal pole in dB,  0 to   80      */
  896.  long A2;               /* Amp of F2 frication in dB,    0 to   80      */
  897.  long A3;               /* Amp of F3 frication in dB,    0 to   80      */
  898.  long A4;               /* Amp of F4 frication in dB,    0 to   80      */
  899.  long A5;               /* Amp of F5 frication in dB,    0 to   80      */
  900.  long A6;               /* Amp of F6 (same as r6p.a),     0 to   80     */
  901.  long AB;               /* Amp of bypass fric. in dB,    0 to   80      */
  902.  /* in sample#/2                                 */
  903.  long Gain0;         /* Overall gain, 60 dB is unity  0 to   60      */
  904.  
  905.  float amp_parF1;
  906.                                 /* A1 converted to linear gain              */
  907.  float amp_parFNP;
  908.                                 /* AP converted to linear gain              */
  909.  float amp_parF2;
  910.                                 /* A2 converted to linear gain              */
  911.  float amp_parF3;
  912.                                 /* A3 converted to linear gain              */
  913.  float amp_parF4;
  914.                                 /* A4 converted to linear gain              */
  915.  float amp_parF5;
  916.                                 /* A5 converted to linear gain              */
  917.  float amp_parF6;
  918.                                 /* A6 converted to linear gain              */
  919.  
  920. #if 0
  921.  show_parms((int *) pars);
  922. #endif
  923.  
  924.  if (initsw == 0)
  925.   {
  926.    unsigned useed = 1;          /* for Lattice random number gen LG */
  927.    /* Initialize speaker definition */
  928.    long FLPhz = (950 * samp_rate) / 10000;
  929.    long BLPhz = (630 * samp_rate) / 10000;
  930.  
  931.    minus_pi_t = -PI / samp_rate;
  932.    two_pi_t = -2.0 * minus_pi_t;
  933.  
  934.    setabc(FLPhz, BLPhz, &rlp);
  935.    srand(useed);      /* Init ran # generator */
  936.    resonclr();        /* LG adds 7/11/87 initializes resonantors for new utterance  */
  937.    nper = 0;          /* LG */
  938.    T0 = 0;            /* LG */
  939.    initsw++;
  940.   }
  941.  
  942.  /*
  943.    Read  speech frame definition into temp store
  944.    and move some parameters into active use immediately
  945.    (voice-excited ones are updated pitch synchronously
  946.    to avoid waveform glitches).
  947.  */
  948.  
  949.  
  950.  F0hz10 = pars->F0hz10;
  951.  AVdb = pars->AVdb - 7;
  952.  if (AVdb < 0)
  953.   AVdb = 0;
  954.  
  955.  F1hz = pars->F1hz;
  956.  B1hz = pars->B1hz;
  957.  
  958.  F2hz = pars->F2hz;
  959.  B2hz = pars->B2hz;
  960.  
  961.  F3hz = pars->F3hz;
  962.  B3hz = pars->B3hz;
  963.  
  964.  F4hz = pars->F4hz;
  965.  B4hz = pars->B4hz;
  966.  
  967.  F5hz = pars->F5hz;
  968.  B5hz = pars->B5hz;
  969.  
  970.  F6hz = pars->F6hz;
  971.                                 /* f  of parallel 6th formant */
  972.  B6hz = pars->B6hz;
  973.                                 /* bw of cascade 6th formant (doesn't exist) */
  974.  
  975.  FNZhz = pars->FNZhz;
  976.  BNZhz = pars->BNZhz;
  977.  
  978.  FNPhz = pars->FNPhz;
  979.  BNPhz = pars->BNPhz;
  980.  
  981.  amp_aspir = DBtoLIN(pars->ASP) * .05;
  982.  Kopen = pars->Kopen;
  983.  
  984.  Aturb = pars->Aturb;
  985.  TLTdb = pars->TLTdb;
  986.  
  987.  AF = pars->AF;
  988.  amp_frica = DBtoLIN(AF) * 0.25;
  989.  Kskew = pars->Kskew;
  990.  AVpdb = pars->AVpdb;
  991.  par_amp_voice = DBtoLIN(AVpdb);
  992.  
  993.  A1 = pars->A1;
  994.  amp_parF1 = DBtoLIN(A1) * 0.4;
  995.  B1phz = pars->B1phz;
  996.  
  997.  A2 = pars->A2;
  998.  amp_parF2 = DBtoLIN(A2) * 0.15;
  999.  B2phz = pars->B2phz;
  1000.  
  1001.  A3 = pars->A3;
  1002.  amp_parF3 = DBtoLIN(A3) * 0.06;
  1003.  B3phz = pars->B3phz;
  1004.  
  1005.  A4 = pars->A4;
  1006.  amp_parF4 = DBtoLIN(A4) * 0.04;
  1007.  B4phz = pars->B4phz;
  1008.  
  1009.  A5 = pars->A5;
  1010.  amp_parF5 = DBtoLIN(A5) * 0.022;
  1011.  B5phz = pars->B5phz;
  1012.  
  1013.  A6 = pars->A6;
  1014.  amp_parF6 = DBtoLIN(A6) * 0.03;
  1015.  B6phz = pars->B6phz;
  1016.  
  1017.  ANP = pars->ANP;
  1018.  amp_parFNP = DBtoLIN(ANP) * 0.6;
  1019.  
  1020.  AB = pars->AB;
  1021.  amp_bypas = DBtoLIN(AB) * 0.05;
  1022.  
  1023.  if (nfcascade >= 8)
  1024.   {
  1025.    /* Inside Nyquist rate ? */
  1026.    if (samp_rate >= 16000)
  1027.     setabc(7500, 600, &r8c);
  1028.    else
  1029.     nfcascade = 6;
  1030.   }
  1031.  
  1032.  if (nfcascade >= 7)
  1033.   {
  1034.    /* Inside Nyquist rate ? */
  1035.    if (samp_rate >= 16000)
  1036.     setabc(6500, 500, &r7c);
  1037.    else
  1038.     nfcascade = 6;
  1039.   }
  1040.  
  1041.  /* Set coefficients of variable cascade resonators */
  1042.  
  1043.  if (nfcascade >= 6)
  1044.   setabc(F6hz, B6hz, &r6c);
  1045.  if (nfcascade >= 5)
  1046.   setabc(F5hz, B5hz, &r5c);
  1047.  
  1048.  setabc(F4hz, B4hz, &r4c);
  1049.  setabc(F3hz, B3hz, &r3c);
  1050.  setabc(F2hz, B2hz, &r2c);
  1051.  setabc(F1hz, B1hz, &r1c);
  1052.  
  1053.  /* Set coeficients of nasal resonator and zero antiresonator */
  1054.  setabc(FNPhz, BNPhz, &rnpc);
  1055.  setzeroabc(FNZhz, BNZhz, &rnz);
  1056.  
  1057.  /* Set coefficients of parallel resonators, and amplitude of outputs */
  1058.  setabc(F1hz, B1phz, &r1p);
  1059.  r1p.a *= amp_parF1;
  1060.  setabc(FNPhz, BNPhz, &rnpp);
  1061.  rnpp.a *= amp_parFNP;
  1062.  setabc(F2hz, B2phz, &r2p);
  1063.  r2p.a *= amp_parF2;
  1064.  setabc(F3hz, B3phz, &r3p);
  1065.  r3p.a *= amp_parF3;
  1066.  setabc(F4hz, B4phz, &r4p);
  1067.  r4p.a *= amp_parF4;
  1068.  setabc(F5hz, B5phz, &r5p);
  1069.  r5p.a *= amp_parF5;
  1070.  setabc(F6hz, B6phz, &r6p);
  1071.  r6p.a *= amp_parF6;
  1072.  
  1073.  /* output low-pass filter - resonator with freq 0 and BW = samp_rate
  1074.     Thus 3db point is samp_rate/2 i.e. Nyquist limit.
  1075.     Only 3db down seems rather mild...
  1076.   */
  1077.  
  1078.  setabc(0L, (long) samp_rate, &rout);
  1079.  /* fold overall gain into output resonator */
  1080.  Gain0 = pars->Gain0 - 3;
  1081.  if (Gain0 <= 0)
  1082.   {
  1083.    Gain0 = 57;
  1084.   }
  1085.  rout.a *= DBtoLIN(Gain0);
  1086.  
  1087.  dispt += nspfr;
  1088.  disptcum += nspfr;
  1089. }
  1090.  
  1091. /*
  1092. function PARWAV
  1093.  
  1094. CONVERT FRAME OF PARAMETER DATA TO A WAVEFORM CHUNK
  1095. Synthesize nspfr samples of waveform and store in jwave[].
  1096. */
  1097.  
  1098. void
  1099. parwav(pars, jwave)
  1100. klatt_ptr pars;
  1101. short int *jwave;
  1102. {
  1103.  long ns;
  1104.  float out = 0.0;
  1105.                                 /* Output of cascade branch, also final output  */
  1106.  
  1107.  /* Initialize synthesizer and get specification for current speech
  1108.     frame from host microcomputer */
  1109.  
  1110.  gethost(pars);
  1111.  
  1112.  if (f0_flutter != 0)
  1113.   {
  1114.    time_count++;   /* used for f0 flutter */
  1115.    flutter(pars); /* add f0 flutter */
  1116.   }
  1117.  
  1118.  /* MAIN LOOP, for each output sample of current frame: */
  1119.  
  1120.  for (ns = 0; ns < nspfr; ns++)
  1121.   {
  1122.    int n4;
  1123.    float noise = gen_noise();
  1124.                                 /* Get low-passed random number for aspiration and frication noise */
  1125.    float sourc;     /* Sound source if all-parallel config used     */
  1126.    float glotout; /* Output of glottal sound source               */
  1127.    float par_glotout;
  1128.                                 /* Output of parallelglottal sound sourc        */
  1129.    float voice;     /* Current sample of voicing waveform           */
  1130.    float frics;     /* Frication sound source                       */
  1131.    float aspiration;
  1132.                                 /* Aspiration sound source                      */
  1133.  
  1134.    /* Amplitude modulate noise (reduce noise amplitude during
  1135.       second half of glottal period) if voicing simultaneously present
  1136.    */
  1137.    if (nper > nmod)
  1138.     {
  1139.      noise *= 0.5;
  1140.     }
  1141.  
  1142.    /* Compute frication noise */
  1143.    frics = amp_frica * noise;
  1144.  
  1145.    /* Compute voicing waveform */
  1146.    /* (run glottal source simulation at 4 times normal sample rate to minimize
  1147.     *    quantization noise in period of female voice) */
  1148.  
  1149.    for (n4 = 0; n4 < 4; n4++)
  1150.     {
  1151.      if (glsource == IMPULSIVE)
  1152.       {
  1153.        /* Use impulsive glottal source */
  1154.        voice = impulsive_source(nper);
  1155.       }
  1156.      else
  1157.       {
  1158.        /* Or use a more-natural-shaped source waveform with excitation
  1159.           occurring both upon opening and upon closure, stronest at closure */
  1160.        voice = natural_source(nper);
  1161.       }
  1162.      /* Reset period when counter 'nper' reaches T0 */
  1163.      if (nper >= T0)
  1164.       {
  1165.        nper = 0;
  1166.        pitch_synch_par_reset(ns);
  1167.       }
  1168.  
  1169.      /* Low-pass filter voicing waveform before downsampling from 4*samp_rate */
  1170.      /* to samp_rate samples/sec.  Resonator f=.09*samp_rate, bw=.06*samp_rate    */
  1171.      voice = resonator(&rlp, voice); /* in=voice, out=voice */
  1172.  
  1173.      /* Increment counter that keeps track of 4*samp_rate samples/sec */
  1174.      nper++;
  1175.     }
  1176.  
  1177.    /* Tilt spectrum of voicing source down by soft low-pass filtering, amount
  1178.       of tilt determined by TLTdb
  1179.    */
  1180.    voice = (voice * onemd) + (vlast * decay);
  1181.    vlast = voice;
  1182.  
  1183.    /* Add breathiness during glottal open phase */
  1184.    if (nper < nopen)
  1185.     {
  1186.      /* Amount of breathiness determined by parameter Aturb */
  1187.      /* Use nrand rather than noise because noise is low-passed */
  1188.      voice += amp_breth * nrand;
  1189.     }
  1190.  
  1191.    /* Set amplitude of voicing */
  1192.    glotout = amp_voice * voice;
  1193.  
  1194.    /* Compute aspiration amplitude and add to voicing source */
  1195.    aspiration = amp_aspir * noise;
  1196.    glotout += aspiration;
  1197.  
  1198.    par_glotout = glotout;
  1199.  
  1200.    if (synthesis_model != ALL_PARALLEL)
  1201.     {
  1202.      /* Cascade vocal tract, excited by laryngeal sources.
  1203.         Nasal antiresonator, then formants FNP, F5, F4, F3, F2, F1
  1204.      */
  1205.      float rnzout = antiresonator(&rnz, glotout);
  1206.                                                   /* Output of cascade nazal zero resonator       */
  1207.      float casc_next_in = resonator(&rnpc, rnzout);
  1208.                                                     /* in=rnzout, out=rnpc.p1 */
  1209.  
  1210.      if (nfcascade >= 8)
  1211.       {
  1212.        casc_next_in = resonator(&r8c, casc_next_in);
  1213.        /* Do not use unless samrat = 16000 */
  1214.       }
  1215.  
  1216.      if (nfcascade >= 7)
  1217.       {
  1218.        casc_next_in = resonator(&r7c, casc_next_in);
  1219.        /* Do not use unless samrat = 16000 */
  1220.       }
  1221.  
  1222.      if (nfcascade >= 6)
  1223.       {
  1224.        casc_next_in = resonator(&r6c, casc_next_in);
  1225.        /* Do not use unless long vocal tract or samrat increased */
  1226.       }
  1227.  
  1228.      if (nfcascade >= 5)
  1229.       casc_next_in = resonator(&r5c, casc_next_in);
  1230.  
  1231.      if (nfcascade >= 4)
  1232.       casc_next_in = resonator(&r4c, casc_next_in);
  1233.  
  1234.      if (nfcascade >= 3)
  1235.       casc_next_in = resonator(&r3c, casc_next_in);
  1236.  
  1237.      if (nfcascade >= 2)
  1238.       casc_next_in = resonator(&r2c, casc_next_in);
  1239.  
  1240.      if (nfcascade >= 1)
  1241.       out = resonator(&r1c, casc_next_in);
  1242.      else
  1243.       out = 0.0;
  1244.      /* Excite parallel F1 and FNP by voicing waveform */
  1245.      /* Source is voicing plus aspiration */
  1246.      /* Add in phase, boost lows for nasalized */
  1247.      out += (resonator(&rnpp, par_glotout) + resonator(&r1p, par_glotout));
  1248.     }
  1249.    else
  1250.     {
  1251.      par_glotout = antiresonator(&rnz, par_glotout);
  1252.      par_glotout = resonator(&rnpc, par_glotout);
  1253.      out = resonator(&r1p, par_glotout);
  1254.     }
  1255.  
  1256.    /* Sound sourc for other parallel resonators is frication plus */
  1257.    /* first difference of voicing waveform */
  1258.  
  1259.    sourc = frics + (par_glotout - glotlast);
  1260.    glotlast = par_glotout;
  1261.  
  1262.    /* Standard parallel vocal tract
  1263.       Formants F6,F5,F4,F3,F2, outputs added with alternating sign
  1264.    */
  1265.    out = resonator(&r6p, sourc) - out; /* in=sourc, out=r6p.p1 */
  1266.    out = resonator(&r5p, sourc) - out; /* in=sourc, out=r5p.p1 */
  1267.    out = resonator(&r4p, sourc) - out; /* in=sourc, out=r4p.p1 */
  1268.    out = resonator(&r3p, sourc) - out; /* in=sourc, out=r3p.p1 */
  1269.    out = resonator(&r2p, sourc) - out; /* in=sourc, out=r2p.p1 */
  1270.  
  1271.    out = amp_bypas * sourc - out;
  1272.  
  1273.    out = resonator(&rout, out);
  1274.    /* Convert back to integer */
  1275.    *jwave++ = clip(out);
  1276.   }
  1277. }
  1278.