/* Spectrum.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/05/22
 * pb 2001/11/30 Spectral moments
 * pb 2002/05/22 changed sign of imaginary part
 * pb 2002/07/16 GPL
 */

#include "Sound_and_Spectrum.h"

#include "oo_READ_ASCII.h"
#include "Spectrum_def.h"
#include "oo_READ_BINARY.h"
#include "Spectrum_def.h"

static void info (I) {
	iam (Spectrum);
	double minimum = 1e300, maximum = -1e300, sum = 0, mean, energy;
	long i;
	Melder_info ("Domain: %.8g Hz ... %.8g Hz", my xmin, my xmax);
	Melder_info ("Sampling:\n   number of frequencies: %ld\n   bin width: %.8g Hz\n   first bin: %.8g Hz",
		my nx, my dx, my x1);
	for (i = 1; i <= my nx; i ++) {
		double re = my z [1] [i], im = my z [2] [i];
		double squaredNorm = re * re + im * im;
		if (i != 1 && i != my nx) squaredNorm += squaredNorm;
		if (squaredNorm < minimum) minimum = squaredNorm;
		if (squaredNorm > maximum) maximum = squaredNorm;
		sum += squaredNorm;
	}
	mean = sum / my nx;
	Melder_info ("Norm:\n   minimum: %.8g\n   maximum: %.8g\n   rms: %.8g",
		sqrt (minimum), sqrt (maximum), sqrt (mean));
	energy = sum * my dx;
	Melder_info ("Total energy: %.8g Pa2 sec", energy);
}

class_methods (Spectrum, Matrix)
	us -> version = 1;   /* Changed sign of imaginary part. */
	class_method (info)
	class_method_local (Spectrum, readAscii)
	class_method_local (Spectrum, readBinary)
class_methods_end

Spectrum Spectrum_create (double fmax, long nf) {
	Spectrum me = new (Spectrum);
	if (! me || ! Matrix_init (me, 0.0, fmax, nf, fmax / (nf - 1), 0.0, 1, 2, 2, 1, 1))
		forget (me);
	return me;
}

int Spectrum_getPowerDensityRange (Spectrum me, double *minimum, double *maximum) {
	long ifreq;
	*minimum = 1e300, *maximum = 0;
	for (ifreq = 1; ifreq <= my nx; ifreq ++) {
		double oneSidedPowerSpectralDensity =   /* Pa2 Hz-2 s-1 */
			2 * (my z [1] [ifreq] * my z [1] [ifreq] + my z [2] [ifreq] * my z [2] [ifreq]) * my dx;
		if (oneSidedPowerSpectralDensity < *minimum) *minimum = oneSidedPowerSpectralDensity;
		if (oneSidedPowerSpectralDensity > *maximum) *maximum = oneSidedPowerSpectralDensity;
	}
	if (*maximum == 0.0) return 0;
	*minimum = 10 * log10 (*minimum / 4.0e-10);
	*maximum = 10 * log10 (*maximum / 4.0e-10);
	return 1;
}

static double Spectrum_sampleToPower (Spectrum me, long ifreq) {
	double re = my z [1] [ifreq], im = my z [2] [ifreq];   /* Pa / Hz */
	double oneSidedPowerSpectralDensity =   /* Pa^2 Hz-2 s-1 */
		2 * (re * re + im * im)
		* my dx;   /* divide by approximate duration */
	if (oneSidedPowerSpectralDensity == 0.0) return -1e30;
	return 10 * log10 (oneSidedPowerSpectralDensity /* Pa2/Hz */ / 4.0e-10 /* Pa2 */);   /* "dB/Hz" */
}

