/* Pitch.cpp
 *
 * Copyright (C) 1992-2009,2011,2012,2014-2018 Paul Boersma
 *
 * This code 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 code 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 work. If not, see <http://www.gnu.org/licenses/>.
 */

#include "Pitch.h"
#include <ctype.h>
#include "Sound_and_Spectrum.h"
#include "Matrix_and_Pitch.h"

#include "oo_DESTROY.h"
#include "Pitch_def.h"
#include "oo_COPY.h"
#include "Pitch_def.h"
#include "oo_EQUAL.h"
#include "Pitch_def.h"
#include "oo_CAN_WRITE_AS_ENCODING.h"
#include "Pitch_def.h"
#include "oo_WRITE_TEXT.h"
#include "Pitch_def.h"
#include "oo_READ_TEXT.h"
#include "Pitch_def.h"
#include "oo_WRITE_BINARY.h"
#include "Pitch_def.h"
#include "oo_READ_BINARY.h"
#include "Pitch_def.h"
#include "oo_DESCRIPTION.h"
#include "Pitch_def.h"

#include "enums_getText.h"
#include "Pitch_enums.h"
#include "enums_getValue.h"
#include "Pitch_enums.h"

Thing_implement (Pitch, Sampled, 1);

int structPitch :: v_getMinimumUnit (integer level) {
	return level == Pitch_LEVEL_FREQUENCY ? (int) kPitch_unit::MIN : Pitch_STRENGTH_UNIT_min;
}

int structPitch :: v_getMaximumUnit (integer level) {
	return level == Pitch_LEVEL_FREQUENCY ? (int) kPitch_unit::MAX : Pitch_STRENGTH_UNIT_max;
}

const char32 * structPitch :: v_getUnitText (integer level, int unit, uint32 flags) {
	if (level == Pitch_LEVEL_FREQUENCY) {
		return
			unit == (int) kPitch_unit::HERTZ ?
				flags & Function_UNIT_TEXT_MENU ? U"Hertz" : U"Hz" :
			unit == (int) kPitch_unit::HERTZ_LOGARITHMIC ?
				flags & Function_UNIT_TEXT_MENU ? U"Hertz (logarithmic)" : (flags & Function_UNIT_TEXT_SHORT) && (flags & Function_UNIT_TEXT_GRAPHICAL) ? U"%%Hz%" : U"Hz" :
			unit == (int) kPitch_unit::MEL ? U"mel" :
			unit == (int) kPitch_unit::LOG_HERTZ ?
				flags & Function_UNIT_TEXT_MENU ? U"logHertz" : U"logHz" :
			unit == (int) kPitch_unit::SEMITONES_1 ?
				flags & Function_UNIT_TEXT_SHORT ? U"st__1_" : flags & Function_UNIT_TEXT_GRAPHICAL ? U"semitones %%re% 1 Hz" : U"semitones re 1 Hz" :
			unit == (int) kPitch_unit::SEMITONES_100 ?
				flags & Function_UNIT_TEXT_SHORT ? U"st__100_" : flags & Function_UNIT_TEXT_GRAPHICAL ? U"semitones %%re% 100 Hz" : U"semitones re 100 Hz" :
			unit == (int) kPitch_unit::SEMITONES_200 ?
				flags & Function_UNIT_TEXT_SHORT ? U"st__200_" : flags & Function_UNIT_TEXT_GRAPHICAL ? U"semitones %%re% 200 Hz" : U"semitones re 200 Hz" :
			unit == (int) kPitch_unit::SEMITONES_440 ?
				flags & Function_UNIT_TEXT_SHORT ? U"st__a_" : flags & Function_UNIT_TEXT_GRAPHICAL ? U"semitones %%re% 440 Hz" : U"semitones re 440 Hz" :
			unit == (int) kPitch_unit::ERB ?
				flags & Function_UNIT_TEXT_SHORT ? U"erb" : U"ERB" :
			U"";
	} else if (level == Pitch_LEVEL_STRENGTH) {
		return
			unit == Pitch_STRENGTH_UNIT_AUTOCORRELATION ? U"" :
			unit == Pitch_STRENGTH_UNIT_NOISE_HARMONICS_RATIO ? U"" :
			unit == Pitch_STRENGTH_UNIT_HARMONICS_NOISE_DB ? U"dB" :
			U"";
	}
	return U"unknown";
}

bool structPitch :: v_isUnitLogarithmic (integer level, int unit) {
	return level == Pitch_LEVEL_FREQUENCY && unit == (int) kPitch_unit::HERTZ_LOGARITHMIC;
}

double structPitch :: v_convertStandardToSpecialUnit (double value, integer level, int unit) {
	if (level == Pitch_LEVEL_FREQUENCY) {
		return
			unit == (int) kPitch_unit::HERTZ ? value :
			unit == (int) kPitch_unit::HERTZ_LOGARITHMIC ? value <= 0.0 ? undefined : log10 (value) :
			unit == (int) kPitch_unit::MEL ? NUMhertzToMel (value) :
			unit == (int) kPitch_unit::LOG_HERTZ ? value <= 0.0 ? undefined : log10 (value) :
			unit == (int) kPitch_unit::SEMITONES_1 ? value <= 0.0 ? undefined : (12.0 / NUMln2) * log (value / 1.0) :
			unit == (int) kPitch_unit::SEMITONES_100 ? value <= 0.0 ? undefined : (12.0 / NUMln2) * log (value / 100.0) :
			unit == (int) kPitch_unit::SEMITONES_200 ? value <= 0.0 ? undefined : (12.0 / NUMln2) * log (value / 200.0) :
			unit == (int) kPitch_unit::SEMITONES_440 ? value <= 0.0 ? undefined : (12.0 / NUMln2) * log (value / 440.0) :
			unit == (int) kPitch_unit::ERB ? NUMhertzToErb (value) :
			undefined;
	} else {
		return
			unit == Pitch_STRENGTH_UNIT_AUTOCORRELATION ? value :
			unit == Pitch_STRENGTH_UNIT_NOISE_HARMONICS_RATIO ?
				value <= 1e-15 ? 1e15 : value > 1.0 - 1e-15 ? 1e-15 : (1.0 - value) / value :   // before losing precision
			unit == Pitch_STRENGTH_UNIT_HARMONICS_NOISE_DB ?
				value <= 1e-15 ? -150.0 : value > 1.0 - 1e-15 ? 150.0 : 10.0 * log10 (value / (1.0 - value)) :   // before losing precision
			undefined;
	}
}

