1 /* PointProcess_and_Sound.cpp 2 * 3 * Copyright (C) 1992-2012,2014-2018,2021 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 #include "PointProcess_and_Sound.h" 20 21 autoSound PointProcess_to_Sound_pulseTrain 22 (PointProcess me, double samplingFrequency, 23 double adaptFactor, double adaptTime, integer interpolationDepth) 24 { 25 try { 26 const integer sound_nt = 1 + Melder_ifloor ((my xmax - my xmin) * samplingFrequency); // >= 1 27 const double dt = 1.0 / samplingFrequency; 28 const double tmid = 0.5 * (my xmin + my xmax); 29 const double t1 = tmid - 0.5 * (sound_nt - 1) * dt; 30 autoSound thee = Sound_create (1, my xmin, my xmax, sound_nt, dt, t1); 31 const VEC sound = thy z.row (1); 32 for (integer it = 1; it <= my nt; it ++) { 33 const double t = my t [it]; 34 double amplitude = 0.9; 35 integer mid = Sampled_xToNearestIndex (thee.get(), t); 36 if (it <= 2 || my t [it - 2] < my t [it] - adaptTime) { 37 amplitude *= adaptFactor; 38 if (it == 1 || my t [it - 1] < my t [it] - adaptTime) 39 amplitude *= adaptFactor; 40 } 41 integer begin = mid - interpolationDepth, end = mid + interpolationDepth; 42 if (begin < 1) begin = 1; ~TargetInstrInfo()43 if (end > thy nx) end = thy nx; 44 double angle = NUMpi * (Sampled_indexToX (thee.get(), begin) - t) / thy dx; 45 double halfampsinangle = 0.5 * amplitude * sin (angle); 46 for (integer j = begin; j <= end; j ++) { getRegClass(const MCInstrDesc & MCID,unsigned OpNum,const TargetRegisterInfo * TRI,const MachineFunction & MF) const47 if (fabs (angle) < 1e-6) 48 sound [j] += amplitude; 49 else if (angle < 0.0) 50 sound [j] += halfampsinangle * 51 (1.0 + cos (angle / (mid - begin + 1))) / angle; 52 else 53 sound [j] += halfampsinangle * 54 (1.0 + cos (angle / (end - mid + 1))) / angle; 55 angle += NUMpi; 56 halfampsinangle = - halfampsinangle; 57 } 58 } 59 return thee; 60 } catch (MelderError) { 61 Melder_throw (me, U": pulse train not synthesized."); 62 } 63 } 64 65 autoSound PointProcess_to_Sound_phonation 66 (PointProcess me, double samplingFrequency, double adaptFactor, double maximumPeriod, insertNoop(MachineBasicBlock & MBB,MachineBasicBlock::iterator MI) const67 double openPhase, double collisionPhase, double power1, double power2) 68 { 69 try { 70 const integer sound_nt = 1 + Melder_ifloor ((my xmax - my xmin) * samplingFrequency); // >= 1 71 const double dt = 1.0 / samplingFrequency; 72 const double tmid = 0.5 * (my xmin + my xmax); 73 const double t1 = tmid - 0.5 * (sound_nt - 1) * dt; 74 const double a = (power1 + power2 + 1.0) / (power2 - power1); 75 autoSound thee = Sound_create (1, my xmin, my xmax, sound_nt, dt, t1); 76 /* 77 Compute "re" by iteration. 78 */ 79 double re = openPhase - collisionPhase; 80 if (collisionPhase <= 0.0) { 81 re = openPhase; 82 } else { 83 const double xmaxFlow = pow (power1 / power2, 1.0 / (power2 - power1)); 84 /* mutable by iteration */ double xleft = xmaxFlow; 85 /* mutable by iteration */ double xright = 1.0; 86 for (int i = 1; i <= 50; i ++) { 87 const double xmid = 0.5 * (xleft + xright); 88 const double gmid = pow (xmid, power1) - pow (xmid, power2); 89 const double gderivmid = power1 * pow (xmid, power1 - 1.0) - power2 * pow (xmid, power2 - 1.0); 90 const double fmid = - gmid / gderivmid; 91 if (fmid > collisionPhase / openPhase) 92 xleft = xmid; 93 else 94 xright = xmid; 95 re = xmid * openPhase; 96 } 97 } 98 /* 99 Cycle through the points. Each will become a period. 100 */ 101 VEC sound = thy z.row (1); 102 for (integer it = 1; it <= my nt; it ++) { 103 const double t = my t [it]; 104 const integer midSample = Sampled_xToNearestIndex (thee.get(), t); 105 /* 106 Determine the period: first look left (because that's where the open phase is), 107 then right. 108 */ 109 double period = undefined; 110 if (it >= 2) { 111 period = my t [it] - my t [it - 1]; 112 if (period > maximumPeriod) 113 period = undefined; 114 } 115 if (isundef (period)) { 116 if (it < my nt) { 117 period = my t [it + 1] - my t [it]; 118 if (period > maximumPeriod) 119 period = undefined; 120 } 121 if (isundef (period)) 122 period = 0.5 * maximumPeriod; // some default value 123 } 124 const double te = re * period; 125 /* 126 Determine the amplitude of this peak. 127 */ 128 const double unadaptedAmplitude = a / (period * openPhase); 129 const double amplitude = unadaptedAmplitude * ( 130 it == 1 || my t [it - 1] < my t [it] - maximumPeriod ? 131 adaptFactor * adaptFactor 132 : it == 2 || my t [it - 2] < my t [it - 1] - maximumPeriod ? 133 adaptFactor 134 : 135 1.0 136 ); 137 /* 138 Fill in the samples to the left of the current point. 139 */ 140 {// scope 141 const integer beginSample = Melder_clippedLeft (1_integer, midSample - Melder_ifloor (te / thy dx)); 142 const integer endSample = Melder_clippedRight (midSample, thy nx); 143 for (integer isamp = beginSample; isamp <= endSample; isamp ++) { 144 const double tsamp = thy x1 + (isamp - 1) * thy dx; 145 const double phase = (tsamp - (t - te)) / (period * openPhase); 146 if (phase > 0.0) 147 sound [isamp] += amplitude * (power1 * pow (phase, power1 - 1.0) - power2 * pow (phase, power2 - 1.0)); 148 } 149 } 150 /* 151 Determine the signal parameters at the current point. 152 */ 153 const double phase = te / (period * openPhase); 154 const double flow = amplitude * (period * openPhase) * (pow (phase, power1) - pow (phase, power2)); 155 /* 156 Fill in the samples to the right of the current point. 157 */ 158 if (flow > 0.0) { 159 const double flowDerivative = amplitude * (power1 * pow (phase, power1 - 1.0) - power2 * pow (phase, power2 - 1.0)); 160 const double ta = - flow / flowDerivative; 161 const double factorPerSample = exp (- thy dx / ta); 162 const integer beginSample = Melder_clippedLeft (1_integer, midSample + 1); 163 const integer endSample = Melder_clippedRight (midSample + Melder_ifloor (20.0 * ta / thy dx), thy nx); 164 double value = flowDerivative * factorPerSample; 165 for (integer isamp = beginSample; isamp <= endSample; isamp ++) { 166 sound [isamp] += value; 167 value *= factorPerSample; 168 } 169 } 170 } 171 Vector_scale (thee.get(), 0.9); 172 return thee; 173 } catch (MelderError) { 174 Melder_throw (me, U": not converted to Sound (phonation)."); 175 } 176 } 177 178 void PointProcess_playPart (PointProcess me, double tmin, double tmax) { 179 try { 180 autoSound sound = PointProcess_to_Sound_pulseTrain (me, 44100.0, 0.7, 0.05, 30); 181 Sound_playPart (sound.get(), tmin, tmax, nullptr, nullptr); 182 } catch (MelderError) { 183 Melder_throw (me, U": not played."); 184 } 185 } 186 187 void PointProcess_play (PointProcess me) { 188 PointProcess_playPart (me, my xmin, my xmax); 189 } 190 191 void PointProcess_hum (PointProcess me, double tmin, double tmax) { 192 static double formant [1 + 6] = { 0, 600.0, 1400.0, 2400.0, 3400.0, 4500.0, 5500.0 }; 193 static double bandwidth [1 + 6] = { 0, 50.0, 100.0, 200.0, 300.0, 400.0, 500.0 }; 194 autoSound sound = PointProcess_to_Sound_pulseTrain (me, 44100, 0.7, 0.05, 30); 195 Sound_filterWithFormants (sound.get(), tmin, tmax, 6, formant, bandwidth); 196 Sound_playPart (sound.get(), tmin, tmax, nullptr, nullptr); 197 } 198 199 autoSound PointProcess_to_Sound_hum (PointProcess me) { 200 static double formant [1 + 6] = { 0, 600.0, 1400.0, 2400.0, 3400.0, 4500.0, 5500.0 }; 201 static double bandwidth [1 + 6] = { 0, 50.0, 100.0, 200.0, 300.0, 400.0, 500.0 }; 202 try { 203 autoSound sound = PointProcess_to_Sound_pulseTrain (me, 44100.0, 0.7, 0.05, 30); 204 Sound_filterWithFormants (sound.get(), my xmin, my xmax, 6, formant, bandwidth); 205 return sound; 206 } catch (MelderError) { 207 Melder_throw (me, U": not converted to Sound (hum)."); 208 } 209 } 210 211 /* End of file PointProcess_and_Sound.cpp */ 212