1 /* FormantTier.cpp 2 * 3 * Copyright (C) 1992-2007,2011,2012,2014-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 2002/07/16 GPL 21 * pb 2006/07/21 made Sound_FormantTier_filter_inplace () accurate for higher numbers of formants 22 * pb 2007/01/27 made compatible with stereo sounds 23 * pb 2007/03/17 domain quantity 24 * pb 2007/10/01 can write as encoding 25 * pb 2011/03/01 moved Formant filtering to FormantGrid (reimplemented) 26 * pb 2011/05/26 C++ 27 */ 28 29 #include "FormantTier.h" 30 #include "AnyTier.h" 31 32 #include "oo_DESTROY.h" 33 #include "FormantTier_def.h" 34 #include "oo_COPY.h" 35 #include "FormantTier_def.h" 36 #include "oo_EQUAL.h" 37 #include "FormantTier_def.h" 38 #include "oo_CAN_WRITE_AS_ENCODING.h" 39 #include "FormantTier_def.h" 40 #include "oo_WRITE_TEXT.h" 41 #include "FormantTier_def.h" 42 #include "oo_READ_TEXT.h" 43 #include "FormantTier_def.h" 44 #include "oo_WRITE_BINARY.h" 45 #include "FormantTier_def.h" 46 #include "oo_READ_BINARY.h" 47 #include "FormantTier_def.h" 48 #include "oo_DESCRIPTION.h" 49 #include "FormantTier_def.h" 50 51 Thing_implement (FormantPoint, Daata, 0); 52 53 autoFormantPoint FormantPoint_create (double time, integer numberOfFormants) { 54 try { 55 autoFormantPoint me = Thing_new (FormantPoint); 56 my number = time; 57 my numberOfFormants = numberOfFormants; 58 my formant = zero_VEC (numberOfFormants); 59 my bandwidth = zero_VEC (numberOfFormants); 60 return me; 61 } catch (MelderError) { 62 Melder_throw (U"Formant point not created."); 63 } 64 } 65 66 Thing_implement (FormantTier, AnyTier, 0); 67 68 autoFormantTier FormantTier_create (double tmin, double tmax) { 69 try { 70 autoFormantTier me = Thing_new (FormantTier); 71 my xmin = tmin; 72 my xmax = tmax; 73 return me; 74 } catch (MelderError) { 75 Melder_throw (U"FormantTier not created."); 76 } 77 } 78 79 double FormantTier_getValueAtTime (FormantTier me, integer iformant, double t) { 80 integer n = my points.size; 81 if (n == 0 || iformant < 1) 82 return undefined; 83 FormantPoint pointRight = my points.at [1]; 84 if (t <= pointRight -> number) { 85 if (iformant > pointRight -> numberOfFormants) 86 return undefined; 87 return pointRight -> formant [iformant]; // constant extrapolation 88 } 89 FormantPoint pointLeft = my points.at [n]; 90 if (t >= pointLeft -> number) { 91 if (iformant > pointLeft -> numberOfFormants) 92 return undefined; 93 return pointLeft -> formant [iformant]; // constant extrapolation 94 } 95 Melder_assert (n >= 2); 96 integer ileft = AnyTier_timeToLowIndex (me->asAnyTier(), t), iright = ileft + 1; 97 Melder_assert (ileft >= 1 && iright <= n); 98 pointLeft = my points.at [ileft]; 99 pointRight = my points.at [iright]; 100 const double tleft = pointLeft -> number; 101 const double fleft = ( iformant > pointLeft -> numberOfFormants ? undefined : pointLeft -> formant [iformant] ); 102 const double tright = pointRight -> number; 103 const double fright = ( iformant > pointRight -> numberOfFormants ? undefined : pointRight -> formant [iformant] ); 104 return isundef (fleft) ? ( isundef (fright) ? undefined : fright ) 105 : isundef (fright) ? fleft 106 : t == tright ? fright // be very accurate 107 : tleft == tright ? 0.5 * (fleft + fright) // unusual, but possible; no preference 108 : fleft + (t - tleft) * (fright - fleft) / (tright - tleft); // linear interpolation 109 } 110 111 double FormantTier_getBandwidthAtTime (FormantTier me, integer iformant, double t) { 112 const integer n = my points.size; 113 if (n == 0 || iformant < 1) 114 return undefined; 115 FormantPoint pointRight = my points.at [1]; 116 if (t <= pointRight -> number) { 117 if (iformant > pointRight -> numberOfFormants) 118 return undefined; 119 return pointRight -> bandwidth [iformant]; // constant extrapolation 120 } 121 FormantPoint pointLeft = my points.at [n]; 122 if (t >= pointLeft -> number) { 123 if (iformant > pointLeft -> numberOfFormants) 124 return undefined; 125 return pointLeft -> bandwidth [iformant]; // constant extrapolation 126 } 127 Melder_assert (n >= 2); 128 const integer ileft = AnyTier_timeToLowIndex (me->asAnyTier(), t), iright = ileft + 1; 129 Melder_assert (ileft >= 1 && iright <= n); 130 pointLeft = my points.at [ileft]; 131 pointRight = my points.at [iright]; 132 const double tleft = pointLeft -> number; 133 const double fleft = iformant > pointLeft -> numberOfFormants ? undefined : pointLeft -> bandwidth [iformant]; 134 const double tright = pointRight -> number; 135 const double fright = iformant > pointRight -> numberOfFormants ? undefined : pointRight -> bandwidth [iformant]; 136 return isundef (fleft) ? ( isundef (fright) ? undefined : fright ) 137 : isundef (fright) ? fleft 138 : t == tright ? fright // be very accurate 139 : tleft == tright ? 0.5 * (fleft + fright) // unusual, but possible; no preference 140 : fleft + (t - tleft) * (fright - fleft) / (tright - tleft); // linear interpolation 141 } 142 143 void FormantTier_speckle (FormantTier me, Graphics g, double tmin, double tmax, double fmax, bool garnish) { 144 Function_unidirectionalAutowindow (me, & tmin, & tmax); 145 Graphics_setWindow (g, tmin, tmax, 0.0, fmax); 146 Graphics_setInner (g); 147 const integer imin = AnyTier_timeToHighIndex (me->asAnyTier(), tmin); 148 const integer imax = AnyTier_timeToLowIndex (me->asAnyTier(), tmax); 149 if (imin > 0) for (integer i = imin; i <= imax; i ++) { 150 FormantPoint point = my points.at [i]; 151 const double t = point -> number; 152 for (integer j = 1; j <= point -> numberOfFormants; j ++) { 153 const double f = point -> formant [j]; 154 if (f <= fmax) 155 Graphics_speckle (g, t, f); 156 } 157 } 158 Graphics_unsetInner (g); 159 if (garnish) { 160 Graphics_drawInnerBox (g); 161 Graphics_textBottom (g, true, U"Time (s)"); 162 Graphics_marksBottom (g, 2, true, true, false); 163 Graphics_marksLeft (g, 2, true, true, false); 164 Graphics_textLeft (g, true, U"Frequency (Hz)"); 165 } 166 } 167 168 autoFormantTier Formant_downto_FormantTier (Formant me) { 169 try { 170 autoFormantTier thee = FormantTier_create (my xmin, my xmax); 171 for (integer i = 1; i <= my nx; i ++) { 172 const Formant_Frame frame = & my frames [i]; 173 /* mutable move */ autoFormantPoint point = FormantPoint_create (Sampled_indexToX (me, i), frame -> numberOfFormants); 174 for (integer j = 1; j <= frame -> numberOfFormants; j ++) { 175 Formant_Formant pair = & frame -> formant [j]; 176 point -> formant [j] = pair -> frequency; 177 point -> bandwidth [j] = pair -> bandwidth; 178 } 179 thy points. addItem_move (point.move()); 180 } 181 return thee; 182 } catch (MelderError) { 183 Melder_throw (me, U": not converted to FormantTier."); 184 } 185 } 186 187 autoFormantTier Formant_PointProcess_to_FormantTier (Formant me, PointProcess pp) { 188 try { 189 const integer maximumNumberOfFormants = Formant_getMaxNumFormants (me); 190 const autoFormantTier temp = Formant_downto_FormantTier (me); 191 /* mutable return */ autoFormantTier thee = FormantTier_create (pp -> xmin, pp -> xmax); 192 for (integer ipoint = 1; ipoint <= pp -> nt; ipoint ++) { 193 const double time = pp -> t [ipoint]; 194 /* mutable move */ autoFormantPoint point = FormantPoint_create (time, maximumNumberOfFormants); 195 /* mutable loop */ integer iformant = 1; 196 for (; iformant <= maximumNumberOfFormants; iformant ++) { 197 const double frequency = FormantTier_getValueAtTime (temp.get(), iformant, time); 198 if (isundef (frequency)) 199 break; 200 point -> formant [iformant] = frequency; 201 const double bandwidth = FormantTier_getBandwidthAtTime (temp.get(), iformant, time); 202 Melder_assert (isdefined (bandwidth)); 203 point -> bandwidth [iformant] = bandwidth; 204 } 205 point -> numberOfFormants = iformant - 1; 206 point -> formant. resize (point -> numberOfFormants); // maintain invariant 207 point -> bandwidth. resize (point -> numberOfFormants); // maintain invariant 208 thy points. addItem_move (point.move()); 209 } 210 return thee; 211 } catch (MelderError) { 212 Melder_throw (me, U" & ", pp, U": not converted to FormantTier."); 213 } 214 } 215 216 integer FormantTier_getMinNumFormants (FormantTier me) { 217 integer minNumFormants = INTEGER_MAX; 218 for (integer ipoint = 1; ipoint <= my points.size; ipoint ++) { 219 const FormantPoint point = my points.at [ipoint]; 220 if (point -> numberOfFormants < minNumFormants) 221 minNumFormants = point -> numberOfFormants; 222 } 223 return minNumFormants; 224 } 225 226 integer FormantTier_getMaxNumFormants (FormantTier me) { 227 integer maxNumFormants = 0; 228 for (integer ipoint = 1; ipoint <= my points.size; ipoint ++) { 229 const FormantPoint point = my points.at [ipoint]; 230 if (point -> numberOfFormants > maxNumFormants) 231 maxNumFormants = point -> numberOfFormants; 232 } 233 return maxNumFormants; 234 } 235 236 autoTableOfReal FormantTier_downto_TableOfReal (FormantTier me, bool includeFormants, bool includeBandwidths) { 237 try { 238 const integer maximumNumberOfFormants = FormantTier_getMaxNumFormants (me); 239 /* mutable return */ autoTableOfReal thee = TableOfReal_create (my points.size, 1 + 240 ( includeFormants ? maximumNumberOfFormants : 0 ) + 241 ( includeBandwidths ? maximumNumberOfFormants : 0 )); 242 TableOfReal_setColumnLabel (thee.get(), 1, U"Time"); 243 for (integer icol = 1, iformant = 1; iformant <= maximumNumberOfFormants; iformant ++) { 244 char32 label [4]; 245 if (includeFormants) { 246 Melder_sprint (label,4, U"F", iformant); 247 TableOfReal_setColumnLabel (thee.get(), ++ icol, label); 248 } 249 if (includeBandwidths) { 250 Melder_sprint (label,4, U"B", iformant); 251 TableOfReal_setColumnLabel (thee.get(), ++ icol, label); 252 } 253 } 254 for (integer ipoint = 1; ipoint <= my points.size; ipoint ++) { 255 const FormantPoint point = my points.at [ipoint]; 256 thy data [ipoint] [1] = point -> number; 257 for (integer icol = 1, iformant = 1; iformant <= maximumNumberOfFormants; iformant ++) { 258 if (includeFormants) 259 thy data [ipoint] [++ icol] = point -> formant [iformant]; 260 if (includeBandwidths) 261 thy data [ipoint] [++ icol] = point -> bandwidth [iformant]; 262 } 263 } 264 return thee; 265 } catch (MelderError) { 266 Melder_throw (me, U": not converted to TableOfReal."); 267 } 268 } 269 270 void Sound_FormantTier_filter_inplace (Sound me, FormantTier formantTier) { 271 const double dt = my dx; 272 const integer maximumNumberOfFormants = FormantTier_getMaxNumFormants (formantTier); 273 if (formantTier -> points.size) for (integer iformant = 1; iformant <= maximumNumberOfFormants; iformant ++) { 274 for (integer isamp = 1; isamp <= my nx; isamp ++) { 275 const double t = my x1 + (isamp - 1) * my dx; 276 /* 277 Compute LP coefficients. 278 */ 279 const double formant = FormantTier_getValueAtTime (formantTier, iformant, t); 280 const double bandwidth = FormantTier_getBandwidthAtTime (formantTier, iformant, t); 281 if (isdefined (formant) && isdefined (bandwidth)) { 282 const double cosomdt = cos (2.0 * NUMpi * formant * dt); 283 const double r = exp (- NUMpi * bandwidth * dt); 284 /* Formants at 0 Hz or the Nyquist are single poles; others are double poles. */ 285 if (fabs (cosomdt) > 0.999999) { // allow for round-off errors 286 /* single pole: D(z) = 1 - r z^-1 */ 287 for (integer channel = 1; channel <= my ny; channel ++) 288 if (isamp > 1) 289 my z [channel] [isamp] += r * my z [channel] [isamp - 1]; 290 } else { 291 /* double pole: D(z) = 1 + p z^-1 + q z^-2 */ 292 const double p = - 2.0 * r * cosomdt; 293 const double q = r * r; 294 for (integer channel = 1; channel <= my ny; channel ++) { 295 if (isamp > 1) 296 my z [channel] [isamp] -= p * my z [channel] [isamp - 1]; 297 if (isamp > 2) 298 my z [channel] [isamp] -= q * my z [channel] [isamp - 2]; 299 } 300 } 301 } 302 } 303 } 304 } 305 306 autoSound Sound_FormantTier_filter (Sound me, FormantTier formantTier) { 307 try { 308 autoSound thee = Data_copy (me); 309 Sound_FormantTier_filter_inplace (thee.get(), formantTier); 310 Vector_scale (thee.get(), 0.99); 311 return thee; 312 } catch (MelderError) { 313 Melder_throw (me, U": not filtered with ", formantTier, U"."); 314 } 315 } 316 317 autoSound Sound_FormantTier_filter_noscale (Sound me, FormantTier formantTier) { 318 try { 319 autoSound thee = Data_copy (me); 320 Sound_FormantTier_filter_inplace (thee.get(), formantTier); 321 return thee; 322 } catch (MelderError) { 323 Melder_throw (me, U": not filtered with ", formantTier, U"."); 324 } 325 } 326 327 /* End of file FormantTier.cpp */ 328