/* Wavelet.c
 *
 * Copyright (C) 1992-2002 David Weenink & 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 1997/05/30
 * dw 2002/06/20 modified selection of k-th smallest (removed NUMselect)
 * pb 2002/07/16 GPL
 */

#include "Wavelet.h"

#include "oo_DESTROY.h"
#include "Wavelet_def.h"
#include "oo_COPY.h"
#include "Wavelet_def.h"
#include "oo_EQUAL.h"
#include "Wavelet_def.h"
#include "oo_WRITE_ASCII.h"
#include "Wavelet_def.h"
#include "oo_READ_ASCII.h"
#include "Wavelet_def.h"
#include "oo_WRITE_BINARY.h"
#include "Wavelet_def.h"
#include "oo_READ_BINARY.h"
#include "Wavelet_def.h"
#include "oo_DESCRIPTION.h"
#include "Wavelet_def.h"

static void info (I) {
	iam (Wavelet);
	Melder_information (
		"Length of Wavelet: %ld\n"
		"Highest frequency component: %.8f Hz\n"
		"Number of filter coefficients: %d",
			my n, my highestFrequency, my ncof); 
}

class_methods (Wavelet, Data)
	class_method_local (Wavelet, destroy)
	class_method_local (Wavelet, copy)
	class_method_local (Wavelet, equal)
	class_method_local (Wavelet, readAscii)
	class_method_local (Wavelet, writeAscii)
	class_method_local (Wavelet, readBinary)
	class_method_local (Wavelet, writeBinary)
	class_method_local (Wavelet, description)
	class_method (info)
class_methods_end

Wavelet Wavelet_create (long n, double highestFrequency, long ncof) {
	Wavelet me = new (Wavelet);
	if (! me || ! (my a = NUMfvector (1, n))) { forget (me); return NULL; }
	my n = n;
	my highestFrequency = highestFrequency;
	my ncof = ncof;
	return me;
}

Wavelet Wavelet_truncateByFraction (Wavelet me, double fraction) {
	Wavelet thee;
	float *buffer, level;
	long i;
	if (! (buffer = NUMfvector (1, my n)) ||
		 ! (thee = Data_copy (me)))
		return Melder_errorp ("Wavelet_truncateByFraction: out of memory.");
	for (i = 1; i <= my n; i ++) buffer [i] = fabs (my a [i]);
	NUMsort_f (my n, buffer);
	i = (int)(fraction * my n);
	level = buffer [i < 1 ? 1 : i];
/*	level = NUMselect ((int) ((1.0 - fraction) * my n), my n, buffer); */
	NUMfvector_free (buffer, 1);
	for (i = 1; i <= my n; i ++) if (fabs (thy a [i]) < level) thy a [i] = 0.0;
	return thee;          
}

Wavelet Wavelet_truncateByLevel (Wavelet me, double level) {
	Wavelet thee;
	long i;
	if (! (thee = Data_copy (me)))
		return Melder_errorp ("Wavelet_truncateByLevel: not enough memory.");
	for (i = 1; i <= my n; i ++) if (fabs (thy a [i]) <= level) thy a[i] = 0.0;
	return thee;          
}

double Wavelet_fractionBelowLevel (Wavelet me, double level) {
	long i, nbelow = 0;
	for (i = 1; i <= my n; i ++) if (fabs (my a [i]) < level) nbelow ++;
	return ((double) nbelow) / my n;
}

