/* PointProcess_and_Sound.c */
/* Paul Boersma, March 9, 2001 */

#include "PointProcess_and_Sound.h"

Sound PointProcess_to_Sound
	(PointProcess me, double samplingFrequency,
	 double adaptFactor, double adaptTime, long interpolationDepth)
{
	Sound thee;
	long it;
	long sound_nt = 1 + floor ((my xmax - my xmin) * samplingFrequency);   /* >= 1 */
	double dt = 1.0 / samplingFrequency;
	double tmid = (my xmin + my xmax) / 2;
	double t1 = tmid - 0.5 * (sound_nt - 1) * dt;
	float *sound;
	thee = Sound_create (my xmin, my xmax, sound_nt, dt, t1);
	if (! thee) return NULL;
	sound = thy z [1];
	for (it = 1; it <= my nt; it ++) {
		double t = my t [it], amplitude = 0.9, angle, halfampsinangle;
		long mid = Sampled_xToNearestIndex (thee, t), j;
		long begin = mid - interpolationDepth, end = mid + interpolationDepth;
		if (it <= 2 || my t [it - 2] < my t [it] - adaptTime) {
			amplitude *= adaptFactor;
			if (it == 1 || my t [it - 1] < my t [it] - adaptTime)
				amplitude *= adaptFactor;
		}
		if (begin < 1) begin = 1;
		if (end > thy nx) end = thy nx;
		angle = NUMpi * (Sampled_indexToX (thee, begin) - t) / thy dx;
		halfampsinangle = 0.5 * amplitude * sin (angle);
		for (j = begin; j <= end; j ++) {
			if (fabs (angle) < 1e-6)
				sound [j] += amplitude;
			else if (angle < 0.0)
				sound [j] += halfampsinangle *
					(1 + cos (angle / (mid - begin + 1))) / angle;
			else
				sound [j] += halfampsinangle *
					(1 + cos (angle / (end - mid + 1))) / angle;
			angle += NUMpi;
			halfampsinangle = - halfampsinangle;
		}
	}
	return thee;
}

int PointProcess_playPart (PointProcess me, double tmin, double tmax) {
	Sound sound = PointProcess_to_Sound (me, 22050, 0.7, 0.05, 30);
	if (! sound) return 0;
	Sound_playPart (sound, tmin, tmax, NULL, NULL);
	forget (sound);
	return 1;
}

int PointProcess_play (PointProcess me) {
	return PointProcess_playPart (me, my xmin, my xmax);
}

int PointProcess_hum (PointProcess me, double tmin, double tmax) {
	static float formant [1 + 6] =
		{ 0, 600, 1400, 2400, 3400, 4500, 5500 };
	static float bandwidth [1 + 6] =
		{ 0, 50, 100, 200, 300, 400, 500 };
	Sound sound = PointProcess_to_Sound (me, 22050, 0.7, 0.05, 30);
	if (! sound) return 0;
	if (! Sound_filterWithFormants (sound, tmin, tmax, 6, formant, bandwidth)) return 0;
	Sound_playPart (sound, tmin, tmax, NULL, NULL);
	forget (sound);
	return 1;
}

Sound PointProcess_to_Sound_hum (PointProcess me) {
	static float formant [1 + 6] =
		{ 0, 600, 1400, 2400, 3400, 4500, 5500 };
	static float bandwidth [1 + 6] =
		{ 0, 50, 100, 200, 300, 400, 500 };
	Sound sound = PointProcess_to_Sound (me, 22050, 0.7, 0.05, 30);
	if (! sound) return 0;
	if (! Sound_filterWithFormants (sound, my xmin, my xmax, 6, formant, bandwidth)) return 0;
	return sound;
}

static double Sound_getRms (Sound me, double tmin, double tmax) {
	long imin, imax, numberOfSamples, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfSamples = Matrix_getWindowSamplesX (me, tmin, tmax, & imin, & imax);
	if (numberOfSamples < 1) return NUMundefined;
	for (i = imin; i <= imax; i ++)
		sum += my z [1] [i] * my z [1] [i];
	sum /= numberOfSamples;
	return sqrt (sum);
}

double PointProcess_Sound_getShimmer (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double numerator = 0.0, denominator = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 3) return NUMundefined;
	for (i = imin + 2; i < imax; i ++) {
		double p1 = my t [i - 1] - my t [i - 2], p2 = my t [i] - my t [i - 1], p3 = my t [i + 1] - my t [i];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax && p3 >= pmin && p3 <= pmax) {
			double e1 = Sound_getRms (thee, my t [i - 2], my t [i - 1]);
			double e2 = Sound_getRms (thee, my t [i - 1], my t [i ]);
			double e3 = Sound_getRms (thee, my t [i ], my t [i + 1]);
			if (NUMdefined (e1) && NUMdefined (e2) && NUMdefined (e3)) {
				numerator += fabs (2 * e2 - e1 - e3);
				denominator += e2;
			}
		}
	}
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}

