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