Drawing trend line respects borders

LPC/CPP.cpp

@@ -30,7 +30,7 @@ static autoVEC CPP_getSoundingValues (CPP me, double tmin, double tmax) {
 	autoVEC soundingValues = raw_VEC (numberOfFrames);
 	integer numberOfSoundingFrames = 0;
 	for (integer iframe = imin; iframe <= imax; iframe ++)
-		if (my z [1] [iframe] > 0.0)
+		if (my voiced [iframe])
 			soundingValues [++ numberOfSoundingFrames] = my z [1] [iframe];
 	if (numberOfSoundingFrames < 1)
 		return autoVEC();
@@ -73,6 +73,7 @@ autoCPP CPP_create (double tmin, double tmax, integer nt, double dt, double t1)
 	try {
 		autoCPP me = Thing_new (CPP);
 		Matrix_init (me.get(), tmin, tmax, nt, dt, t1, 1.0, 1.0, 1, 1.0, 1.0);
+		my voiced = zero_BOOLVEC (nt);
 		return me;
 	} catch (MelderError) {
 		Melder_throw (U"CPP not created.");
@@ -113,19 +114,18 @@ double CPP_getQuantile (CPP me, double quantile) {
 	return NUMquantile (soundingValues.get(), quantile);
 }
 
-autoCPP CPP_and_Pitch_to_CPP_markUnvoiced (CPP me, Pitch thee) {
+autoCPP CPP_and_Pitch_to_CPP_markVoicing (CPP me, Pitch thee) {
 	Melder_require (my xmin == thy xmin && my xmax == thy xmax,
 		U"The domains of the CPP and the Pitch should be equal.");
 	try {
 		autoCPP him = Data_copy (me);
 		for (integer i = 1; i <= my nx; i ++) {
-			const double time = my x1 + (i - 1) * my dx;
-			if (! Pitch_isVoiced_t (thee, time))
-				his z [1] [i] = 0.0;
+			const double time = Sampled_indexToX (me, i);
+			my voiced [i] = ( Pitch_isVoiced_t (thee, time) ? true : false );
 		}
 		return him;
 	} catch (MelderError) {
-		Melder_throw (me, U"cannot convert to object with marked unvoiced frames.");
+		Melder_throw (me, U"cannot convert to object with marked (un)voiced frames.");
 	}
 }
 

LPC/CPP.h

@@ -22,10 +22,14 @@
 #include "Pitch.h"
 
 Thing_define (CPP, Vector) {
+
+	autoBOOLVEC voiced;   // from Sound
+
 	void v1_info ()
 		override;
 	int v_domainQuantity () const
 		override { return MelderQuantity_TIME_SECONDS; }
+
 };
 
 /* 
@@ -47,8 +51,9 @@ Thing_define (CPP, Vector) {
 		ny = 1
 		z [1] [1..nt]
 		TODO
-			The Cepstral Prominence Peak value, a real number between 0 dB and > 100 dB:
-			0 dB means no periodicity at all
+			The Cepstral Prominence Peak value, a real number between 0 dB and > 100 dB if background
+			estimation was done by a straight line. 
+			For exponential decay it might incidentally be negative.
 			Normal values for speech are between 10 dB for non-voiced parts and 60 dB for monotone perfectly voiced
 */
 
@@ -71,6 +76,6 @@ autoMatrix CPP_to_Matrix (CPP me);
 
 autoCPP Matrix_to_CPP (Matrix me);
 
-autoCPP CPP_and_Pitch_to_CPP_markUnvoiced (CPP me, Pitch thee);
+autoCPP CPP_and_Pitch_to_CPP_markVoicing (CPP me, Pitch thee);
 
 #endif // _CPP_h_

LPC/PowerCepstrum.cpp

@@ -260,7 +260,7 @@ void structPowerCepstrumWorkspace :: initWorkspace (constPowerCepstrum me, doubl
 {
 	powerCepstrum = me;
 	Function_unidirectionalAutowindow (powerCepstrum, & qminFit, & qmaxFit);
-	Melder_require (qminFit >= powerCepstrum ->xmin && qmaxFit <= powerCepstrum ->xmax,
+	Melder_require (qminFit >= powerCepstrum -> xmin && qmaxFit <= powerCepstrum -> xmax,
 		U"Your quefrency range is outside the domain.");
 	(void) Matrix_getWindowSamplesX (powerCepstrum, qminFit, qmaxFit, & imin, & imax);
 	Melder_clipLeft (2_integer, & imin); // never use q=0 in fitting
@@ -276,7 +276,7 @@ void structPowerCepstrumWorkspace :: initWorkspace (constPowerCepstrum me, doubl
 			powerCepstrum -> ymin, powerCepstrum -> ymax, powerCepstrum -> ny, powerCepstrum -> dy, powerCepstrum ->y1);
 	maximumNumberOfRhamonics = 15;
 	rhamonics = raw_MAT (maximumNumberOfRhamonics, 5_integer);
-	slopeSelector = SlopeSelector_create (x.get(), y.get()); // only reference to the x and y vectors
+	slopeSelector = SlopeSelector_create (x.get(), y.get()); // argumenys are references to the x and y vectors
 	newData (powerCepstrum); // new xp and yp reference
 }
 
@@ -333,14 +333,22 @@ void PowerCepstrum_drawTrendLine (PowerCepstrum me, Graphics g, double qmin, dou
 	double qstartFit, double qendFit, kCepstrum_trendType lineType, kCepstrum_trendFit fitMethod)
 {
 	Function_unidirectionalAutowindow (me, & qmin, & qmax);
-	Function_intersectRangeWithDomain (me, & qstartFit, & qendFit);
+	Melder_require (Function_intersectRangeWithDomain (me, & qstartFit, & qendFit),
+		U"The fit region should be within the domain.");
+	const double qminLine = std::max (std::max (qmin, qstartFit), 0.1 * my dx);
+	const integer i1Line = Sampled_xToHighIndex (me, qminLine);
+	const double qmaxLine = std::min (qmax, qendFit);
+	const integer i2Line = Sampled_xToLowIndex (me, qmaxLine);
+	const integer numberOfQPoints = i2Line - i1Line;
+	if (numberOfQPoints < 1)
+		return; // nothing to do, interval too short
 	autoPowerCepstrumWorkspace workspace = PowerCepstrumWorkspace_create (me, qstartFit, qendFit, lineType, fitMethod);
 	workspace -> getSlopeAndIntercept ();
 
 	if (dBminimum >= dBmaximum) {   // autoscaling
 		MelderExtremaWithInit extrema_db;
 		for (integer i = workspace -> imin; i <= workspace -> imax; i ++)
-			extrema_db.update (workspace -> y[i]);
+			extrema_db.update (workspace -> y [i]);
 		dBmaximum = extrema_db.max;
 		dBminimum = extrema_db.min;
 	}
@@ -355,30 +363,24 @@ void PowerCepstrum_drawTrendLine (PowerCepstrum me, Graphics g, double qmin, dou
 	Graphics_setLineWidth (g, 2);
 	if (lineType == kCepstrum_trendType::EXPONENTIAL_DECAY ) {
 		integer n = 500;
-		const double dq = (qendFit - qstartFit) / (n + 1);
-		const double q1 = qstartFit;
-		if (qstartFit <= 0.0) {
-			qstartFit = 0.1 * dq;   // some small offset to avoid log(0)
-			n --;
-		}
+		const double dq = (qmaxLine - qminLine) / (n + 1);
 		autoVEC y = raw_VEC (n);
-
 		for (integer i = 1; i <= n; i ++) {
-			const double q = q1 + (i - 1) * dq;
+			const double q = qminLine + (i - 1) * dq;
 			y [i] = slope * log (q) + intercept;
 		}
-		Graphics_function (g, y.asArgumentToFunctionThatExpectsOneBasedArray(), 1, n, qstartFit, qendFit);
+		Graphics_function (g, y.asArgumentToFunctionThatExpectsOneBasedArray(), 1, n, qminLine, qmaxLine);
 	} else {
-		const double y1 = slope * qstartFit + intercept;
-		const double y2 = slope * qendFit + intercept;
+		const double y1 = slope * qminLine + intercept;
+		const double y2 = slope * qmaxLine + intercept;
 		if (y1 >= dBminimum && y2 >= dBminimum) {
-			Graphics_line (g, qstartFit, y1, qendFit, y2);
+			Graphics_line (g, qminLine, y1, qmaxLine, y2);
 		} else if (y1 < dBminimum) {
-			qstartFit = (dBminimum - intercept) / slope;
-			Graphics_line (g, qstartFit, dBminimum, qendFit, y2);
+			const double qstart = (dBminimum - intercept) / slope;
+			Graphics_line (g, qstart, dBminimum, qmaxLine, y2);
 		} else if (y2 < dBminimum) {
-			qendFit = (dBminimum - intercept) / slope;
-			Graphics_line (g, qstartFit, y1, qendFit, dBminimum);
+			const double qend = (dBminimum - intercept) / slope;
+			Graphics_line (g, qminLine, y1, qend, dBminimum);
 		} else {
 			// don't draw anything below lower limit
 		}
@@ -424,7 +426,7 @@ static void PowerCepstrum_smooth_inplaceRectangular (mutablePowerCepstrum me, do
 			for (integer k = 1; k <= numberOfIterations; k ++) {
 				for (integer isamp = 1; isamp <= my nx; isamp ++) {
 					const double xmid = Sampled_indexToX (me, isamp);
-					qout [isamp] = Sampled_getMean (me, xmid - halfWindwow, xmid + halfWindwow, 1, 0, true);
+					qout [isamp] = Sampled_getMean (me, xmid - halfWindwow, xmid + halfWindwow, 1, 0, true); // smooth underlying values not dB's'
 				}
 				my z.row (1)  <<=  qout.all();
 			}

LPC/praat_LPC_init.cpp

@@ -163,7 +163,7 @@ DO
 
 DIRECT (CONVERT_ONE_AND_ONE_TO_ONE__CPP_Pitch_markUnvoiced) {
 	CONVERT_ONE_AND_ONE_TO_ONE (CPP, Pitch)
-		autoCPP result = CPP_and_Pitch_to_CPP_markUnvoiced (me, you);
+		autoCPP result = CPP_and_Pitch_to_CPP_markVoicing (me, you);
 	CONVERT_ONE_AND_ONE_TO_ONE_END (my name.get(), U"_uv")
 }
 
@@ -1808,7 +1808,7 @@ void praat_uvafon_LPC_init () {
 		praat_addAction1 (classCPP, 0, U"Formula...",
 			nullptr, 0, MODIFY_CPP_formula);
 
-	praat_addAction2 (classCPP, 1, classPitch, 1, U"Mark unvoiced frames", nullptr, 0, 
+	praat_addAction2 (classCPP, 1, classPitch, 1, U"Mark unvoiced CPP frames", nullptr, 0, 
 			CONVERT_ONE_AND_ONE_TO_ONE__CPP_Pitch_markUnvoiced);
 
 	praat_addAction1 (classCepstrumc, 0, U"Analyse", nullptr, 0, nullptr);

dwsys/Graphics_extensions.cpp

@@ -263,77 +263,4 @@ bool LineSegmentClipper_clip (LineSegmentClipper me, double& x1, double& y1, dou
 	return my clip (x1, y1, x2, y2);
 }
 
-struct structLineClipper {
-private:
-	typedef struct structHPoint {
-		double x, y, w;
-	} *HPoint;
-	/*
-		The four corners of the rectangle in homogeneous coordinates (x, y, w)
-		Order BottomLeft, BottomRight, TopRight, TopLeft
-		The edges from c[1] to c[2], c[2] to c[3], c[3] to c[4], c[4] to c[1].
-	*/
-	autovector<structHPoint> corners;
-	autovector<structHPoint> edges;
-	autoINTVEC tab1, tab2;
-	struct structHPoint from = {0.0, 0.0, 1.0}, to = {0.0, 0.0, 1.0};
-
-	void cross (HPoint a, HPoint b, HPoint p) {  // p = a x b
-		p -> x = a -> y * b -> w - b -> y * a -> w;
-		p -> y = b -> x * a -> w - a -> x * b -> w;
-		p -> w = a -> x * b -> y - b -> x * a -> y;
-	}
-
-	double in (HPoint a, HPoint b) {
-		return  a -> x * b -> x + a -> y * b -> y + a -> w * b -> w;
-	}
-
-	bool clipLine () {
-		struct structHPoint p;
-		cross (& from, & to, & p);
-		int c = 0;
-		for (int i = 1; i <= 4; i ++)
-			if (in (& p, & corners [i]) >= 0)
-				c |= 1 << (i - 1);
-		if (c == 0 || c == 15) // case 6, Fig 5:
-			return false;
-		Melder_assert (c % 5 != 0);
-		const int i = tab1 [c + 1], j = tab2 [c + 1]; // count from 1
-		cross (& p, & edges [i], & from);
-		cross (& p, & edges [j], & to);
-		return true;
-	}
-
-public:
-
-	void init (double xL, double xR, double yB, double yT) {
-		Melder_assert (xL < xR && yB < yT);
-		corners = {{xL, yB, 1.0}, {xR, yB, 1.0}, {xR, yT, 1.0}, {xL, yT, 1.0}};
-		edges = newvectorzero<structHPoint> (4);
-		for (integer i = 1; i <= 4; i ++) {
-			const integer ip1 = i % 4 + 1;
-			cross (& corners [i], & corners[ip1], & edges[i]);
-		}
-		/*
-			Table 1 in Skala (2012)
-		*/
-		//tab1 = {-1, 3, 0, 3, 1, -2, 0, 3, 2, 2, -2, 2, 1, 1, 0, -1};
-		tab1 = {-1, 4, 1, 4, 2, -2, 1, 4, 3, 3, -2, 3, 2, 2, 1, -1}; // count from 1
-		//tab2 = {-1, 0, 1, 1, 2, -2, 2, 2, 3, 0, -2, 1, 3, 0, 3, -1};
-		tab2 = {-1, 1, 2, 2, 3, -2, 3, 3, 4, 1, -2, 2, 4, 1, 4, -1}; // count from 1
-	}
-
-	bool clip (double& x1, double& y1, double& x2, double& y2) {
-		bool result = false;
-		from.x = x1; from.y = y1; from.w = 1.0;
-		to.x = x2; to.y = y2; to.w = 1.0;
-		if (clipLine ()) {
-			x1 = from.x / from.w; y1 = from.y / from.w;
-			x2 = to.x / to.w;     y2 = to.y / to.w;
-			result = true;
-		}
-		return result;
-	}
-};
-
 /* End of file Graphics_extensions.cpp */

dwsys/Graphics_extensions.h

@@ -125,7 +125,7 @@ struct structLineSegmentClipper {
 		for (int i = 1; i <= 4; i ++)
 			if (in (& p, & corners [i]) >= 0)
 				c |= 1 << (i - 1);
-		c += 1;// we count from 1
+		c += 1;	// we count from 1 instead of 0!
 		if (c == 1 || c == 16) // case 6, Fig 5:
 			return false;
 		Melder_assert (c % 5 != 1);
@@ -148,6 +148,24 @@ struct structLineSegmentClipper {
 		}
 		return true;
 	}
+
+	void init_common () {
+		edges = newvectorzero<structHPoint> (4);
+		for (integer i = 1; i <= 4; i ++) {
+			const integer ip1 = i % 4 + 1;
+			cross (& corners [i], & corners[ip1], & edges[i]);
+		}
+		/*
+			Table 1 in Skala (2005) for indices 0,1,2,3
+			tab1 = {-1, 0, 0, 1, 1, -2, 0, 2, 2, 0, -2, 1, 1, 0, 0, -1};
+			tab2 = {-1, 3, 1, 3, 2, -2, 2, 3, 3, 2, -2, 2, 3, 1, 3, -1};
+			Our indices are 1,2,3,4 so we add 1 to these tables.
+		*/
+		tab1 = {-1, 1, 1, 2, 2, -2, 1, 3, 3, 1, -2, 2, 2, 1, 1, -1}; // count from 1
+		tab2 = {-1, 4, 2, 4, 3, -2, 3, 4, 4, 3, -2, 3, 4, 2, 4, -1};  // count from 1
+		mask = {-1, 4, 4, 2, 2, -2, 4, 8, 8, 4, -2, 2, 2, 4, 4, -1};
+	}
+
 public:
 
 	/*
@@ -164,29 +182,56 @@ struct structLineSegmentClipper {
 			std::swap (y1, y2);
 		//            BL             BR               TR            TL
 		corners = {{x1, y1, 1.0}, {x2, y1, 1.0}, {x2, y2, 1.0}, {x1, y2, 1.0}};
-		edges = newvectorzero<structHPoint> (4);
+		init_common ();
+	}
+
+	/*
+		For a clipping window that consists of 4 points that define a convex quadrilateral.
+		Precondition:
+			The four points need to form a counter clockwise convex quadrilateral
+			No three points on a line.
+	*/
+	void init4 (VEC const& x, VEC const& y) {
+		Melder_assert (x.size == 4 && y.size == 4);
+		/*
+			Check convexity.
+		*/
+		integer cpm1;
 		for (integer i = 1; i <= 4; i ++) {
+			const integer im1 = (i + 2) % 4 + 1;
 			const integer ip1 = i % 4 + 1;
-			cross (& corners [i], & corners[ip1], & edges[i]);
+			const double cp = x [im1] * (y [ip1] - y [i]) - x [i] * (y [ip1] - y [im1]) + x [ip1] * (y [i] - y [im1]);
+			if (i > 1) {
+				Melder_assert ((cp > 0 && cpm1 > 0) || (cp < 0 && cpm1 < 0));
+			}
+			cpm1 = cp;
 		}
-		/*
-			Table 1 in Skala (2005) for indices 0,1,2,3
-			tab1 = {-1, 0, 0, 1, 1, -2, 0, 2, 2, 0, -2, 1, 1, 0, 0, -1};
-			tab2 = {-1, 3, 1, 3, 2, -2, 2, 3, 3, 2, -2, 2, 3, 1, 3, -1};
-			Our indices are 1,2,3,4 so we add 1 to these tables.
-		*/
-		tab1 = {-1, 1, 1, 2, 2, -2, 1, 3, 3, 1, -2, 2, 2, 1, 1, -1}; // count from 1
-		tab2 = {-1, 4, 2, 4, 3, -2, 3, 4, 4, 3, -2, 3, 4, 2, 4, -1};  // count from 1
-		mask = {-1, 4, 4, 2, 2, -2, 4, 8, 8, 4, -2, 2, 2, 4, 4, -1};
+		for (integer i = 1; i <= 4; i ++) {
+			corners [i] = {x [i], y [i], 1.0};
+		}
+		init_common ();
 	}
 
+	/*
+		Check if line from (x1,y1) to (x2,y2) needs to be clipped.
+		if yes, return true and the line to be drawn;
+		E.g if the clipping quadrilateral is given by the four corners ( 0,0), (1,0), (1,1), (0,1):
+			line from (-0.5,0.5) to ( 0.5, 0.5) we return (0.0,0.5) to (0.5,0.5);
+			line from (-0.5,0.5) to ( 1.5, 0.5) we return (0.0,0.5) to (1.0,0.5);
+			line from ( 1.5,0.5) to (-0.5, 0.5) we return (1.0,0.5) to (0.0,0.5);
+			line from ( 0.5,1.5) to ( 0.5, 0.5) we return (0.5,1.0) to (0.5,0.5);
+						 x1 y1         x2 y2                x1 y1        x2 y2
+		if no, return false and arguments are not changed.
+	*/
 	bool clip (double& x1, double& y1, double& x2, double& y2) {
 		bool result = false;
 		from = { x1, y1, 1.0 };
 		to   = { x2, y2, 1.0 };
 		if (clipLine ()) {
-			x1 = from.x / from.w; y1 = from.y / from.w;
-			x2 = to.x / to.w;     y2 = to.y / to.w;
+			x1 = from.x / from.w;
+			y1 = from.y / from.w;
+			x2 =   to.x / to.w;
+			y2 =   to.y / to.w;
 			result = true;
 		}
 		return result;

dwsys/NUMselect.cpp

@@ -22,7 +22,7 @@
 void timeMedian () {
     try {
         Melder_clearInfo ();
-        autoINTVEC numberOfElements {11_integer, 101_integer, 1001_integer, 10001_integer, 100001_integer};
+        autoINTVEC numberOfElements {11_integer, 101_integer, 501, 1001_integer, 10001_integer, 100001_integer};
 		MelderInfo_writeLine (U"Old: O(n log(n)); sort, NUMquantile(0.5)\n"
 			"New: Alexandrescu (2017) O(n)\n"
 		);

dwsys/SlopeSelector.cpp

@@ -428,7 +428,7 @@ static void oneSpecial () { // gave a wrong slope (e-316)
 	const double slope4 = sls -> slopeQuantile_TheilSen (factor);
 }
 
-static void special50 () { // gets stucck because needs > maximumNumberOfTries
+static void special50 () { // gets stuck because needs > maximumNumberOfTries
 	//const integer n = 50;
 	autoVEC x {0, 0.20408163265306123, 0.40816326530612246, 0.61224489795918369, 0.81632653061224492, 1.0204081632653061, 1.2244897959183674, 1.4285714285714286, 1.6326530612244898, 1.8367346938775511, 2.0408163265306123, 2.2448979591836737, 2.4489795918367347, 2.6530612244897958, 2.8571428571428572, 3.0612244897959187, 3.2653061224489797, 3.4693877551020407, 3.6734693877551021, 3.8775510204081636, 4.0816326530612246, 4.2857142857142856, 4.4897959183673475, 4.6938775510204085, 4.8979591836734695, 5.1020408163265305, 5.3061224489795915, 5.5102040816326534, 5.7142857142857144, 5.9183673469387754, 6.1224489795918373, 6.3265306122448983, 6.5306122448979593, 6.7346938775510203, 6.9387755102040813, 7.1428571428571432, 7.3469387755102042, 7.5510204081632653, 7.7551020408163271, 7.9591836734693882, 8.1632653061224492, 8.3673469387755102, 8.5714285714285712, 8.7755102040816322, 8.979591836734695, 9.183673469387756, 9.387755102040817, 9.591836734693878, 9.795918367346939, 10};
 	autoVEC y {0, 0.92803998653746245, 1.8560799730749249, 2.7841199596123873, 3.7121599461498498, 4.6401999326873122, 5.5682399192247747, 6.4962799057622371, 7.4243198922996996, 8.3523598788371629, 9.2803998653746245, 10.208439851912088, 11.136479838449549, 12.064519824987011, 12.992559811524474, 13.920599798061938, 14.848639784599399, 15.776679771136861, 16.704719757674326, 17.632759744211789, 18.560799730749249, 11.892870479075036, 20.416879703824176, 21.344919690361639, 22.272959676899099, 23.200999663436559, 24.129039649974022, 25.057079636511489, 25.985119623048949, 26.913159609586408, 27.841199596123875, 28.769239582661339, 29.697279569198798, 30.625319555736258, 31.553359542273721, 32.481399528811188, 33.409439515348652, 34.337479501886108, 35.265519488423578, 36.193559474961035, 52.020555528667835, 38.049639448035961, 38.977679434573417, 39.905719421110881, 40.833759407648351, 48.501531333854921, 42.689839380723278, 43.617879367260734, 44.545919353798197, 35.700411620170151};