/* Sound_to_Intensity.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/07/03
 * pb 2002/07/16 GPL
 * pb 2003/05/20 default time step is four times oversampling
 */

#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 = 0.8 / minimumPitch;   /* Default: four times oversampling Hanning-wise. */
	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 */