void Spectrum_drawInside (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum) {
	float *yWC = NULL;
	long ifreq, ifmin, ifmax;
	int autoscaling = minimum >= maximum;

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

	yWC = NUMfvector (ifmin, ifmax); cherror

	/*
	 * First pass: compute power density.
	 */
	if (autoscaling) maximum = -1e30;
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		double y = Spectrum_sampleToPower (me, ifreq);
		if (autoscaling && y > maximum) maximum = y;
		yWC [ifreq] = y;
		/*yWC [ifreq] = Spectrum_sampleToPower (me, ifreq);
		if (autoscaling && yWC [ifreq] > maximum) maximum = yWC [ifreq];*/
	}
	if (autoscaling) minimum = maximum - 60;   /* Default dynamic range is 60 dB. */

	/*
	 * Second pass: clip.
	 */
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		if (yWC [ifreq] < minimum) yWC [ifreq] = minimum;
		else if (yWC [ifreq] > maximum) yWC [ifreq] = maximum;
	}

	Graphics_setWindow (g, fmin, fmax, minimum, maximum);
	Graphics_function (g, yWC, ifmin, ifmax, Matrix_columnToX (me, ifmin), Matrix_columnToX (me, ifmax));
end:
	NUMfvector_free (yWC, ifmin);
	Melder_clearError ();
}

void Spectrum_draw (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum, int garnish) {
	Graphics_setInner (g);
	Spectrum_drawInside (me, g, fmin, fmax, minimum, maximum);
	Graphics_unsetInner (g);
	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, 1, "Frequency (Hz)");
		Graphics_marksBottom (g, 2, TRUE, TRUE, FALSE);
		Graphics_textLeft (g, 1, "Sound pressure level (dB/Hz)");
		Graphics_marksLeftEvery (g, 1.0, 20.0, TRUE, TRUE, FALSE);
	}
}

void Spectrum_drawLogFreq (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum, int garnish) {
	float *xWC = NULL, *yWC = NULL;
	long ifreq, ifmin, ifmax;
	int autoscaling = minimum >= maximum;
	if (fmax <= fmin) { fmin = my xmin; fmax = my xmax; }
	if (! Matrix_getWindowSamplesX (me, fmin, fmax, & ifmin, & ifmax)) return;
if(ifmin==1)ifmin=2;  /* BUG */
	xWC = NUMfvector (ifmin, ifmax); cherror
	yWC = NUMfvector (ifmin, ifmax); cherror

	/*
	 * First pass: compute power density.
	 */
	if (autoscaling) maximum = -1e6;
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		xWC [ifreq] = log10 (my x1 + (ifreq - 1) * my dx);
		yWC [ifreq] = Spectrum_sampleToPower (me, ifreq);
		if (autoscaling && yWC [ifreq] > maximum) maximum = yWC [ifreq];
	}
	if (autoscaling) minimum = maximum - 60;   /* Default dynamic range is 60 dB. */

	/*
	 * Second pass: clip.
	 */
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		if (yWC [ifreq] < minimum) yWC [ifreq] = minimum;
		else if (yWC [ifreq] > maximum) yWC [ifreq] = maximum;
	}

	Graphics_setInner (g);
	Graphics_setWindow (g, log10 (fmin), log10 (fmax), minimum, maximum);
	Graphics_polyline (g, ifmax - ifmin + 1, & xWC [ifmin], & yWC [ifmin]);
	Graphics_unsetInner (g);
	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, 1, "Frequency (Hz)");
		Graphics_marksBottomLogarithmic (g, 3, TRUE, TRUE, FALSE);
		Graphics_textLeft (g, 1, "Sound pressure level (dB/Hz)");
		Graphics_marksLeftEvery (g, 1.0, 20.0, TRUE, TRUE, FALSE);
	}
end:
	NUMfvector_free (xWC, ifmin);
	NUMfvector_free (yWC, ifmin);
	Melder_clearError ();
}

Spectrum Matrix_to_Spectrum (Matrix me) {
	Spectrum thee = NULL;
	if (my ny != 2)
		return Melder_errorp ("(Matrix_to_Spectrum:) Matrix must have exactly 2 rows.");
	if (! (thee = Data_copy (me))) return NULL;
	Thing_overrideClass (thee, classSpectrum);
	return thee;
}

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