double PointProcess_Sound_getShimmer_vs (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double sum = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 2) return NUMundefined;
	numberOfPeriods = 0;
	for (i = imin + 1; i < imax; i ++) {
		double p1 = my t [i] - my t [i - 1], p2 = my t [i + 1] - my t [i];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax) {
			double e1 = Sound_getRms (thee, my t [i - 1], my t [i]);
			double e2 = Sound_getRms (thee, my t [i ], my t [i + 1]);
			if (NUMdefined (e1) && NUMdefined (e2) && e1 > 0.0 && e2 > 0.0) {
				sum += fabs (log10 (e1 / e2));
				numberOfPeriods ++;
			}
		}
	}
	if (numberOfPeriods < 1) return NUMundefined;
	return 20.0 * sum / numberOfPeriods;
}

double PointProcess_Sound_getShimmer_apq (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	long imin, imax, numberOfPeriods, i;
	double numerator = 0.0, denominator = 0.0;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	numberOfPeriods = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax) - 1;
	if (numberOfPeriods < 11) return NUMundefined;
	numberOfPeriods = 0;
	for (i = imin + 5; i < imax - 5; i ++) {
		double p1 = my t [i - 4] - my t [i - 5], p2 = my t [i - 3] - my t [i - 4], p3 = my t [i - 2] - my t [i - 3],
			p4 = my t [i - 1] - my t [i - 2], p5 = my t [i] - my t [i - 1], p6 = my t [i + 1] - my t [i],
			p7 = my t [i + 2] - my t [i + 1], p8 = my t [i + 3] - my t [i + 2], p9 = my t [i + 4] - my t [i + 3],
			p10 = my t [i + 5] - my t [i + 4], p11 = my t [i + 6] - my t [i + 5];
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax
			&& p3 >= pmin && p3 <= pmax && p4 >= pmin && p4 <= pmax && p5 >= pmin && p5 <= pmax
			&& p6 >= pmin && p6 <= pmax && p7 >= pmin && p7 <= pmax && p8 >= pmin && p8 <= pmax
			&& p9 >= pmin && p9 <= pmax && p10 >= pmin && p10 <= pmax && p11 >= pmin && p11 <= pmax)
		{
			double e1 = Sound_getRms (thee, my t [i - 5], my t [i - 4]);
			double e2 = Sound_getRms (thee, my t [i - 4], my t [i - 3]);
			double e3 = Sound_getRms (thee, my t [i - 3], my t [i - 2]);
			double e4 = Sound_getRms (thee, my t [i - 2], my t [i - 1]);
			double e5 = Sound_getRms (thee, my t [i - 1], my t [i]);
			double e6 = Sound_getRms (thee, my t [i], my t [i + 1]);
			double e7 = Sound_getRms (thee, my t [i + 1], my t [i + 2]);
			double e8 = Sound_getRms (thee, my t [i + 2], my t [i + 3]);
			double e9 = Sound_getRms (thee, my t [i + 3], my t [i + 4]);
			double e10 = Sound_getRms (thee, my t [i + 4], my t [i + 5]);
			double e11 = Sound_getRms (thee, my t [i + 5], my t [i + 6]);
			if (NUMdefined (e1) && NUMdefined (e2) && NUMdefined (e3) && NUMdefined (e4)
				&& NUMdefined (e5) && NUMdefined (e6) && NUMdefined (e7) && NUMdefined (e8)
				&& NUMdefined (e9) && NUMdefined (e10) && NUMdefined (e11))
			{
				numerator += fabs (e6 - (e1 + e2 + e3 + e4 + e5 + e6 + e7 + e8 + e9 + e10 + e11) / 11.0);
				numberOfPeriods ++;
			}
		}
	}
	if (numberOfPeriods < 1) return NUMundefined;
	numerator /= numberOfPeriods;
	numberOfPeriods = 0;
	for (i = imin; i < imax; i ++) {
		double period = my t [i + 1] - my t [i];
		if (pmin == pmax || period >= pmin && period <= pmax) {
			double amplitude = Sound_getRms (thee, my t [i], my t [i + 1]);
			if (NUMdefined (amplitude)) {
				denominator += amplitude;
				numberOfPeriods ++;
			}
		}
	}
	denominator /= numberOfPeriods;
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}

/* End of file PointProcess_and_Sound.c */