double structPitch :: v_convertSpecialToStandardUnit (double value, integer level, int unit) {
	if (level == Pitch_LEVEL_FREQUENCY) {
		return
			unit == (int) kPitch_unit::HERTZ ? value :
			unit == (int) kPitch_unit::HERTZ_LOGARITHMIC ? pow (10.0, value) :
			unit == (int) kPitch_unit::MEL ? NUMmelToHertz (value) :
			unit == (int) kPitch_unit::LOG_HERTZ ? pow (10.0, value) :
			unit == (int) kPitch_unit::SEMITONES_1 ? 1.0 * exp (value * (NUMln2 / 12.0)):
			unit == (int) kPitch_unit::SEMITONES_100 ? 100.0 * exp (value * (NUMln2 / 12.0)):
			unit == (int) kPitch_unit::SEMITONES_200 ? 200.0 * exp (value * (NUMln2 / 12.0)):
			unit == (int) kPitch_unit::SEMITONES_440 ? 440.0 * exp (value * (NUMln2 / 12.0)):
			unit == (int) kPitch_unit::ERB ? NUMerbToHertz (value) :
			undefined;
	} else {
		return undefined;
	}
}

#define doesUnitAllowNegativeValues(unit)  \
	( (unit) == kPitch_unit::HERTZ_LOGARITHMIC || (unit) == kPitch_unit::LOG_HERTZ ||  \
	  (unit) == kPitch_unit::SEMITONES_1 || (unit) == kPitch_unit::SEMITONES_100 ||  \
	  (unit) == kPitch_unit::SEMITONES_200 || (unit) == kPitch_unit::SEMITONES_440 )

double structPitch :: v_getValueAtSample (integer iframe, integer ilevel, int unit) {
	double f = our frame [iframe]. candidate [1]. frequency;
	if (! Pitch_util_frequencyIsVoiced (f, our ceiling)) return undefined;
	return v_convertStandardToSpecialUnit (ilevel == Pitch_LEVEL_FREQUENCY ? f : our frame [iframe]. candidate [1]. strength, ilevel, unit);
}

bool Pitch_isVoiced_i (Pitch me, integer iframe) {
	return isdefined (Sampled_getValueAtSample (me, iframe, Pitch_LEVEL_FREQUENCY, (int) kPitch_unit::HERTZ));
}

bool Pitch_isVoiced_t (Pitch me, double time) {
	return isdefined (Sampled_getValueAtX (me, time, Pitch_LEVEL_FREQUENCY, (int) kPitch_unit::HERTZ, false));
}

double Pitch_getValueAtTime (Pitch me, double time, kPitch_unit unit, bool interpolate) {
	return Sampled_getValueAtX (me, time, Pitch_LEVEL_FREQUENCY, (int) unit, interpolate);
}

double Pitch_getStrengthAtTime (Pitch me, double time, kPitch_unit unit, bool interpolate) {
	return Sampled_getValueAtX (me, time, Pitch_LEVEL_STRENGTH, (int) unit, interpolate);
}

integer Pitch_countVoicedFrames (Pitch me) {
	return Sampled_countDefinedSamples (me, Pitch_LEVEL_FREQUENCY, (int) kPitch_unit::HERTZ);
}

double Pitch_getMean (Pitch me, double tmin, double tmax, kPitch_unit unit) {
	return Sampled_getMean (me, tmin, tmax, Pitch_LEVEL_FREQUENCY, (int) unit, true);
}

double Pitch_getMeanStrength (Pitch me, double tmin, double tmax, int strengthUnit) {
	return Sampled_getMean (me, tmin, tmax, Pitch_LEVEL_STRENGTH, strengthUnit, true);
}

double Pitch_getQuantile (Pitch me, double tmin, double tmax, double quantile, kPitch_unit unit) {
	double value = Sampled_getQuantile (me, tmin, tmax, quantile, Pitch_LEVEL_FREQUENCY, (int) unit);
	if (value <= 0.0 && ! doesUnitAllowNegativeValues (unit)) {
		value = undefined;
	}
	return value;
}

double Pitch_getStandardDeviation (Pitch me, double tmin, double tmax, kPitch_unit unit) {
	return Sampled_getStandardDeviation (me, tmin, tmax, Pitch_LEVEL_FREQUENCY, (int) unit, true);
}

void Pitch_getMaximumAndTime (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate,
	double *return_maximum, double *return_timeOfMaximum)
{
	Sampled_getMaximumAndX (me, tmin, tmax, Pitch_LEVEL_FREQUENCY, (int) unit, interpolate, return_maximum, return_timeOfMaximum);
	if (! doesUnitAllowNegativeValues (unit) && return_maximum && *return_maximum <= 0.0)
	{
		*return_maximum = undefined;   // unlikely
	}
}

double Pitch_getMaximum (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate) {
	double maximum;
	Pitch_getMaximumAndTime (me, tmin, tmax, unit, interpolate, & maximum, nullptr);
	return maximum;
}

double Pitch_getTimeOfMaximum (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate) {
	double time;
	Pitch_getMaximumAndTime (me, tmin, tmax, unit, interpolate, nullptr, & time);
	return time;
}

