/* Ltas.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 2002/06/04
 * pb 2002/07/16 GPL
 * pb 2003/04/16 Ltases_merge
 * pb 2003/06/19 Ltas_computeTrendLine, Ltas_subtractTrendLine
 */

#include "Ltas.h"

class_methods (Ltas, Vector)
class_methods_end

Ltas Ltas_create (long nx, double dx) {
	Ltas me = new (Ltas);
	if (! me || ! Matrix_init (me, 0, nx * dx, nx, dx, 0.5 * dx, 1, 1, 1, 1, 1)) return NULL;
	return me;
}

void Ltas_draw (Ltas me, Graphics g, double fmin, double fmax, double minimum, double maximum, int garnish) {
	long ifreq, ifmin, ifmax;

	if (fmax <= fmin) { fmin = my xmin; fmax = my xmax; }
	if (! Matrix_getWindowSamplesX (me, fmin, fmax, & ifmin, & ifmax)) return;

	Graphics_setWindow (g, fmin, fmax, minimum, maximum);
	Graphics_setInner (g);
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		double yWC = my z [1] [ifreq];
		double f = Sampled_indexToX (me, ifreq);
		double left = f - 0.5 * my dx, right = f + 0.5 * my dx;
		if (yWC < minimum) yWC = minimum;
		if (left < fmin) left = fmin;
		if (right > fmax) right = fmax;
		Graphics_line (g, left, yWC, right, yWC);
		Graphics_line (g, left, yWC, left, minimum);
		Graphics_line (g, right, yWC, right, minimum);
	}
	Graphics_unsetInner (g);

	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, 1, "Frequency (Hz)");
		Graphics_marksBottom (g, 2, 1, 1, 0);
		Graphics_textLeft (g, 1, "Sound pressure level (dB/Hz)");
		Graphics_marksLeft (g, 2, 1, 1, 0);
	}
}

Matrix Ltas_to_Matrix (Ltas me) {
	Matrix thee = Data_copy (me);
	if (! thee) return NULL;
	Thing_overrideClass (thee, classMatrix);
	return thee;
}

Ltas Matrix_to_Ltas (Matrix me) {
	Ltas thee = Data_copy (me);
	if (! thee) return NULL;
	Thing_overrideClass (thee, classLtas);
	return thee;
}

Ltas Ltases_merge (Collection ltases) {
	Ltas me, thee = NULL;
	long iband, ispec;
	if (ltases -> size < 1) { Melder_error ("Cannot merge zero Ltas objects."); goto end; }
	me = ltases -> item [1];
	thee = Data_copy (me); cherror
	/*
	 * Convert to energy.
	 */
	for (iband = 1; iband <= thy nx; iband ++) {
		thy z [1] [iband] = pow (10.0, thy z [1] [iband] / 10.0);
	}
	for (ispec = 2; ispec <= ltases -> size; ispec ++) {
		Ltas him = ltases -> item [ispec];
		if (his xmin != thy xmin || his xmax != thy xmax) { Melder_error ("Frequency domains do not match."); goto end; }
		if (his dx != thy dx) { Melder_error ("Bandwidths do not match."); goto end; }
		if (his nx != thy nx || his x1 != thy x1) { Melder_error ("Frequency bands do not match."); goto end; }
		/*
		 * Add band energies.
		 */
		for (iband = 1; iband <= thy nx; iband ++) {
			thy z [1] [iband] += pow (10.0, his z [1] [iband] / 10.0);
		}
	}
	/*
	 * Convert back to dB.
	 */
	for (iband = 1; iband <= thy nx; iband ++) {
		thy z [1] [iband] = (10.0 * NUMlog10e) * log (thy z [1] [iband]);
	}
end:
	iferror forget (thee);
	return thee;
}

Ltas Ltas_computeTrendLine (Ltas me, double fmin, double fmax) {
	Ltas thee = NULL;
	long imin, imax, n, i;
	double sum = 0.0, amean, fmean, numerator = 0.0, denominator = 0.0, slope;
	/*
	 * Find the first and last bin.
	 */
	if ((n = Sampled_getWindowSamples (me, fmin, fmax, & imin, & imax)) < 2) {
		return Melder_errorp ("(Ltas_computeTrendLine:) Number of bins too low (%ld). Should be at least 2.", n);
	}
	thee = Data_copy (me);
	/*
	 * Compute average amplitude and frequency.
	 */
	for (i = imin; i <= imax; i ++) {
		sum += thy z [1] [i];
	}
	amean = sum / n;
	fmean = thy x1 + (0.5 * (imin + imax) - 1) * thy dx;
	/*
	 * Compute slope.
	 */
	for (i = imin; i <= imax; i ++) {
		double da = thy z [1] [i] - amean, df = thy x1 + (i - 1) * thy dx - fmean;
		numerator += da * df;
		denominator += df * df;
	}
	slope = numerator / denominator;
	/*
	 * Modify bins.
	 */
	for (i = 1; i <= thy nx; i ++) {
		double df = thy x1 + (i - 1) * thy dx - fmean;
		thy z [1] [i] = amean + slope * df;
	}
	return thee;
}

Ltas Ltas_subtractTrendLine (Ltas me, double fmin, double fmax) {
	Ltas thee = NULL;
	long imin, imax, n, i;
	double sum = 0.0, amean, fmean, numerator = 0.0, denominator = 0.0, slope;
	/*
	 * Find the first and last bin.
	 */
	if ((n = Sampled_getWindowSamples (me, fmin, fmax, & imin, & imax)) < 2) {
		return Melder_errorp ("(Ltas_subtractTrendLine:) Number of bins too low (%ld). Should be at least 2.", n);
	}
	thee = Data_copy (me);
	/*
	 * Compute average amplitude and frequency.
	 */
	for (i = imin; i <= imax; i ++) {
		sum += thy z [1] [i];
	}
	amean = sum / n;
	fmean = thy x1 + (0.5 * (imin + imax) - 1) * thy dx;
	/*
	 * Compute slope.
	 */
	for (i = imin; i <= imax; i ++) {
		double da = thy z [1] [i] - amean, df = thy x1 + (i - 1) * thy dx - fmean;
		numerator += da * df;
		denominator += df * df;
	}
	slope = numerator / denominator;
	/*
	 * Modify bins.
	 */
	for (i = 1; i < imin; i ++) {
		thy z [1] [i] = 0.0;
	}
	for (i = imin; i <= imax; i ++) {
		double df = thy x1 + (i - 1) * thy dx - fmean;
		thy z [1] [i] -= amean + slope * df;
	}
	for (i = imax + 1; i <= thy nx; i ++) {
		thy z [1] [i] = 0.0;
	}
	return thee;
}

/* End of file Ltas.c */
