/* Spectrum_to_Ltas.c
 *
 * Copyright (C) 1992-2002 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 1998/01/05
 * pb 2002/07/16 GPL
 */

#include "Spectrum_to_Ltas.h"

Ltas Spectrum_to_Ltas (Spectrum me, double bandWidth) {
	Ltas thee = NULL;
	long numberOfBands = ceil ((my xmax - my xmin) / bandWidth);
	long i, j;
	if (bandWidth <= my dx) {
		Melder_error ("Bandwidth must be greater than %f.", my dx);
		goto error;
	}
	if (! (thee = new (Ltas)) ||
		! Matrix_init (thee, my xmin, my xmax, numberOfBands,
			bandWidth, my xmin + 0.5 * bandWidth, 1, 1, 1, 1, 1))
		goto error;
	for (i = 1; i <= numberOfBands; i ++) {
		double fmin = my xmin + (i - 1) * bandWidth;
		double fmax = my xmin + i * bandWidth;
		long jmin = Sampled_xToHighIndex (me, fmin);
		long jmax = Sampled_xToLowIndex (me, fmax);
		if (jmin < 1) jmin = 1; else if (jmin > my nx) jmin = my nx;
		if (jmax < 1) jmax = 1; else if (jmax > my nx) jmax = my nx;
		for (j = jmin; j <= jmax; j ++)
			thy z [1] [i] += my z [1] [j] * my z [1] [j] + my z [2] [j] * my z [2] [j];
		if (thy z [1] [i] != 0.0) {
			double oneSidedPowerSpectralDensity = 2 * thy z [1] [i] / (jmax - jmin + 1) * my dx;
			thy z [1] [i] = 10 * log10 (oneSidedPowerSpectralDensity / 4.0e-10);
		}
	}
	return thee;
error:
	forget (thee);
	return Melder_errorp ("Spectrum_to_Ltas: not performed.");
}

Ltas Spectrum_to_Ltas_1to1 (Spectrum me) {
	long i;
	Ltas thee = new (Ltas); cherror
	Matrix_init (thee, my xmin, my xmax, my nx, my dx, my x1, 1, 1, 1, 1, 1); cherror
	for (i = 1; i <= my nx; i ++) {
		double value = my z [1] [i] * my z [1] [i] + my z [2] [i] * my z [2] [i];
		if (value != 0.0) {
			double oneSidedPowerSpectralDensity = 2 * value * my dx;
			value = 10 * log10 (oneSidedPowerSpectralDensity / 4.0e-10);
		}
		thy z [1] [i] = value;
	}
end:
	iferror forget (thee);
	return thee;
}

/* End of file Spectrum_to_Ltas.c */