void Pitch_getMinimumAndTime (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate,
	double *return_minimum, double *return_timeOfMinimum)
{
	Sampled_getMinimumAndX (me, tmin, tmax, Pitch_LEVEL_FREQUENCY, (int) unit, interpolate, return_minimum, return_timeOfMinimum);
	if (! doesUnitAllowNegativeValues (unit) && return_minimum && *return_minimum <= 0.0)
	{
		*return_minimum = undefined;   // not so unlikely
	}
}

double Pitch_getMinimum (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate) {
	double minimum;
	Pitch_getMinimumAndTime (me, tmin, tmax, unit, interpolate, & minimum, nullptr);
	return minimum;
}

double Pitch_getTimeOfMinimum (Pitch me, double tmin, double tmax, kPitch_unit unit, bool interpolate) {
	double time;
	Pitch_getMinimumAndTime (me, tmin, tmax, unit, interpolate, nullptr, & time);
	return time;
}

static integer Pitch_getMeanAbsoluteSlope (Pitch me,
	double *out_hertz, double *out_mel, double *out_semitones, double *out_erb, double *out_withoutOctaveJumps)
{
	integer firstVoicedFrame = 0, lastVoicedFrame = 0, numberOfVoicedFrames = 0;
	autoNUMvector <double> frequencies (1, my nx);
	for (integer i = 1; i <= my nx; i ++) {
		double frequency = my frame [i]. candidate [1]. frequency;
		if (Pitch_util_frequencyIsVoiced (frequency, my ceiling)) {
			frequencies [i] = frequency;
			numberOfVoicedFrames ++;
		} else {
			frequencies [i] = 0.0;
		}
	}
	for (integer i = 1; i <= my nx; i ++)   // look for first voiced frame
		if (frequencies [i] != 0.0) { firstVoicedFrame = i; break; }
	for (integer i = my nx; i >= 1; i --)   // look for last voiced frame
		if (frequencies [i] != 0.0) { lastVoicedFrame = i; break; }
	if (numberOfVoicedFrames > 1) {
		integer ilast = firstVoicedFrame;
		double span = (lastVoicedFrame - firstVoicedFrame) * my dx, flast = frequencies [ilast];
		double slopeHz = 0.0, slopeMel = 0.0, slopeSemitones = 0.0, slopeErb = 0.0, slopeRobust = 0.0;
		for (integer i = firstVoicedFrame + 1; i <= lastVoicedFrame; i ++) if (frequencies [i] != 0.0) {
			double localStepSemitones = fabs (NUMhertzToSemitones (frequencies [i]) - NUMhertzToSemitones (flast));
			slopeHz += fabs (frequencies [i] - flast);
			slopeMel += fabs (NUMhertzToMel (frequencies [i]) - NUMhertzToMel (flast));
			slopeSemitones += localStepSemitones;
			slopeErb += fabs (NUMhertzToErb (frequencies [i]) - NUMhertzToErb (flast));
			while (localStepSemitones >= 12.0) localStepSemitones -= 12.0;   // kill octave jumps
			if (localStepSemitones > 6.0) localStepSemitones = 12.0 - localStepSemitones;
			slopeRobust += localStepSemitones;
			ilast = i;
			flast = frequencies [ilast];
		}
		if (out_hertz) *out_hertz = slopeHz / span;
		if (out_mel) *out_mel = slopeMel / span;
		if (out_semitones) *out_semitones = slopeSemitones / span;
		if (out_erb) *out_erb = slopeErb / span;
		if (out_withoutOctaveJumps) *out_withoutOctaveJumps = slopeRobust / span;
	} else {
		if (out_hertz) *out_hertz = undefined;
		if (out_mel) *out_mel = undefined;
		if (out_semitones) *out_semitones = undefined;
		if (out_erb) *out_erb = undefined;
		if (out_withoutOctaveJumps) *out_withoutOctaveJumps = undefined;
	}
	return numberOfVoicedFrames;
}

integer Pitch_getMeanAbsSlope_hertz (Pitch me, double *slope) {
	return Pitch_getMeanAbsoluteSlope (me, slope, nullptr, nullptr, nullptr, nullptr);
}

integer Pitch_getMeanAbsSlope_mel (Pitch me, double *slope) {
	return Pitch_getMeanAbsoluteSlope (me, nullptr, slope, nullptr, nullptr, nullptr);
}

integer Pitch_getMeanAbsSlope_semitones (Pitch me, double *slope) {
	return Pitch_getMeanAbsoluteSlope (me, nullptr, nullptr, slope, nullptr, nullptr);
}

integer Pitch_getMeanAbsSlope_erb (Pitch me, double *slope) {
	return Pitch_getMeanAbsoluteSlope (me, nullptr, nullptr, nullptr, slope, nullptr);
}

integer Pitch_getMeanAbsSlope_noOctave (Pitch me, double *slope) {
	return Pitch_getMeanAbsoluteSlope (me, nullptr, nullptr, nullptr, nullptr, slope);
}

