/* Sound_and_Spectrum.c */
/* Paul Boersma, May 21, 2002 */

#include "Sound_and_Spectrum.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);	 	
	
	NUMmemcof (sound -> z [1], sound -> nx, numberOfCoefficients, & gain, a);
	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;
	Sound thee = NULL;
	spec = Sound_to_Spectrum_fft (me); cherror
	Spectrum_passHannBand (spec, fmin, fmax, smooth);
	thee = Spectrum_to_Sound_fft (spec); cherror
end:
	forget (spec);
	iferror forget (thee);
	return thee;
}

Sound Sound_filter_stopHannBand (Sound me, double fmin, double fmax, double smooth) {
	Spectrum spec;
	Sound thee = NULL;
	spec = Sound_to_Spectrum_fft (me); cherror
	Spectrum_stopHannBand (spec, fmin, fmax, smooth);
	thee = Spectrum_to_Sound_fft (spec); cherror
end:
	forget (spec);
	iferror forget (thee);
	return thee;
}

/* End of file Sound_and_Spectrum.c */
