/* Sound_to_Spectrogram8.c */
/* Paul Boersma, June 5, 1995 */

#include "Sound_to_Spectrogram8.h"

Spectrogram8 Sound_to_Spectrogram8 (Sound me, double analysisWidth, double fmax,
	double timeOversampling, double frequencyOversampling, int windowType)
{
	Spectrogram8 thee = NULL;
	double timeWidth = analysisWidth / sqrt (NUMpi);
	double frequencyWidth = 1 / timeWidth, scaling;
	double timeStep = timeWidth / timeOversampling, freqStep = frequencyWidth / frequencyOversampling;
	long nBands = (long) floor (fmax / freqStep);
	long nsamp_window, halfnsamp_window, nFrames, nsampFFT = 1, half_nsampFFT;
	long iframe, iband, i, j;
	float *amplitude = my z [1], *frame, *window, *oneByBinWidth;
	long *lofreqi;
	double duration = my dx * (double) my nx, t1;
	double actualAnalysisWidth = analysisWidth;
	if (nBands < 1) return NULL;
	if (windowType == 1) actualAnalysisWidth *= 2.0;   /* Gaussian window. */
	nsamp_window = (long) floor (actualAnalysisWidth / my dx);
	halfnsamp_window = nsamp_window / 2 - 1;
	nsamp_window = halfnsamp_window * 2;
	if (actualAnalysisWidth > duration) return NULL;
	nFrames = 1 + (long) floor ((duration - actualAnalysisWidth) / timeStep);   /* >= 1 */
	t1 = my x1 + 0.5 * ((double) (my nx - 1) * my dx - (double) (nFrames - 1) * timeStep);
		/* Centre of first frame. */

	if (! (thee = Spectrogram8_create (my xmin, my xmax, nFrames, timeStep, t1,
			0.0, fmax, nBands, freqStep, 0.5 * freqStep)))
		return NULL;

	nsampFFT = 1; while (nsampFFT < nsamp_window || nsampFFT < 2 * nBands) nsampFFT *= 2;
	half_nsampFFT = nsampFFT / 2;

	if (! (frame = NUMfvector (1, nsampFFT)) ||
		 ! (window = NUMfvector (1, nsamp_window)) ||
		 ! (lofreqi = NUMlvector (1, nBands + 1)) ||
		 ! (oneByBinWidth = NUMfvector (1, nBands))) return NULL;

	for (i = 1; i <= nsamp_window; i ++) {
		double nSamplesPerWindow_f = actualAnalysisWidth / my dx;
		double phase = (double) i / nSamplesPerWindow_f;   /* 0 .. 1 */
		double value;
		switch (windowType) {
			case 1: /* Gaussian. */
			{
				double imid = 0.5 * (double) (nsamp_window + 1), edge = exp (-12.0);
				phase = ((double) i - imid) / nSamplesPerWindow_f;   /* -0.5 .. +0.5 */
				value = (exp (-48.0 * phase * phase) - edge) / (1.0 - edge);
				break;
			}
			case 2: /* Hamming (raised sine-squared). */
				value = 0.54 - 0.46 * cos (2.0 * NUMpi * phase); break;
			case 3: /* Hanning (sine-squared). */
				value = 0.5 * (1.0 - cos (2.0 * NUMpi * phase)); break;
			case 4: /* Welch (parabolic). */
				value = 1.0 - (2.0 * phase - 1.0) * (2.0 * phase - 1.0); break;
			case 5: /* Bartlett (triangular). */
				value = 1.0 - fabs ((2.0 * phase - 1.0)); break;
			case 6: /* None (rectangular). */
				value = 1.0; break;
		}
		window [i] = (float) value;
	}
	/* Binning frame [1..nBands]. */
	for (iband = 1; iband <= nBands + 1; iband ++) {
		double lofreq = (double) (iband - 1) * freqStep;
		lofreqi [iband] = (long) floor (lofreq * my dx * (double) nsampFFT) + 1;
		if (lofreqi [iband] > half_nsampFFT) lofreqi [iband] = half_nsampFFT;
	}
	scaling = my dx * my dx   /* See Sound_and_Spectrum.c */
		* analysisWidth / 4.0e-10;   /* Apx. twice mean duration; see Spectrum_draw.c */
	for (iband = 1; iband <= nBands; iband ++) {
		oneByBinWidth [iband] = scaling *
			(float) lofreqi [iband + 1] / nsampFFT * 1000 * my dx /   /* Preemphasis. */
			(float) (lofreqi [iband + 1] - lofreqi [iband] + 1) *   /* Undersampling correction. */
				(freqStep * my dx * nsampFFT);
	}
	Melder_progress (0.0, "Sound to Spectrogram8...");   /* 20010216 */
	for (iframe = 1; iframe <= nFrames; iframe ++) {
		double t = Sampled_indexToX (thee, iframe);
		float tmp;
		long leftSample = Sampled_xToLowIndex (me, t), rightSample = leftSample + 1;
		long startSample = rightSample - halfnsamp_window;
		long endSample = leftSample + halfnsamp_window;
		Melder_assert (startSample >= 1);
		Melder_assert (endSample <= my nx);
		if (iframe % 10 == 1 && ! Melder_progress ((double) iframe / nFrames,
			"From Sound to Spectrogram8: %.4g / %.4g seconds", t, my xmax)) break;
		for (j = 1, i = startSample; j <= nsamp_window; j ++)
			frame [j] = amplitude [i ++] * window [j];
		for (j = nsamp_window + 1; j <= nsampFFT; j ++) frame [j] = 0.0f;

		/* Compute Fast Fourier Transform of the frame. */

		NUMrealft (frame, nsampFFT, 1);   /* Complex spectrum. */

		/* Put power spectrum in frame [1..half_nsampFFT + 1]. */

		frame [1] *= frame [1];   /* DC component. */
		tmp = frame [2];   /* Nyquist frequency. */
		for (i = 2; i <= half_nsampFFT; i ++)
			frame [i] = frame [i + i - 1] * frame [i + i - 1] + frame [i + i] * frame [i + i];
		frame [half_nsampFFT + 1] = tmp * tmp;   /* ? */

		/* Bin into frame [1..nBands]. */
		for (iband = 1; iband <= nBands; iband ++) {
			double power = 0.0;
			for (i = lofreqi [iband]; i <= lofreqi [iband + 1]; i ++) power += frame [i];
			power = 20 * log10 (power * oneByBinWidth [iband] + 1e-30) + 200;
			thy z [iband] [iframe] = power < 0 ? 0 : power > 255 ? 255 : power;
		}
	}
	Melder_progress (1.0, NULL);
	NUMfvector_free (frame, 1, nsampFFT);
	NUMfvector_free (window, 1, nsamp_window);
	NUMlvector_free (lofreqi, 1, nBands + 1);
	NUMfvector_free (oneByBinWidth, 1, nBands);
	return thee;
}

/* End of file Sound_and_Spectrogram8.c */