void structPitch :: v_info () {
	integer nVoiced;
	autoNUMvector <double> frequencies (Sampled_getSortedValues (this, Pitch_LEVEL_FREQUENCY, (int) kPitch_unit::HERTZ, & nVoiced), 1);
	structDaata :: v_info ();
	MelderInfo_writeLine (U"Time domain:");
	MelderInfo_writeLine (U"   Start time: ", xmin, U" seconds");
	MelderInfo_writeLine (U"   End time: ", xmax, U" seconds");
	MelderInfo_writeLine (U"   Total duration: ", xmax - xmin, U" seconds");
	MelderInfo_writeLine (U"Time sampling:");
	MelderInfo_writeLine (U"   Number of frames: ", nx, U" (", nVoiced, U" voiced)");
	MelderInfo_writeLine (U"   Time step: ", dx, U" seconds");
	MelderInfo_writeLine (U"   First frame centred at: ", x1, U" seconds");
	MelderInfo_writeLine (U"Ceiling at: ", ceiling, U" Hz");

	if (nVoiced >= 1) {   // quantiles
		double quantile10, quantile16, quantile50, quantile84, quantile90;
		quantile10 = NUMquantile (nVoiced, frequencies.peek(), 0.10);
		quantile16 = NUMquantile (nVoiced, frequencies.peek(), 0.16);
		quantile50 = NUMquantile (nVoiced, frequencies.peek(), 0.50);   // median
		quantile84 = NUMquantile (nVoiced, frequencies.peek(), 0.84);
		quantile90 = NUMquantile (nVoiced, frequencies.peek(), 0.90);
		MelderInfo_writeLine (U"\nEstimated quantiles:");
		MelderInfo_write (U"   10% = ", Melder_single (quantile10), U" Hz = ", Melder_single (NUMhertzToMel (quantile10)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (quantile10)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (quantile10)), U" ERB");
		MelderInfo_write (U"   16% = ", Melder_single (quantile16), U" Hz = ", Melder_single (NUMhertzToMel (quantile16)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (quantile16)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (quantile16)), U" ERB");
		MelderInfo_write (U"   50% = ", Melder_single (quantile50), U" Hz = ", Melder_single (NUMhertzToMel (quantile50)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (quantile50)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (quantile50)), U" ERB");
		MelderInfo_write (U"   84% = ", Melder_single (quantile84), U" Hz = ", Melder_single (NUMhertzToMel (quantile84)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (quantile84)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (quantile84)), U" ERB");
		MelderInfo_write (U"   90% = ", Melder_single (quantile90), U" Hz = ", Melder_single (NUMhertzToMel (quantile90)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (quantile90)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (quantile90)), U" ERB");
		if (nVoiced > 1) {
			double corr = sqrt (nVoiced / (nVoiced - 1.0));
			MelderInfo_writeLine (U"\nEstimated spreading:");
			MelderInfo_write (U"   84%-median = ", Melder_half ((quantile84 - quantile50) * corr), U" Hz = ", Melder_half ((NUMhertzToMel (quantile84) - NUMhertzToMel (quantile50)) * corr), U" Mel = ");
			MelderInfo_writeLine (Melder_half ((NUMhertzToSemitones (quantile84) - NUMhertzToSemitones (quantile50)) * corr), U" semitones = ", Melder_half ((NUMhertzToErb (quantile84) - NUMhertzToErb (quantile50)) * corr), U" ERB");
			MelderInfo_write (U"   median-16% = ", Melder_half ((quantile50 - quantile16) * corr), U" Hz = ", Melder_half ((NUMhertzToMel (quantile50) - NUMhertzToMel (quantile16)) * corr), U" Mel = ");
			MelderInfo_writeLine (Melder_half ((NUMhertzToSemitones (quantile50) - NUMhertzToSemitones (quantile16)) * corr), U" semitones = ", Melder_half ((NUMhertzToErb (quantile50) - NUMhertzToErb (quantile16)) * corr), U" ERB");
			MelderInfo_write (U"   90%-10% = ", Melder_half ((quantile90 - quantile10) * corr), U" Hz = ", Melder_half ((NUMhertzToMel (quantile90) - NUMhertzToMel (quantile10)) * corr), U" Mel = ");
			MelderInfo_writeLine (Melder_half ((NUMhertzToSemitones (quantile90) - NUMhertzToSemitones (quantile10)) * corr), U" semitones = ", Melder_half ((NUMhertzToErb (quantile90) - NUMhertzToErb (quantile10)) * corr), U" ERB");
		}
	}
	if (nVoiced >= 1) {   // extrema, range, mean and standard deviation
		double minimum = Pitch_getMinimum (this, xmin, xmax, kPitch_unit::HERTZ, false);
		double maximum = Pitch_getMaximum (this, xmin, xmax, kPitch_unit::HERTZ, false);
		double meanHertz, meanMel, meanSemitones, meanErb;
		MelderInfo_write (U"\nMinimum ", Melder_single (minimum), U" Hz = ", Melder_single (NUMhertzToMel (minimum)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (minimum)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (minimum)), U" ERB");
		MelderInfo_write (U"Maximum ", Melder_single (maximum), U" Hz = ", Melder_single (NUMhertzToMel (maximum)), U" Mel = ");
		MelderInfo_writeLine (Melder_single (NUMhertzToSemitones (maximum)), U" semitones above 100 Hz = ", Melder_single (NUMhertzToErb (maximum)), U" ERB");
		MelderInfo_write (U"Range ", Melder_half (maximum - minimum), U" Hz = ", Melder_single (NUMhertzToMel (maximum) - NUMhertzToMel (minimum)), U" Mel = ");
		MelderInfo_writeLine (Melder_half (NUMhertzToSemitones (maximum) - NUMhertzToSemitones (minimum)), U" semitones = ", Melder_half (NUMhertzToErb (maximum) - NUMhertzToErb (minimum)), U" ERB");
		meanHertz = Pitch_getMean (this, 0, 0, kPitch_unit::HERTZ);
		meanMel = Pitch_getMean (this, 0, 0, kPitch_unit::MEL);
		meanSemitones = Pitch_getMean (this, 0, 0, kPitch_unit::SEMITONES_100);
		meanErb = Pitch_getMean (this, 0, 0, kPitch_unit::ERB);
		MelderInfo_write (U"Average: ", Melder_single (meanHertz), U" Hz = ", Melder_single (meanMel), U" Mel = ");
		MelderInfo_writeLine (Melder_single (meanSemitones), U" semitones above 100 Hz = ", Melder_single (meanErb), U" ERB");
		if (nVoiced >= 2) {
			double stdevHertz = Pitch_getStandardDeviation (this, 0, 0, kPitch_unit::HERTZ);
			double stdevMel = Pitch_getStandardDeviation (this, 0, 0, kPitch_unit::MEL);
			double stdevSemitones = Pitch_getStandardDeviation (this, 0, 0, kPitch_unit::SEMITONES_100);
			double stdevErb = Pitch_getStandardDeviation (this, 0, 0, kPitch_unit::ERB);
			MelderInfo_write (U"Standard deviation: ", Melder_half (stdevHertz), U" Hz = ", Melder_half (stdevMel), U" Mel = ");
			MelderInfo_writeLine (Melder_half (stdevSemitones), U" semitones = ", Melder_half (stdevErb), U" ERB");
		}
	}
	if (nVoiced > 1) {   // variability: mean absolute slope
		double slopeHertz, slopeMel, slopeSemitones, slopeErb, slopeWithoutOctaveJumps;
		Pitch_getMeanAbsoluteSlope (this, & slopeHertz, & slopeMel, & slopeSemitones, & slopeErb, & slopeWithoutOctaveJumps);
		MelderInfo_write (U"\nMean absolute slope: ", Melder_half (slopeHertz), U" Hz/s = ", Melder_half (slopeMel), U" Mel/s = ");
		MelderInfo_writeLine (Melder_half (slopeSemitones), U" semitones/s = ", Melder_half (slopeErb), U" ERB/s");
		MelderInfo_writeLine (U"Mean absolute slope without octave jumps: ", Melder_half (slopeWithoutOctaveJumps), U" semitones/s");
	}
}

