/* PointProcess_and_Sound.c
 *
 * Copyright (C) 1992-2003 Paul Boersma
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 * pb 2001/03/09
 * pb 2002/07/16 GPL
 * pb 2003/04/15 improved handling of edges in Sound_getRms
 */

#include "PointProcess_and_Sound.h"
#include "RealTier.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_getPeak (Sound me, double tmin, double tmax) {
	double minimum, timeOfMinimum, maximum, timeOfMaximum;
	float *y = my z [1];
	long i, imin, imax, sampleOfMinimum, sampleOfMaximum;
	if (Sampled_getWindowSamples (me, tmin, tmax, & imin, & imax) < 3) return NUMundefined;
	maximum = minimum = y [imin];
	sampleOfMaximum = sampleOfMinimum = imin;
	for (i = imin + 1; i <= imax; i ++) {
		if (y [i] < minimum) { minimum = y [i]; sampleOfMinimum = i; }
		if (y [i] > maximum) { maximum = y [i]; sampleOfMaximum = i; }
	}
	timeOfMinimum = my x1 + (sampleOfMinimum - 1) * my dx;
	timeOfMaximum = my x1 + (sampleOfMaximum - 1) * my dx;
	Vector_getMinimumAndX (me, timeOfMinimum - my dx, timeOfMinimum + my dx, NUM_PEAK_INTERPOLATE_SINC700, & minimum, & timeOfMinimum);
	Vector_getMaximumAndX (me, timeOfMaximum - my dx, timeOfMaximum + my dx, NUM_PEAK_INTERPOLATE_SINC700, & maximum, & timeOfMaximum);
	return maximum - minimum;
}
static RealTier PointProcess_Sound_getPeaks (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier him = NULL;
	long i, imin, imax, numberOfPeaks = PointProcess_getWindowPoints (me, tmin, tmax, & imin, & imax);
	if (numberOfPeaks < 3) return NULL;
	him = RealTier_create (my xmin, my xmax);
	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 peak = Sound_getPeak (thee, 0.5 * (my t [i - 1] + my t [i]), 0.5 * (my t [i] + my t [i + 1]));
			if (NUMdefined (peak)) RealTier_addPoint (him, my t [i], peak);
		}
	}
	return him;
}
static double RealTier_getShimmer_local (RealTier me, double pmin, double pmax) {
	long i, numberOfPeaks = 0;
	double numerator = 0.0, denominator = 0.0;
	RealPoint *points = (RealPoint *) my points -> item;
	for (i = 2; i <= my points -> size; i ++) {
		double p = points [i] -> time - points [i - 1] -> time;
		if (pmin == pmax || p >= pmin && p <= pmax) {
			numerator += fabs (points [i] -> value - points [i - 1] -> value);
			numberOfPeaks ++;
		}
	}
	if (numberOfPeaks < 1) return NUMundefined;
	numerator /= numberOfPeaks;
	numberOfPeaks = 0;
	for (i = 1; i < my points -> size; i ++) {
		denominator += points [i] -> value;
		numberOfPeaks ++;
	}
	denominator /= numberOfPeaks;
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}
static double RealTier_getShimmer_local_dB (RealTier me, double pmin, double pmax) {
	long i, numberOfPeaks = 0;
	double result = 0.0;
	RealPoint *points = (RealPoint *) my points -> item;
	for (i = 2; i <= my points -> size; i ++) {
		double p = points [i] -> time - points [i - 1] -> time;
		if (pmin == pmax || p >= pmin && p <= pmax) {
			result += fabs (log10 ((points [i] -> value / points [i - 1] -> value)));
			numberOfPeaks ++;
		}
	}
	if (numberOfPeaks < 1) return NUMundefined;
	result /= numberOfPeaks;
	return 20.0 * result;
}
static double RealTier_getShimmer_apq3 (RealTier me, double pmin, double pmax) {
	long i, numberOfPeaks = 0;
	double numerator = 0.0, denominator = 0.0;
	RealPoint *points = (RealPoint *) my points -> item;
	for (i = 2; i <= my points -> size - 1; i ++) {
		double
			p1 = points [i] -> time - points [i - 1] -> time,
			p2 = points [i + 1] -> time - points [i] -> time;
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax)
		{
			double threePointAverage = (points [i - 1] -> value + points [i] -> value + points [i + 1] -> value) / 3.0;
			numerator += fabs (points [i] -> value - threePointAverage);
			numberOfPeaks ++;
		}
	}
	if (numberOfPeaks < 1) return NUMundefined;
	numerator /= numberOfPeaks;
	numberOfPeaks = 0;
	for (i = 1; i < my points -> size; i ++) {
		denominator += points [i] -> value;
		numberOfPeaks ++;
	}
	denominator /= numberOfPeaks;
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}
static double RealTier_getShimmer_apq5 (RealTier me, double pmin, double pmax) {
	long i, numberOfPeaks = 0;
	double numerator = 0.0, denominator = 0.0;
	RealPoint *points = (RealPoint *) my points -> item;
	for (i = 3; i <= my points -> size - 2; i ++) {
		double
			p1 = points [i - 1] -> time - points [i - 2] -> time,
			p2 = points [i] -> time - points [i - 1] -> time,
			p3 = points [i + 1] -> time - points [i] -> time,
			p4 = points [i + 2] -> time - points [i + 1] -> time;
		if (pmin == pmax || p1 >= pmin && p1 <= pmax && p2 >= pmin && p2 <= pmax
			&& p3 >= pmin && p3 <= pmax && p4 >= pmin && p4 <= pmax)
		{
			double fivePointAverage = (points [i - 2] -> value + points [i - 1] -> value + points [i] -> value +
				points [i + 1] -> value + points [i + 2] -> value) / 5.0;
			numerator += fabs (points [i] -> value - fivePointAverage);
			numberOfPeaks ++;
		}
	}
	if (numberOfPeaks < 1) return NUMundefined;
	numerator /= numberOfPeaks;
	numberOfPeaks = 0;
	for (i = 1; i < my points -> size; i ++) {
		denominator += points [i] -> value;
		numberOfPeaks ++;
	}
	denominator /= numberOfPeaks;
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}
static double RealTier_getShimmer_apq11 (RealTier me, double pmin, double pmax) {
	long i, numberOfPeaks = 0;
	double numerator = 0.0, denominator = 0.0;
	RealPoint *points = (RealPoint *) my points -> item;
	for (i = 6; i <= my points -> size - 5; i ++) {
		double
			p1 = points [i - 4] -> time - points [i - 5] -> time,
			p2 = points [i - 3] -> time - points [i - 4] -> time,
			p3 = points [i - 2] -> time - points [i - 3] -> time,
			p4 = points [i - 1] -> time - points [i - 2] -> time,
			p5 = points [i] -> time - points [i - 1] -> time,
			p6 = points [i + 1] -> time - points [i] -> time,
			p7 = points [i + 2] -> time - points [i + 1] -> time,
			p8 = points [i + 3] -> time - points [i + 2] -> time,
			p9 = points [i + 4] -> time - points [i + 3] -> time,
			p10 = points [i + 5] -> time - points [i + 4] -> time;
		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)
		{
			double elevenPointAverage = (points [i - 5] -> value + points [i - 4] -> value + points [i - 3] -> value +
				points [i - 2] -> value + points [i - 1] -> value + points [i] -> value + points [i + 1] -> value +
				points [i + 2] -> value + points [i + 3] -> value + points [i + 4] -> value + points [i + 5] -> value) / 11.0;
			numerator += fabs (points [i] -> value - elevenPointAverage);
			numberOfPeaks ++;
		}
	}
	if (numberOfPeaks < 1) return NUMundefined;
	numerator /= numberOfPeaks;
	numberOfPeaks = 0;
	for (i = 1; i < my points -> size; i ++) {
		denominator += points [i] -> value;
		numberOfPeaks ++;
	}
	denominator /= numberOfPeaks;
	if (denominator == 0.0) return NUMundefined;
	return numerator / denominator;
}

