1 /* VocalTract_to_Spectrum.cpp
2 *
3 * Copyright (C) 1991-2005,2008,2011,2015-2018,2020 Paul Boersma
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 * See the GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /*
20 * pb 1991/02/28 Pascal and Fortran versions
21 * pb 1994/03/23 C version
22 * pb 2002/06/04
23 * pb 2002/07/16 GPL
24 * pb 2008/01/19 double
25 * pb 2011/06/12 C++
26 */
27
28 #include "VocalTract_to_Spectrum.h"
29
30 #define cc 353.0
31 #define rho 1.14
32 #define mu 1.86e-5
33 #define lambda 0.0055
34 #define cp 240.0
35 #define eta 1.4
36 #define shapeFactor 2.0
37
TUBE_transfer(double area[],integer numberOfSections,double sectionLength,double frequency,double * re,double * im,double glottalDamping,bool hasRadiationDamping,bool hasInternalDamping)38 static void TUBE_transfer (double area [], integer numberOfSections, double sectionLength, double frequency,
39 double *re, double *im, /* Output. */
40 double glottalDamping, /* 0 or 0.1 */
41 bool hasRadiationDamping, bool hasInternalDamping)
42 /* (Re,Im) := (air current at lips) / (air current at glottis) */
43 {
44 double omega = 2.0 * NUMpi * frequency;
45 if (hasInternalDamping) {
46 double bareResistance = sqrt (mu * rho / 2.0) * sqrt (omega);
47 double bareConductance = (eta - 1.0) / (rho * cc * cc) *
48 sqrt (lambda / (2.0 * cp * rho)) * sqrt (omega);
49 dcomplex c { glottalDamping * area [1] / (rho * cc), 0.0 };
50 dcomplex d { 1.0, 0.0 };
51 for (integer section = 1; section <= numberOfSections; section ++) {
52 double a = area [section];
53 double perimeter = shapeFactor * 2.0 * sqrt (NUMpi * a);
54 double perimeter_by_a2 = perimeter / (a * a);
55 double inertance = rho / a;
56 double compliance = a / (rho * cc * cc);
57 dcomplex cascade { perimeter_by_a2 * bareResistance, omega * inertance + perimeter_by_a2 * bareResistance };
58 dcomplex parallel { perimeter * bareConductance, omega * compliance };
59 dcomplex gamma = dcomplex_sqrt (dcomplex_mul (cascade, parallel));
60 dcomplex impedance = dcomplex_div (gamma, parallel);
61 gamma = dcomplex_rmul (sectionLength, gamma);
62 dcomplex help = dcomplex_rmul (0.5, dcomplex_exp (gamma));
63 dcomplex help1 = dcomplex_div ({ 0.25, 0.0 }, help);
64 dcomplex sinhg = dcomplex_sub (help, help1);
65 dcomplex coshg = dcomplex_add (help, help1);
66 help = dcomplex_add (dcomplex_mul (c, coshg), dcomplex_div (dcomplex_mul (d, sinhg), impedance));
67 d = dcomplex_add (dcomplex_mul (d, coshg), dcomplex_mul (dcomplex_mul (impedance, c), sinhg));
68 c = help;
69 }
70 if (hasRadiationDamping) {
71 double ka = omega * sqrt (area [numberOfSections] / NUMpi) / cc;
72 double z = rho * cc / area [numberOfSections];
73 double radiationResistance = z * ka * ka / 2.0;
74 double radiationReactance = z * 8.0 * ka / 3.0 / NUMpi;
75 *re = d.real() + c.real() * radiationResistance - c.imag() * radiationReactance;
76 *im = d.imag() + c.imag() * radiationResistance + c.real() * radiationReactance;
77 } else { *re = d.real(); *im = d.imag(); };
78 } else {
79 double c_re, c_im;
80 double angle = omega * sectionLength / cc, cosAngle = cos (angle), sinAngle = sin (angle);
81 /* The parallel conductance of the glottis is due to air escaping */
82 /* through it. It is approximated as 0.1 A / rho C */
83 c_re = glottalDamping * cosAngle;
84 c_im = sinAngle;
85 *re = cosAngle;
86 *im = glottalDamping * sinAngle;
87 for (integer section = 1; section < numberOfSections; section ++) {
88 double k = area [section] / area [section + 1];
89 double reDummy = c_re * k * cosAngle - *im * sinAngle;
90 double imDummy = c_im * k * cosAngle + *re * sinAngle;
91 *re = *re * cosAngle - c_im * k * sinAngle;
92 *im = *im * cosAngle + c_re * k * sinAngle;
93 c_re = reDummy;
94 c_im = imDummy;
95 }
96 /* Radiation impedance at the lips. */
97 /* Energy radiates from the mouth. This loses energy to the formants. */
98 if (hasRadiationDamping) {
99 double ka = omega * sqrt (area [numberOfSections] / NUMpi) / cc;
100 double radiationResistance = ka * ka / 2;
101 double radiationReactance = 8.0 * ka / 3.0 / NUMpi;
102 *re += c_re * radiationResistance - c_im * radiationReactance;
103 *im += c_re * radiationReactance + c_im * radiationResistance;
104 }
105 }
106 /* Divide into 1.0. BUG: this is not division into 1. */
107 {
108 double power = *re * *re + *im * *im;
109 if (power != 0.0) { *re = *re / power; *im = *im / power; }
110 }
111 }
112
VocalTract_to_Spectrum(VocalTract me,integer numberOfFrequencies,double maximumFrequency,double glottalDamping,bool hasRadiationDamping,bool hasInternalDamping)113 autoSpectrum VocalTract_to_Spectrum
114 (VocalTract me, integer numberOfFrequencies, double maximumFrequency,
115 double glottalDamping, bool hasRadiationDamping, bool hasInternalDamping)
116 {
117 try {
118 autoSpectrum thee = Spectrum_create (maximumFrequency, numberOfFrequencies);
119 for (integer ifreq = 1; ifreq <= numberOfFrequencies; ifreq ++) {
120 TUBE_transfer (& my z [1] [0], my nx, my dx,
121 (ifreq - 0.9999) * maximumFrequency / (numberOfFrequencies - 1),
122 & thy z [1] [ifreq], & thy z [2] [ifreq],
123 glottalDamping, hasRadiationDamping, hasInternalDamping);
124 thy z [1] [ifreq] *= 0.02;
125 thy z [2] [ifreq] *= 0.02;
126 }
127 return thee;
128 } catch (MelderError) {
129 Melder_throw (me, U": not converted to Spectrum.");
130 }
131 }
132
133 /* End of file VocalTract_to_Spectrum.cpp */
134