void Pitch_Frame_init (Pitch_Frame me, int nCandidates) {
	/*
	 * Create without change.
	 */
	autoNUMvector <structPitch_Candidate> candidate (1, nCandidates);
	/*
	 * Change without error.
	 */
	NUMvector_free (my candidate, 1);
	my candidate = candidate.transfer();
	my nCandidates = nCandidates;
}

autoPitch Pitch_create (double tmin, double tmax, integer nt, double dt, double t1,
	double ceiling, int maxnCandidates)
{
	try {
		autoPitch me = Thing_new (Pitch);
		Sampled_init (me.get(), tmin, tmax, nt, dt, t1);
		my ceiling = ceiling;
		my maxnCandidates = maxnCandidates;
		my frame = NUMvector <structPitch_Frame> (1, nt);

		/* Put one candidate in every frame (unvoiced, silent). */
		for (integer it = 1; it <= nt; it ++) {
			Pitch_Frame_init (& my frame [it], 1);
		}

		return me;
	} catch (MelderError) {
		Melder_throw (U"Pitch not created.");
	}
}

void Pitch_setCeiling (Pitch me, double ceiling) {
	my ceiling = ceiling;
}

int Pitch_getMaxnCandidates (Pitch me) {
	int result = 0;
	for (integer i = 1; i <= my nx; i ++) {
		int nCandidates = my frame [i]. nCandidates;
		if (nCandidates > result) result = nCandidates;
	}
	return result;
}