void Wavelet_draw (Wavelet me, Graphics graphics,
	double tmin, double tmax, double fmin, double fmax)
{
	long itime, itmin, itmax, ifreq, ifmin, ifmax;
	long offset, npsd, nguess, numberOfFrequencies;
	double dflog, dflin;
	double dt, t = my n / (2.0 * my highestFrequency); 
	float *tmp, bias, weight, x90, x10, grey = 0.8;

	if (tmax <= tmin) { tmin = 0; tmax = t; }
	if (fmax <= fmin) { fmin = 0; fmax = my highestFrequency; }
	/* For linear frequency scale: 
	 *     df == FN/(2^m-1),
	 *     FN == Nyquist frequency 
	 *     m  == number of frequency bands (scales)
	 *     i-th frequency band starts at (2^(i-1)-1)*df
	 *     width of i-th band is 2^(i-1)*df
	 * For log frequency scale:
	 *     df == FN/ m;
	 *     i-th frequency band starts at (i-1)*df
	 *     width of i-th band is df
	 */
	/* Calculate lowest and highest frequency to display */
	numberOfFrequencies = floor (log (my n - 2) / log (2.0));  
	dflin = my highestFrequency / ((1 << numberOfFrequencies) - 1.0);
	dflog = my highestFrequency / numberOfFrequencies;
	ifmin = 1; ifmax = numberOfFrequencies;
	for (ifreq = 1; ifreq <= numberOfFrequencies; ifreq ++)
	{
		double fb = ((1 << (ifreq - 1)) - 1) * dflin;
		double fe = ((1 << ifreq) - 1) * dflin;
		if (fmin >= fb && fmin <  fe ) ifmin = ifreq;
		if (fmax >  fb && fmax <= fe ) ifmax = ifreq;
	} 
	/* temporary storage for finding 10% and 90% values */
	nguess = 2 * (1 << ifmax);
	if (nguess > my n) nguess = my n;
	tmp = NUMfvector (1, nguess);
	if (! tmp) { Melder_flushError (NULL); return; }     
	for (npsd = 0, ifreq = ifmin; ifreq <= ifmax; ifreq ++) {   /* 0 --> 1; ppgb 19970530 */
		dt = t / (1 << ifreq); 
		offset = (1 << ifreq);
		itmin = offset + 1 + floor (tmin / dt);
		itmax = offset + floor (tmax / dt);
		for (itime = itmin; itime <= itmax; itime ++) {
			Melder_assert (itime <= my n);
			tmp [++npsd] = my a [itime];
		}
	}
	NUMsort_f (npsd, tmp);
	x90 = tmp [(int)(0.9 * npsd)];
	x10 = tmp [(int)(0.1 * npsd)];
	/*x90 = NUMselect (0.95 * npsd, npsd, tmp);
	x10 = NUMselect (0.1 * npsd / 0.95, 0.95 * npsd, tmp);*/
	NUMfvector_free (tmp, 1);
	if (x10 != x90) {
		bias = - (x10 + x90) / (x10 - x90) * log ((1.0 - grey) / grey);
		weight = 2.0 / (x10 - x90) * log ((1.0 - grey) / grey);
	}
	Graphics_setInner (graphics);
	Graphics_setWindow (graphics, tmin, tmax, fmin, fmax);
	for (ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
		double t1WC, t2WC, f1WC, f2WC;
		dt = t / (1 << ifreq); 
		offset = (1 << ifreq);
		itmin = 1 + floor (tmin / dt);
		itmax = floor (tmax / dt);
		f1WC = (ifreq - 1) * dflog;
		f2WC = f1WC + dflog;
		if (f1WC < fmin) f1WC = fmin;
		if (f2WC > fmax) f2WC = fmax;
		for (itime = itmin; itime <= itmax; itime ++) {
			t1WC = (itime - 1) * dt;
			t2WC = t1WC + dt;
			if (t1WC < tmin) t1WC = tmin;
			if (t2WC > tmax) t2WC = tmax;
			Graphics_setGrey (graphics, NUMsigmoid (- weight * my a [offset + itime] - bias));
			Graphics_fillRectangle (graphics, t1WC, t2WC, f1WC, f2WC);
		}
	}
	Graphics_setGrey (graphics, 0);
	Graphics_unsetInner (graphics);

	Graphics_drawInnerBox (graphics);
	Graphics_textBottom (graphics, 1, "Time (s)");
	Graphics_marksBottom (graphics, 2, 1, 1, 0);
	Graphics_textLeft (graphics, 1, "Log frequency (Hz)");
	Graphics_marksLeft (graphics, 2, 1, 1, 0);    
}

/* End of file Wavelet.c */