Spectrum Spectrum_cepstralSmoothing (Spectrum me, double bandWidth) {
	Spectrum dBspectrum = NULL, thee = NULL;
	Sound cepstrum = NULL;
	float *re, *im;
	double factor = - bandWidth * bandWidth;
	long i;

	/*
	 * dB-spectrum is log (power).
	 */
	if (! (dBspectrum = Data_copy (me))) goto cleanUp;
	re = dBspectrum -> z [1], im = dBspectrum -> z [2];
	for (i = 1; i <= dBspectrum -> nx; i ++)
		{ re [i] = log (re [i] * re [i] + im [i] * im [i] + 1e-300); im [i] = 0; }

	/*
	 * Cepstrum is Fourier transform of dB-spectrum.
	 */
	if (! (cepstrum = Spectrum_to_Sound_fft (dBspectrum))) goto cleanUp;

	/*
	 * Multiply cepstrum by a Gaussian.
	 */
	for (i = 1; i <= cepstrum -> nx; i ++) {
		double t = cepstrum -> x1 + (i - 1) * cepstrum -> dx;
		cepstrum -> z [1] [i] *= exp (factor * t * t);
	}

	/*
	 * Smoothed power spectrum is original power spectrum convolved with a Gaussian.
	 */
	if (! (thee = Sound_to_Spectrum_fft (cepstrum))) goto cleanUp;

	/*
	 * Convert power spectrum back into a "complex" spectrum without phase information.
	 */
	re = thy z [1], im = thy z [2];
	for (i = 1; i <= thy nx; i ++)
		{ re [i] = exp (0.5 * re [i]); im [i] = 0; }   /* I.e., sqrt (exp (re [i])). */

cleanUp:
	forget (dBspectrum);
	forget (cepstrum);
	if (Melder_hasError ()) { forget (thee); Melder_error ("(Spectrum_cepstralSmoothing:) Not performed."); }
	return thee;
}

void Spectrum_passHannBand (Spectrum me, double fmin, double fmax0, double smooth) {
	long i;
	double fmax = fmax0 == 0.0 ? my xmax : fmax0;
	double f1 = fmin - smooth, f2 = fmin + smooth, f3 = fmax - smooth, f4 = fmax + smooth;
	double halfpibysmooth = smooth != 0.0 ? NUMpi / (2 * smooth) : 0.0;
	float *re = my z [1], *im = my z [2];
	for (i = 1; i <= my nx; i ++) {
		double frequency = my x1 + (i - 1) * my dx;
		if (frequency < f1 || frequency > f4) re [i] = im [i] = 0.0;
		if (frequency < f2 && fmin > 0.0) {
			double factor = 0.5 - 0.5 * cos (halfpibysmooth * (frequency - f1));
			re [i] *= factor;
			im [i] *= factor;
		} else if (frequency > f3 && fmax < my xmax) {
			double factor = 0.5 + 0.5 * cos (halfpibysmooth * (frequency - f3));
			re [i] *= factor;
			im [i] *= factor;
		}
	}
}

void Spectrum_stopHannBand (Spectrum me, double fmin, double fmax0, double smooth) {
	long i;
	double fmax = fmax0 == 0.0 ? my xmax : fmax0;
	double f1 = fmin - smooth, f2 = fmin + smooth, f3 = fmax - smooth, f4 = fmax + smooth;
	double halfpibysmooth = smooth != 0.0 ? NUMpi / (2 * smooth) : 0.0;
	float *re = my z [1], *im = my z [2];
	for (i = 1; i <= my nx; i ++) {
		double frequency = my x1 + (i - 1) * my dx;
		if (frequency < f1 || frequency > f4) continue;
		if (frequency < f2 && fmin > 0.0) {
			double factor = 0.5 + 0.5 * cos (halfpibysmooth * (frequency - f1));
			re [i] *= factor;
			im [i] *= factor;
		} else if (frequency > f3 && fmax < my xmax) {
			double factor = 0.5 - 0.5 * cos (halfpibysmooth * (frequency - f3));
			re [i] *= factor;
			im [i] *= factor;
		} else re [i] = im [i] = 0.0;
	}
}

double Spectrum_getBandDensity (Spectrum me, double fmin, double fmax) {
	long ifmin, ifmax, i;
	double sum = 0.0;
	if (fmax <= fmin) { fmin = my xmin; fmax = my xmax; }
	if (! Matrix_getWindowSamplesX (me, fmin, fmax, & ifmin, & ifmax)) return NUMundefined;
	for (i = ifmin; i <= ifmax; i ++) {
		double re = my z [1] [i], im = my z [2] [i];
		double squaredNorm = re * re + im * im;
		if (i != 1 && i != my nx) squaredNorm += squaredNorm;
		sum += squaredNorm;
	}
	/* energy = sum * dx; more accurately: energy = sum * (fmax - fmin) / (ifmax - ifmin + 1); */
	/* density = energy / (fmax - fmin); */
	return sum / (ifmax - ifmin + 1);
}

