/* Sound_to_Intensity.c */
/* Paul Boersma, July 3, 2001 */

#include "Sound_to_Intensity.h"

static double Sound_smoothAt (Sound me, double midTime, double windowDuration) {
	double halfWindowDuration = 0.5 * windowDuration;
	double leftTime = midTime - halfWindowDuration, rightTime = midTime + halfWindowDuration;
	long i, leftSample = Sampled_xToHighIndex (me, leftTime), rightSample = Sampled_xToLowIndex (me, rightTime);
	float *amplitude = my z [1], sumxw = 0.0, sumw = 0.0;

	if (leftSample < 1) leftSample = 1;
	if (rightSample > my nx) rightSample = my nx;
	for (i = leftSample; i <= rightSample; i ++) {
		double t = Sampled_indexToX (me, i), x = (midTime - t) / halfWindowDuration, root = 1 - x * x;
		double kaiser = root <= 0.0 ? 0.0 : NUMbesselI (0, (2 * NUMpi * NUMpi + 0.5) * sqrt (root));
		sumxw += amplitude [i] * kaiser;
		sumw += kaiser;
	}
	return sumxw / sumw;
}

Intensity Sound_to_Intensity (Sound me, double minimumPitch, double timeStep) {
	Sound up;
	Intensity smooth;
	float *amplitude;
	long i, iframe, numberOfFrames;
	double windowDuration, thyFirstTime;

	/* Step 1: upsample by a factor of two. */
	up = /*Sound_upsample (me);*/ Data_copy (me);
	if (! up) return Melder_errorp ("Intensity analysis not performed.");

	/* Step 2: square all samples (this doubles frequency content). */
	amplitude = up -> z [1];
	for (i = 1; i <= up -> nx; i ++)
		amplitude [i] *= amplitude [i];

	/* Step 3: smoothe and resample. */
	Melder_assert (minimumPitch > 0.0);
	windowDuration = 6.4 / minimumPitch;
	if (timeStep <= 0.0) timeStep = 1 / minimumPitch;   /* Default. */
	if (! Sampled_shortTermAnalysis (up, windowDuration, timeStep, & numberOfFrames, & thyFirstTime)) {
		forget (up);
		return Melder_errorp ("Sound should be at least 6.4 / minimumPitch long.");
	}
	if (! (smooth = Intensity_create (my xmin, my xmax, numberOfFrames, timeStep, thyFirstTime))) {
		forget (up);
		return Melder_errorp ("Intensity analysis not performed.");
	}
	for (iframe = 1; iframe <= numberOfFrames; iframe ++) {
		double midTime = Sampled_indexToX (smooth, iframe);
		double intensity = Sound_smoothAt (up, midTime, windowDuration);
		if (intensity != 0.0) intensity /= 4e-10;
		smooth -> z [1] [iframe] = intensity < 1e-30 ? -300 : 10 * log10 (intensity);
	}

	/* Clean up and return. */
	forget (up);
	return smooth;
}

IntensityTier Sound_to_IntensityTier (Sound me, double minimumPitch, double timeStep) {
	Intensity intensity = Sound_to_Intensity (me, minimumPitch, timeStep);
	IntensityTier thee;
	if (! intensity) return NULL;
	thee = Intensity_downto_IntensityTier (intensity);
	forget (intensity);
	return thee;
}

/* End of file Sound_to_Intensity.c */