void Pitch_pathFinder (Pitch me, double silenceThreshold, double voicingThreshold,
	double octaveCost, double octaveJumpCost, double voicedUnvoicedCost,
	double ceiling, int pullFormants)
{
	if (Melder_debug == 33)
		Melder_casual (U"Pitch path finder:"
			U"\nSilence threshold = ", silenceThreshold,
			U"\nVoicing threshold = ", voicingThreshold,
			U"\nOctave cost = ", octaveCost,
			U"\nOctave jump cost = ", octaveJumpCost,
			U"\nVoiced/unvoiced cost = ", voicedUnvoicedCost,
			U"\nCeiling = ", ceiling,
			U"\nPull formants = ", pullFormants);
	try {
		integer maxnCandidates = Pitch_getMaxnCandidates (me);
		integer place;
		volatile double maximum, value;
		double ceiling2 = pullFormants ? 2.0 * ceiling : ceiling;
		/* Next three lines 20011015 */
		double timeStepCorrection = 0.01 / my dx;
		octaveJumpCost *= timeStepCorrection;
		voicedUnvoicedCost *= timeStepCorrection;

		my ceiling = ceiling;
		autoNUMmatrix <double> delta (1, my nx, 1, maxnCandidates);
		autoNUMmatrix <integer> psi (1, my nx, 1, maxnCandidates);

		for (integer iframe = 1; iframe <= my nx; iframe ++) {
			Pitch_Frame frame = & my frame [iframe];
			double unvoicedStrength = silenceThreshold <= 0 ? 0.0 :
				2.0 - frame -> intensity / (silenceThreshold / (1.0 + voicingThreshold));
			unvoicedStrength = voicingThreshold + (unvoicedStrength > 0.0 ? unvoicedStrength : 0.0);
			for (integer icand = 1; icand <= frame -> nCandidates; icand ++) {
				Pitch_Candidate candidate = & frame -> candidate [icand];
				bool voiceless = ! Pitch_util_frequencyIsVoiced (candidate -> frequency, ceiling2);
				delta [iframe] [icand] = voiceless ? unvoicedStrength :
					candidate -> strength - octaveCost * NUMlog2 (ceiling / candidate -> frequency);
			}
		}

		/* Look for the most probable path through the maxima. */
		/* There is a cost for the voiced/unvoiced transition, */
		/* and a cost for a frequency jump. */

		for (integer iframe = 2; iframe <= my nx; iframe ++) {
			Pitch_Frame prevFrame = & my frame [iframe - 1], curFrame = & my frame [iframe];
			double *prevDelta = delta [iframe - 1], *curDelta = delta [iframe];
			integer *curPsi = psi [iframe];
			for (integer icand2 = 1; icand2 <= curFrame -> nCandidates; icand2 ++) {
				double f2 = curFrame -> candidate [icand2]. frequency;
				maximum = -1e30;
				place = 0;
				for (integer icand1 = 1; icand1 <= prevFrame -> nCandidates; icand1 ++) {
					double f1 = prevFrame -> candidate [icand1]. frequency;
					double transitionCost;
					bool previousVoiceless = ! Pitch_util_frequencyIsVoiced (f1, ceiling2);
					bool currentVoiceless = ! Pitch_util_frequencyIsVoiced (f2, ceiling2);
					if (currentVoiceless) {
						if (previousVoiceless) {
							transitionCost = 0.0;   // both voiceless
						} else {
							transitionCost = voicedUnvoicedCost;   // voiced-to-unvoiced transition
						}
					} else {
						if (previousVoiceless) {
							transitionCost = voicedUnvoicedCost;   // unvoiced-to-voiced transition
							if (Melder_debug == 30) {
								/*
									Try to take into account a frequency jump across a voiceless stretch.
								*/
								integer place1 = icand1;
								for (integer jframe = iframe - 2; jframe >= 1; jframe --) {
									place1 = psi [jframe + 1] [place1];
									f1 = my frame [jframe]. candidate [place1]. frequency;
									if (Pitch_util_frequencyIsVoiced (f1, ceiling)) {
										transitionCost += octaveJumpCost * fabs (NUMlog2 (f1 / f2)) / (iframe - jframe);
										break;
									}
								}
							}
						} else {
							transitionCost = octaveJumpCost * fabs (NUMlog2 (f1 / f2));   // both voiced
						}
					}
					value = prevDelta [icand1] - transitionCost + curDelta [icand2];
					//if (Melder_debug == 33) Melder_casual (U"Frame ", iframe, U", current candidate ", icand2,
					//  U" (delta ", curDelta [icand2], U"), previous candidate ", icand1, U" (delta ", prevDelta [icand1], U"), ",
					//	U"transition cost ", transitionCost, U", value ", value, U", maximum ", maximum);
					if (value > maximum) {
						maximum = value;
						place = icand1;
					} else if (value == maximum) {
						if (Melder_debug == 33)
							Melder_casual (
								U"A tie in frame ", iframe,
								U", current candidate ", icand2,
								U", previous candidate ", icand1
							);
					}
				}
				curDelta [icand2] = maximum;
				curPsi [icand2] = place;
			}
		}

		/* Find the end of the most probable path. */

		place = 1;
		maximum = delta [my nx] [place];
		for (integer icand = 2; icand <= my frame [my nx]. nCandidates; icand ++) {
			if (delta [my nx] [icand] > maximum) {
				place = icand;
				maximum = delta [my nx] [place];
			}
		}

		/* Backtracking: follow the path backwards. */

		for (integer iframe = my nx; iframe >= 1; iframe --) {
			if (Melder_debug == 33)
				Melder_casual (
					U"Frame ", iframe, U":",
					U" swapping candidates 1 and ", place
				);
			Pitch_Frame frame = & my frame [iframe];
			structPitch_Candidate help = frame -> candidate [1];
			frame -> candidate [1] = frame -> candidate [place];
			frame -> candidate [place] = help;
			place = psi [iframe] [place];   // This assignment is challenging to CodeWarrior 11.
		}

		/* Pull formants: devoice frames with frequencies between ceiling and ceiling2. */

		if (ceiling2 > ceiling) {
			if (Melder_debug == 33)
				Melder_casual (U"Pulling formants...");
			for (integer iframe = my nx; iframe >= 1; iframe --) {
				Pitch_Frame frame = & my frame [iframe];
				Pitch_Candidate winner = & frame -> candidate [1];
				double f = winner -> frequency;
				if (f > ceiling && f < ceiling2) {
					for (integer icand = 2; icand <= frame -> nCandidates; icand ++) {
						Pitch_Candidate loser = & frame -> candidate [icand];
						if (loser -> frequency == 0.0) {
							structPitch_Candidate help = * winner;
							* winner = * loser;
							* loser = help;
							break;
						}
					}
				}
			}
		}
	} catch (MelderError) {
		Melder_throw (me, U": path not found.");
	}
}

void Pitch_drawInside (Pitch me, Graphics g, double xmin, double xmax, double fmin, double fmax, bool speckle, kPitch_unit unit) {
	Sampled_drawInside (me, g, xmin, xmax, fmin, fmax, speckle, Pitch_LEVEL_FREQUENCY, (int) unit);
}

void Pitch_draw (Pitch me, Graphics g, double tmin, double tmax, double fmin, double fmax, bool garnish, bool speckle, kPitch_unit unit) {
	Graphics_setInner (g);
	Pitch_drawInside (me, g, tmin, tmax, fmin, fmax, speckle, unit);
	Graphics_unsetInner (g);
	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, true, U"Time (s)");
		Graphics_marksBottom (g, 2, true, true, false);
		Graphics_textLeft (g, true, Melder_cat (U"Pitch (", Function_getUnitText (me, Pitch_LEVEL_FREQUENCY, (int) unit, Function_UNIT_TEXT_GRAPHICAL), U")"));
		if (Function_isUnitLogarithmic (me, Pitch_LEVEL_FREQUENCY, (int) unit)) {
			Graphics_marksLeftLogarithmic (g, 6, true, true, false);
		} else {
			Graphics_marksLeft (g, 2, true, true, false);
		}
	}
}