double Spectrum_getBandEnergy (Spectrum me, double fmin, double fmax) {
	double power;
	if (fmax <= fmin) { fmin = my xmin; fmax = my xmax; }
	power = Spectrum_getBandDensity (me, fmin, fmax);
	if (power == NUMundefined) return NUMundefined;
	return power * (fmax - fmin);
}

double Spectrum_getBandDensityDifference (Spectrum me, double lowBandMin, double lowBandMax, double highBandMin, double highBandMax) {
	double lowBandDensity = Spectrum_getBandDensity (me, lowBandMin, lowBandMax);
	double highBandDensity = Spectrum_getBandDensity (me, highBandMin, highBandMax);
	if (lowBandDensity == NUMundefined || highBandDensity == NUMundefined) return NUMundefined;
	if (lowBandDensity == 0.0 || highBandDensity == 0.0) return NUMundefined;
	return 10 * log10 (highBandDensity / lowBandDensity);
}

double Spectrum_getBandEnergyDifference (Spectrum me, double lowBandMin, double lowBandMax, double highBandMin, double highBandMax) {
	double lowBandEnergy = Spectrum_getBandEnergy (me, lowBandMin, lowBandMax);
	double highBandEnergy = Spectrum_getBandEnergy (me, highBandMin, highBandMax);
	if (lowBandEnergy == NUMundefined || highBandEnergy == NUMundefined) return NUMundefined;
	if (lowBandEnergy == 0.0 || highBandEnergy == 0.0) return NUMundefined;
	return 10 * log10 (highBandEnergy / lowBandEnergy);
}

double Spectrum_getCentreOfGravity (Spectrum me, double power) {
	long i;
	double halfpower = 0.5 * power, sumenergy = 0.0, sumfenergy = 0.0;
	for (i = 1; i <= my nx; i ++) {
		double re = my z [1] [i], im = my z [2] [i], energy = re * re + im * im;
		double f = my x1 + (i - 1) * my dx;
		if (halfpower != 1.0) energy = pow (energy, halfpower);
		sumenergy += energy;
		sumfenergy += f * energy;
	}
	return sumenergy == 0.0 ? NUMundefined : sumfenergy / sumenergy;
}

double Spectrum_getCentralMoment (Spectrum me, double moment, double power) {
	long i;
	double halfpower = 0.5 * power, sumenergy = 0.0, sumfenergy = 0.0;
	double fmean = Spectrum_getCentreOfGravity (me, power);
	if (fmean == NUMundefined) return NUMundefined;
	for (i = 1; i <= my nx; i ++) {
		double re = my z [1] [i], im = my z [2] [i], energy = re * re + im * im;
		double f = my x1 + (i - 1) * my dx;
		if (halfpower != 1.0) energy = pow (energy, halfpower);
		sumenergy += energy;
		sumfenergy += pow (f - fmean, moment) * energy;
	}
	return sumfenergy / sumenergy;
}

double Spectrum_getStandardDeviation (Spectrum me, double power) {
	return sqrt (Spectrum_getCentralMoment (me, 2.0, power));
}

double Spectrum_getSkewness (Spectrum me, double power) {
	double m2 = Spectrum_getCentralMoment (me, 2.0, power);
	double m3 = Spectrum_getCentralMoment (me, 3.0, power);
	if (m2 == NUMundefined || m3 == NUMundefined || m2 == 0) return NUMundefined;
	return m3 / (m2 * sqrt (m2));
}

double Spectrum_getKurtosis (Spectrum me, double power) {
	double m2 = Spectrum_getCentralMoment (me, 2.0, power);
	double m4 = Spectrum_getCentralMoment (me, 4.0, power);
	if (m2 == NUMundefined || m4 == NUMundefined || m2 == 0) return NUMundefined;
	return m4 / (m2 * m2) - 3;
}

/* End of file Spectrum.c */