double PointProcess_Sound_getShimmer_local (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double result;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	result = RealTier_getShimmer_local (peaks, pmin, pmax);
	forget (peaks);
	return result;
}

double PointProcess_Sound_getShimmer_local_dB (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double result;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	result = RealTier_getShimmer_local_dB (peaks, pmin, pmax);
	forget (peaks);
	return result;
}

double PointProcess_Sound_getShimmer_apq3 (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double result;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	result = RealTier_getShimmer_apq3 (peaks, pmin, pmax);
	forget (peaks);
	return result;
}

double PointProcess_Sound_getShimmer_apq5 (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double result;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	result = RealTier_getShimmer_apq5 (peaks, pmin, pmax);
	forget (peaks);
	return result;
}

double PointProcess_Sound_getShimmer_apq11 (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double result;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	result = RealTier_getShimmer_apq11 (peaks, pmin, pmax);
	forget (peaks);
	return result;
}

double PointProcess_Sound_getShimmer_dda (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax) {
	RealTier peaks;
	double apq3;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	apq3 = RealTier_getShimmer_apq3 (peaks, pmin, pmax);
	forget (peaks);
	return NUMdefined (apq3) ? 3.0 * apq3 : NUMundefined;
}

void PointProcess_Sound_getShimmer_multi (PointProcess me, Sound thee, double tmin, double tmax, double pmin, double pmax,
	double *local, double *local_dB, double *apq3, double *apq5, double *apq11, double *dda)
{
	RealTier peaks;
	if (tmax <= tmin) tmin = my xmin, tmax = my xmax;   /* Autowindowing. */
	peaks = PointProcess_Sound_getPeaks (me, thee, tmin, tmax, pmin, pmax);
	if (local) *local = RealTier_getShimmer_local (peaks, pmin, pmax);
	if (local_dB) *local_dB = RealTier_getShimmer_local_dB (peaks, pmin, pmax);
	if (apq3) *apq3 = RealTier_getShimmer_apq3 (peaks, pmin, pmax);
	if (apq5) *apq5 = RealTier_getShimmer_apq5 (peaks, pmin, pmax);
	if (apq11) *apq11 = RealTier_getShimmer_apq11 (peaks, pmin, pmax);
	if (dda) *dda = 3.0 * RealTier_getShimmer_apq3 (peaks, pmin, pmax);
	forget (peaks);
}

/* End of file PointProcess_and_Sound.c */