void Pitch_difference (Pitch me, Pitch thee) {
	integer nuvtov = 0, nvtouv = 0, ndfdown = 0, ndfup = 0;
	if (my nx != thy nx || my dx != thy dx || my x1 != thy x1) {
		Melder_flushError (U"Pitch_difference: these Pitches are not aligned.");
		return;
	}
	for (integer i = 1; i <= my nx; i ++) {
		double myf = my frame [i]. candidate [1]. frequency, thyf = thy frame [i]. candidate [1]. frequency;
		int myUnvoiced = ! Pitch_util_frequencyIsVoiced (myf, my ceiling);
		int thyUnvoiced = ! Pitch_util_frequencyIsVoiced (thyf, thy ceiling);
		double t = Sampled_indexToX (me, i);
		if (myUnvoiced && ! thyUnvoiced) {
			Melder_casual (
				U"Frame ", i,
				U" time ", t,
				U": unvoiced to voiced."
			);
			nuvtov ++;
		} else if (! myUnvoiced && thyUnvoiced) {
			Melder_casual (
				U"Frame ", i,
				U" time ", t,
				U": voiced to unvoiced."
			);
			nvtouv ++;
		} else if (! myUnvoiced && ! thyUnvoiced) {
			if (myf > thyf) {
				//Melder_casual ("Frame %ld time %f: downward frequency jump from %.5g Hz to %.5g Hz.", i, t, myf, thyf);
				ndfdown ++;
			} else if (myf < thyf) {
				//Melder_casual ("Frame %ld time %f: upward frequency jump from %.5g Hz to %.5g Hz.", i, t, myf, thyf);
				ndfup ++;
			}
		}
	}
	MelderInfo_open ();
	MelderInfo_writeLine (U"Difference between two Pitches:");
	MelderInfo_writeLine (U"Unvoiced to voiced: ", nuvtov, U" frames.");
	MelderInfo_writeLine (U"Voiced to unvoiced: ", nvtouv, U" frames.");
	MelderInfo_writeLine (U"Downward frequency jump: ", ndfdown, U" frames.");
	MelderInfo_writeLine (U"Upward frequency jump: ", ndfup, U" frames.");
	MelderInfo_close ();
}

autoPitch Pitch_killOctaveJumps (Pitch me) {
	try {
		autoPitch thee = Pitch_create (my xmin, my xmax, my nx, my dx, my x1, my ceiling, 2);
		integer nVoiced = 0, nUp = 0;
		double lastFrequency = 0.0;
		for (integer i = 1; i <= my nx; i ++) {
			double f = my frame [i]. candidate [1]. frequency;
			thy frame [i]. candidate [1]. strength = my frame [i]. candidate [1]. strength;
			if (Pitch_util_frequencyIsVoiced (f, my ceiling)) {
				nVoiced ++;
				if (lastFrequency != 0.0) {
					double fmin = lastFrequency * 0.7071, fmax = 2.0 * fmin;
					while (f < fmin) { f *= 2.0; nUp ++; }
					while (f > fmax) { f *= 0.5; nUp --; }
				}
				lastFrequency = thy frame [i]. candidate [1]. frequency = f;
			}
		}
		thy ceiling *= 2.0;   // make room for some octave jumps
		while (nUp > nVoiced / 2) {
			for (integer i = 1; i <= thy nx; i ++)
				thy frame [i]. candidate [1]. frequency *= 0.5;
			nUp -= nVoiced;
		}
		while (nUp < - nVoiced / 2) {
			for (integer i = 1; i <= thy nx; i ++)
				thy frame [i]. candidate [1]. frequency *= 2.0;
			nUp += nVoiced;
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": octave jumps not killed.");
	}
}

