/* Sound_and_Spectrum.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/05/21
 * pb 2002/07/16 GPL
 * pb 2003/03/09 shorter sounds from Hann band filtering
 * pb 2003/05/15 replaced memcof with NUMburg
 */

#include "Sound_and_Spectrum.h"
#include "NUM2.h"

Spectrum Sound_to_Spectrum_fft (Sound me) {
	Spectrum thee = NULL;
	long nfft = 2, i, halfnfft;
	float *data = NULL, scaling, *re, *im;
	while (nfft < my nx) nfft *= 2;
	halfnfft = nfft / 2;
	data = NUMfvector (1, nfft); cherror
	for (i = 1; i <= my nx; i ++) data [i] = my z [1] [i];
	NUMrealft (data, nfft, 1);
	thee = Spectrum_create (0.5 / my dx, halfnfft + 1); cherror
	re = thy z [1]; im = thy z [2];
	scaling = my dx;
	re [1] = data [1] * scaling;
	for (i = 2; i <= halfnfft; i ++) {
		re [i] = data [i + i - 1] * scaling;
		im [i] = data [i + i] * scaling;
	}
	re [halfnfft + 1] = data [2] * scaling;
end:
	NUMfvector_free (data, 1);
	iferror forget (thee);
	return thee;
}

Sound Spectrum_to_Sound_fft (Spectrum me) {
	Sound thee = NULL;
	long nfft, i, halfnfft = 1;
	float *amp, scaling, *re = my z [1], *im = my z [2];
	if (my x1 != 0.0) {
		Melder_error ("(Spectrum_to_Sound:) An FFTable Spectrum must have a first frequency of 0 Hz, not %.8g Hz.", my x1);
		goto end;
	}
	while (halfnfft < my nx - 1) halfnfft *= 2;
	nfft = halfnfft * 2;
	thee = Sound_createSimple (1 / my dx, nfft * my dx); cherror
	amp = thy z [1];
	scaling = 2 * my dx;
	for (i = 1; i < my nx; i ++) {
		amp [i + i - 1] = re [i] * scaling;
		amp [i + i] = im [i] * scaling;
	}
	if (halfnfft == my nx - 1)   /* Probably produced by forward FFT. */
		amp [2] = re [my nx] * scaling;
	else {
		i = my nx;
		amp [i + i - 1] = re [i] * scaling;
		amp [i + i] = im [i] * scaling;
	}
	NUMrealft (amp, nfft, -1);
end:
	iferror forget (thee);
	return thee;
}

Spectrum Spectrum_lpcSmoothing (Spectrum me, int numberOfPeaks, double preemphasisFrequency) {
	Sound sound = NULL;
	Spectrum thee = NULL;
	float gain, a [100], *data = NULL, *re, *im;
	long i, numberOfCoefficients = 2 * numberOfPeaks, nfft, halfnfft, ndata;
	double scale;

	sound = Spectrum_to_Sound_fft (me); cherror
	NUMpreemphasize_f (sound -> z [1], sound -> nx, sound -> dx, preemphasisFrequency);	 	
	
	NUMburg (sound -> z [1], sound -> nx, a, numberOfCoefficients, & gain);
	for (i = 1; i <= numberOfCoefficients; i ++) a [i] = - a [i];
	thee = Data_copy (me); cherror

	nfft = 2 * (thy nx - 1);
	ndata = numberOfCoefficients < nfft ? numberOfCoefficients : nfft - 1;
	scale = 10 * (gain > 0 ? sqrt (gain) : 1) / numberOfCoefficients;
	data = NUMfvector (1, nfft); cherror
	data [1] = 1;
	for (i = 1; i <= ndata; i ++)
		data [i + 1] = a [i];
	NUMrealft (data, nfft, 1);
	re = thy z [1], im = thy z [2];
	re [1] = scale / data [1];
	im [1] = 0.0;
	halfnfft = nfft / 2;
	for (i = 2; i <= halfnfft; i ++) {
		double real = data [i + i - 1], imag = data [i + i];
		re [i] = scale / sqrt (real * real + imag * imag) / (1 + thy dx * (i - 1) / preemphasisFrequency);
		im [i] = 0;
	}
	re [halfnfft + 1] = scale / data [2] / (1 + thy dx * halfnfft / preemphasisFrequency);
	im [halfnfft + 1] = 0.0;
end:
	forget (sound);
	NUMfvector_free (data, 1);
	iferror forget (thee);
	return thee;
}

Sound Sound_filter_formula (Sound me, const char *formula) {
	Sound filtered = NULL, result = NULL;
	Spectrum spec;
	long i;
	spec = Sound_to_Spectrum_fft (me); cherror
	Matrix_formula ((Matrix) spec, formula, 0); cherror
	filtered = Spectrum_to_Sound_fft (spec); cherror
	/* The filtered signal may be much longer than the original, so: */
	result = Data_copy (me); cherror
	for (i = 1; i <= my nx; i ++)
		result -> z [1] [i] = filtered -> z [1] [i];
end:
	forget (filtered);
	forget (spec);
	iferror forget (result);
	return result;
}

Sound Sound_filter_passHannBand (Sound me, double fmin, double fmax, double smooth) {
	Spectrum spec = NULL;
	Sound thee = Data_copy (me), him;
	spec = Sound_to_Spectrum_fft (me); cherror
	Spectrum_passHannBand (spec, fmin, fmax, smooth);
	him = Spectrum_to_Sound_fft (spec); cherror
	NUMfvector_copyElements (his z [1], thy z [1], 1, thy nx);
end:
	forget (spec);
	forget (him);
	iferror forget (thee);
	return thee;
}

Sound Sound_filter_stopHannBand (Sound me, double fmin, double fmax, double smooth) {
	Spectrum spec = NULL;
	Sound thee = Data_copy (me), him;
	spec = Sound_to_Spectrum_fft (me); cherror
	Spectrum_stopHannBand (spec, fmin, fmax, smooth);
	him = Spectrum_to_Sound_fft (spec); cherror
end:
	NUMfvector_copyElements (his z [1], thy z [1], 1, thy nx);
	forget (spec);
	forget (him);
	iferror forget (thee);
	return thee;
}

/* End of file Sound_and_Spectrum.c */