autoPitch Pitch_interpolate (Pitch me) {
	try {
		autoPitch thee = Pitch_create (my xmin, my xmax, my nx, my dx, my x1, my ceiling, 2);
		for (integer i = 1; i <= my nx; i ++) {
			double f = my frame [i]. candidate [1]. frequency;
			thy frame [i]. candidate [1]. strength = 0.9;
			if (Pitch_util_frequencyIsVoiced (f, my ceiling)) {
				thy frame [i]. candidate [1]. frequency = f;
			} else {
				integer left, right;
				double fleft = 0.0, fright = 0.0;
				for (left = i - 1; left >= 1 && fleft == 0.0; left --) {
					fleft = my frame [left]. candidate [1]. frequency;
					if (fleft >= my ceiling) fleft = 0.0;
				}
				for (right = i + 1; right <= my nx && fright == 0.0; right ++) {
					fright = my frame [right]. candidate [1]. frequency;
					if (fright >= my ceiling) fright = 0.0;
				}
				if (fleft != 0.0 && fright != 0.0)
					thy frame [i]. candidate [1]. frequency =
						((i - left) * fright + (right - i) * fleft) / (right - left);
			}
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not interpolated.");
	}
}

autoPitch Pitch_subtractLinearFit (Pitch me, kPitch_unit unit) {
	try {
		autoPitch thee = Pitch_interpolate (me);
		/*
			Find the first and last voiced frame.
		*/
		integer imin = thy nx + 1, imax = 0;
		for (integer i = 1; i <= my nx; i ++)
			if (Pitch_isVoiced_i (thee.get(), i)) { imin = i; break; }
		for (integer i = imin + 1; i <= my nx; i ++)
			if (! Pitch_isVoiced_i (thee.get(), i)) { imax = i - 1; break; }
		integer n = imax - imin + 1;
		if (n <= 1) {
			/*
				If there is only one point, there is no slope to subtract.
				The result is the original frequency if the original was voiced, else 0.0 Hz.
				The interpolated Pitch already has these properties set correctly.
			*/
			return thee;
		}
		/*
			Compute average pitch and time.
		*/
		double sum = 0.0;
		for (integer i = imin; i <= imax; i ++) {
			sum += Sampled_getValueAtSample (thee.get(), i, Pitch_LEVEL_FREQUENCY, (int) unit);
		}
		double fmean = sum / n;
		double tmean = thy x1 + (0.5 * (imin + imax) - 1) * thy dx;
		/*
			Compute slope.
		*/
		double numerator = 0.0, denominator = 0.0;
		for (integer i = imin; i <= imax; i ++) {
			double t = thy x1 + (i - 1) * thy dx - tmean;
			double f = Sampled_getValueAtSample (thee.get(), i, Pitch_LEVEL_FREQUENCY, (int) unit) - fmean;
			numerator += f * t;
			denominator += t * t;
		}
		double slope = numerator / denominator;
		/*
			Modify frequencies.
		*/
		for (integer i = imin; i <= imax; i ++) {
			Pitch_Frame myFrame = & my frame [i], thyFrame = & thy frame [i];
			Pitch_Candidate myCandidate = & myFrame -> candidate [1], thyCandidate = & thyFrame -> candidate [1];
			double t = thy x1 + (i - 1) * thy dx - tmean, myFreq = myCandidate -> frequency;
			double f = Sampled_getValueAtSample (thee.get(), i, Pitch_LEVEL_FREQUENCY, (int) unit);
			f -= slope * t;
			if (Pitch_util_frequencyIsVoiced (myFreq, my ceiling)) {
				thyCandidate -> frequency = Function_convertSpecialToStandardUnit (me, f, Pitch_LEVEL_FREQUENCY, (int) unit);
			} else {
				thyCandidate -> frequency = 0.0;
			}
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": linear fit not subtracted.");
	}
}

autoPitch Pitch_smooth (Pitch me, double bandWidth) {
	try {
		autoPitch interp = Pitch_interpolate (me);
		autoMatrix matrix1 = Pitch_to_Matrix (interp.get());
		double logicalDuration = matrix1 -> xmax - matrix1 -> xmin;
		double physicalDuration = matrix1 -> nx * matrix1 -> dx;
		autoSound sound1 = Sound_create (
			1,   // one channel
			matrix1 -> xmin - logicalDuration,   // triply wide logical time domain (ignored)
			matrix1 -> xmax + logicalDuration,   // triply wide logical time domain (ignored)
			3 * matrix1 -> nx,   // triply wide physical time domain
			matrix1 -> dx,
			matrix1 -> x1 - physicalDuration   // (ignored)
		);

		integer firstVoiced = 0, lastVoiced = 0;
		for (integer i = 1; i <= matrix1 -> nx; i ++) {
			double f = matrix1 -> z [1] [i];
			if (f != 0.0) {
				if (! firstVoiced) firstVoiced = i;
				lastVoiced = i;
				sound1 -> z [1] [i + matrix1 -> nx] = f;
			}
		}

		/* Extrapolate. */
		double fextrap = matrix1 -> z [1] [firstVoiced];
		firstVoiced += matrix1 -> nx;
		for (integer i = 1; i < firstVoiced; i ++)
			sound1 -> z [1] [i] = fextrap;
		fextrap = matrix1 -> z [1] [lastVoiced];
		lastVoiced += matrix1 -> nx;
		for (integer i = lastVoiced + 1; i <= sound1 -> nx; i ++)
			sound1 -> z [1] [i] = fextrap;

		/* Smooth. */
		autoSpectrum spectrum = Sound_to_Spectrum (sound1.get(), true);
		for (integer i = 1; i <= spectrum -> nx; i ++) {
			double f = (i - 1) * spectrum -> dx, fT = f / bandWidth, factor = exp (- fT * fT);
			spectrum -> z [1] [i] *= factor;
			spectrum -> z [2] [i] *= factor;
		}
		autoSound sound2 = Spectrum_to_Sound (spectrum.get());

		autoMatrix matrix2 = Matrix_create (my xmin, my xmax, my nx, my dx, my x1, 1.0, 1.0, 1, 1.0, 1.0);
		for (integer i = 1; i <= my nx; i ++) {
			double originalF0 = my frame [i]. candidate [1]. frequency;
			matrix2 -> z [1] [i] = ( Pitch_util_frequencyIsVoiced (originalF0, my ceiling) ?
				sound2 -> z [1] [i + matrix2 -> nx] : 0.0 );
		}
		autoPitch thee = Matrix_to_Pitch (matrix2.get());
		thy ceiling = my ceiling;
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not smoothed.");
	}
}

void Pitch_step (Pitch me, double step, double precision, double tmin, double tmax) {
	Melder_assert (precision >= 0.0 && precision < 1.0);
	integer imin, imax;
	if (! Sampled_getWindowSamples (me, tmin, tmax, & imin, & imax)) return;
	for (integer i = imin; i <= imax; i ++) {
		Pitch_Frame frame = & my frame [i];
		double currentFrequency = frame -> candidate [1]. frequency;
		if (Pitch_util_frequencyIsVoiced (currentFrequency, my ceiling)) {
			double targetFrequency = currentFrequency * step;
			double fmin = (1.0 - precision) * targetFrequency;
			double fmax = (1.0 + precision) * targetFrequency;
			int nearestCandidate = 0;
			double nearestDistance = my ceiling;
			if (fmax > my ceiling) fmax = my ceiling;
			for (int icand = 2; icand <= frame -> nCandidates; icand ++) {
				double f = frame -> candidate [icand]. frequency;
				if (f > fmin && f < fmax) {
					double localDistance = fabs (f - targetFrequency);
					if (localDistance < nearestDistance) {
						nearestCandidate = icand;
						nearestDistance = localDistance;
					}
				}
			}
			if (nearestCandidate) {   // swap candidates
				struct structPitch_Candidate candidate = frame -> candidate [nearestCandidate];
				frame -> candidate [nearestCandidate] = frame -> candidate [1];
				frame -> candidate [1] = candidate;
			}
		}
	}
}

/* End of file Pitch.cpp */
