/* Libvisual - The audio visualisation framework.
 *
 * Copyright (C) 2004, 2005, 2006 Dennis Smit <ds@nerds-incorporated.org>
 *
 * The FFT implementation found in this file is based upon the NULLSOFT
 * Milkdrop FFT implementation.
 *
 * Authors: Dennis Smit <ds@nerds-incorporated.org>
 *          Chong Kai Xiong <descender@phreaker.net>
 *
 * $Id: lv_fourier.c,v 1.15 2006/02/13 20:54:08 synap Exp $
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation; either version 2.1
 * 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

#define FOURIER_PI 3.141592653589793238462643383279502884197169399f

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <string.h>

#include "lv_cache.h"
#include "lv_utils.h"
#include "lv_math.h"
#include "lv_fourier.h"

/* Log scale settings */
#define AMP_LOG_SCALE_THRESHOLD0	0.001f
#define AMP_LOG_SCALE_DIVISOR		6.908f	/* divisor = -log threshold */
#define FREQ_LOG_SCALE_BASE		2.0f

#define DFT_CACHE_ENTRY(obj)				(VISUAL_CHECK_CAST ((obj), DFTCacheEntry))
#define LOG_SCALE_CACHE_ENTRY(obj)			(VISUAL_CHECK_CAST ((obj), LogScaleCacheEntry))

typedef struct _DFTCacheEntry DFTCacheEntry;
typedef struct _LogScaleCacheEntry LogScaleCacheEntry;

struct _DFTCacheEntry {
	VisObject	 object;

	int		 spectrum_size;

	float		*bitrevtable;

	float		*sintable;
	float		*costable;
};

struct _LogScaleCacheEntry {
	VisObject	 object;

	float		*range;
};

static VisCache __lv_dft_cache;
static VisCache __lv_log_scale_cache;
static int __lv_fourier_initialized = FALSE;


static int dft_dtor (VisObject *object);

static void fft_table_bitrev_init (DFTCacheEntry *fcache, VisDFT *fourier);
static void fft_table_cossin_init (DFTCacheEntry *fcache, VisDFT *fourier);
static void dft_table_cossin_init (DFTCacheEntry *fcache, VisDFT *fourier);
static void range_table_init (LogScaleCacheEntry *lcache, int size);

static int dft_cache_destroyer (VisObject *object);
static DFTCacheEntry *dft_cache_get (VisDFT *dft);

static int log_scale_cache_destroy (VisObject *object);
static LogScaleCacheEntry *log_scale_cache_get (int size);

static void perform_dft_brute_force (VisDFT *fourier, float *output, float *input);
static void perform_fft_radix2_dit (VisDFT *fourier, float *output, float *input);

static int dft_dtor (VisObject *object)
{
	VisDFT *dft = VISUAL_DFT (object);

	/* FIXME not all object elements are freed, destroyed, whatever */

	if (dft->real != NULL)
		visual_mem_free (dft->real);

	if (dft->imag != NULL)
		visual_mem_free (dft->imag);

	dft->real = NULL;
	dft->imag = NULL;

	return VISUAL_OK;
}

static void fft_table_bitrev_init (DFTCacheEntry *fcache, VisDFT *fourier)
{
	unsigned int i, m, temp;
	unsigned int j = 0;

	fcache->bitrevtable = visual_mem_malloc0 (sizeof (unsigned int) * fourier->spectrum_size);

	for (i = 0; i < fourier->spectrum_size; i++)
		fcache->bitrevtable[i] = i;

	for (i = 0; i < fourier->spectrum_size; i++) {
		if (j > i) {
			temp = fcache->bitrevtable[i];
			fcache->bitrevtable[i] = fcache->bitrevtable[j];
			fcache->bitrevtable[j] = temp;
		}

		m = fourier->spectrum_size >> 1;

		while (m >= 1 && j >= m) {
			j -= m;
			m >>= 1;
		}

		j += m;
	}
}

static void fft_table_cossin_init (DFTCacheEntry *fcache, VisDFT *fourier)
{
	unsigned int i, dftsize, tabsize;
	float theta;

	dftsize = 2;
	tabsize = 0;
	while (dftsize <= fourier->spectrum_size) {
		tabsize++;

		dftsize <<= 1;
	}

	fcache->sintable = visual_mem_malloc0 (sizeof (float) * tabsize);
	fcache->costable = visual_mem_malloc0 (sizeof (float) * tabsize);

	dftsize = 2;
	i = 0;
	while (dftsize <= fourier->spectrum_size) {
		theta = (float) (-2.0f * FOURIER_PI / (float) dftsize);

		fcache->costable[i] = (float) cosf (theta);
		fcache->sintable[i] = (float) sinf (theta);

		i++;

		dftsize <<= 1;
	}
}

static void dft_table_cossin_init (DFTCacheEntry *fcache, VisDFT *fourier)
{
	unsigned int i, tabsize;
	float theta;

	tabsize = fourier->spectrum_size / 2;
	fcache->sintable = visual_mem_malloc0 (sizeof (float) * tabsize);
	fcache->costable = visual_mem_malloc0 (sizeof (float) * tabsize);

	for (i = 0; i < tabsize; i++) {
		theta = (-2.0f * FOURIER_PI * i) / fourier->spectrum_size;

		fcache->costable[i] = cosf (theta);
		fcache->sintable[i] = sinf (theta);
	}
}

static void range_table_init (LogScaleCacheEntry *lcache, int size)
{
	unsigned int i, tabsize;
	float factor, factor_scale;

	tabsize = size;// / 2 + 1;

	lcache->range = visual_mem_malloc0 (sizeof (float) * tabsize);

	factor = 1.0f;
	factor_scale = 1.0f / FREQ_LOG_SCALE_BASE;

	for (i = tabsize; i >= 1; i--)
	{
		lcache->range[i] = (unsigned int) (factor * tabsize);
		factor *= factor_scale;
	}

	/* cache->range[0] is 0.0 */
}

static int dft_cache_destroyer (VisObject *object)
{
	DFTCacheEntry *fcache = DFT_CACHE_ENTRY (object);

	if (fcache->bitrevtable != NULL)
		visual_mem_free (fcache->bitrevtable);

	if (fcache->sintable != NULL)
		visual_mem_free (fcache->sintable);

	if (fcache->costable != NULL)
		visual_mem_free (fcache->costable);

	fcache->bitrevtable = NULL;
	fcache->sintable = NULL;
	fcache->costable = NULL;

	return VISUAL_OK;
}

static DFTCacheEntry *dft_cache_get (VisDFT *fourier)
{
	DFTCacheEntry *fcache;
	char key[16];

	visual_log_return_val_if_fail (__lv_fourier_initialized == TRUE, NULL);

	snprintf (key, 16, "%d", fourier->spectrum_size);
	fcache = visual_cache_get (&__lv_dft_cache, key);

	if (fcache == NULL) {
		fcache = visual_mem_new0 (DFTCacheEntry, 1);

		visual_object_initialize (VISUAL_OBJECT (fcache), TRUE, dft_cache_destroyer);

		if (fourier->brute_force) {
			dft_table_cossin_init (fcache, fourier);
		} else {
			fft_table_bitrev_init (fcache, fourier);
			fft_table_cossin_init (fcache, fourier);
		}

		visual_cache_put (&__lv_dft_cache, key, fcache);
	}

	return fcache;
}

static int log_scale_cache_destroyer (VisObject *object)
{
	LogScaleCacheEntry *lcache = LOG_SCALE_CACHE_ENTRY (object);

	if (lcache->range != NULL)
		visual_mem_free (lcache->range);

	lcache->range = NULL;

	return VISUAL_OK;
}

static LogScaleCacheEntry *log_scale_cache_get (int size)
{
	LogScaleCacheEntry *lcache;
	char key[16];

	visual_log_return_val_if_fail (__lv_fourier_initialized == TRUE, NULL);

	snprintf (key, 16, "%d", size);
	lcache = visual_cache_get (&__lv_log_scale_cache, key);

	if (lcache == NULL) {
		lcache = visual_mem_new0 (LogScaleCacheEntry, 1);

		visual_object_initialize (VISUAL_OBJECT (lcache), TRUE, log_scale_cache_destroyer);

		range_table_init (lcache, size);

		visual_cache_put (&__lv_log_scale_cache, key, lcache);
	}

	return lcache;
}


/**
 * @defgroup VisDFT VisDFT
 * @{
 */

int visual_fourier_initialize ()
{
	visual_cache_init (&__lv_dft_cache, visual_object_collection_destroyer, 50, NULL, TRUE);
	visual_cache_init (&__lv_log_scale_cache, visual_object_collection_destroyer, 50, NULL, TRUE);

	__lv_fourier_initialized = TRUE;

	return VISUAL_OK;
}

int visual_fourier_is_initialized ()
{
	return __lv_fourier_initialized;
}

int visual_fourier_deinitialize ()
{
	if (__lv_fourier_initialized == FALSE)
		return -VISUAL_ERROR_FOURIER_NOT_INITIALIZED;

	visual_object_unref (VISUAL_OBJECT (&__lv_dft_cache));
	visual_object_unref (VISUAL_OBJECT (&__lv_log_scale_cache));

	__lv_fourier_initialized = FALSE;

	return VISUAL_OK;
}

/**
 * Function to create a new VisDFT Discrete Fourier Transform context used
 * to calculate amplitude spectrums over audio data.
 *
 * \note For optimal performance, use a power-of-2 spectrum size. The
 * current implementation does not use the Fast Fourier Transform for
 * non powers of 2.
 *
 * \note If samples_in is smaller than 2 * samples_out, the input will be padded
 * with zeroes.
 *
 * @param samples_in The number of samples provided to every call to
 * visual_dft_perform() as input.
 * @param samples_out Size of output spectrum (number of output samples).
 *
 * @return A newly created VisDFT.
 */
VisDFT *visual_dft_new (unsigned int samples_out, unsigned int samples_in)
{
	VisDFT *dft;

	dft = visual_mem_new0 (VisDFT, 1);

	visual_dft_init (dft, samples_in, samples_out);

	/* Do the VisObject initialization */
	visual_object_set_allocated (VISUAL_OBJECT (dft), TRUE);
	visual_object_ref (VISUAL_OBJECT (dft));

	return dft;
}

int visual_dft_init (VisDFT *dft, unsigned int samples_out, unsigned int samples_in)
{
	visual_log_return_val_if_fail (dft != NULL, -VISUAL_ERROR_FOURIER_NULL);

	/* Do the VisObject initialization */
	visual_object_clear (VISUAL_OBJECT (dft));
	visual_object_set_dtor (VISUAL_OBJECT (dft), dft_dtor);
	visual_object_set_allocated (VISUAL_OBJECT (dft), FALSE);

	/* Set the VisDFT data */
	dft->samples_in = samples_in;
	dft->spectrum_size = samples_out * 2;
	dft->brute_force = !visual_utils_is_power_of_2 (dft->spectrum_size);

	/* Initialize the VisDFT */
	dft_cache_get (dft);

	dft->real = visual_mem_malloc0 (sizeof (float) * dft->spectrum_size);
	dft->imag = visual_mem_malloc0 (sizeof (float) * dft->spectrum_size);

	return VISUAL_OK;
}

static void perform_dft_brute_force (VisDFT *dft, float *output, float *input)
{
	DFTCacheEntry *fcache;
	unsigned int i, j;
	float xr, xi, wr, wi, wtemp;

	fcache = dft_cache_get (dft);
	visual_object_ref (VISUAL_OBJECT (fcache));

	for (i = 0; i < dft->spectrum_size / 2; i++) {
		xr = 0.0f;
		xi = 0.0f;

		wr = 1.0f;
		wi = 0.0f;

		for (j = 0; j < dft->spectrum_size; j++) {
			xr += input[j] * wr;
			xi += input[j] * wi;

			wtemp = wr;
			wr = wr * fcache->costable[i] - wi * fcache->sintable[i];
			wi = wtemp * fcache->sintable[i] + wi * fcache->costable[i];
		}

		dft->real[i] = xr;
		dft->imag[i] = xi;
	}

	visual_object_unref (VISUAL_OBJECT (fcache));
}

static void perform_fft_radix2_dit (VisDFT *dft, float *output, float *input)
{
	DFTCacheEntry *fcache;
	unsigned int j, m, i, dftsize, hdftsize, t;
	float wr, wi, wpi, wpr, wtemp, tempr, tempi;

	fcache = dft_cache_get (dft);
	visual_object_ref (VISUAL_OBJECT (fcache));

	for (i = 0; i < dft->spectrum_size; i++) {
		unsigned int idx = fcache->bitrevtable[i];

		if (idx < dft->samples_in)
			dft->real[i] = input[idx];
		else
			dft->real[i] = 0;
	}

	visual_mem_set (dft->imag, 0, sizeof (float) * dft->spectrum_size);

	dftsize = 2;
	t = 0;
	while (dftsize <= dft->spectrum_size) {
		wpr = fcache->costable[t];
		wpi = fcache->sintable[t];

		wr = 1.0f;
		wi = 0.0f;

		hdftsize = dftsize >> 1;

		for (m = 0; m < hdftsize; m += 1) {
			for (i = m; i < dft->spectrum_size; i+=dftsize) {
				j = i + hdftsize;

				tempr = wr * dft->real[j] - wi * dft->imag[j];
				tempi = wr * dft->imag[j] + wi * dft->real[j];

				dft->real[j] = dft->real[i] - tempr;
				dft->imag[j] = dft->imag[i] - tempi;

				dft->real[i] += tempr;
				dft->imag[i] += tempi;
			}

			wr = (wtemp = wr) * wpr - wi * wpi;
			wi = wi * wpr + wtemp * wpi;
		}

		dftsize <<= 1;
		t++;
	}

	visual_object_unref (VISUAL_OBJECT (fcache));
}


/**
 * Function to perform a Discrete Fourier Transform over a set of data.
 *
 * \note Output samples are normalised to [0.0, 1.0] by dividing with the
 * spectrum size.
 *
 * @param dft Pointer to the VisDFT context for this transform.
 * @param output Array of output samples
 * @param input Array of input samples with values in [-1.0, 1.0]
 *
 * @return VISUAL_OK on succes, -VISUAL_ERROR_FOURIER_NULL or -VISUAL_ERROR_NULL on failure.
 */
int visual_dft_perform (VisDFT *dft, float *output, float *input)
{
	unsigned int i;

	visual_log_return_val_if_fail (dft != NULL, -VISUAL_ERROR_FOURIER_NULL);
	visual_log_return_val_if_fail (output != NULL, -VISUAL_ERROR_NULL);
	visual_log_return_val_if_fail (input != NULL, -VISUAL_ERROR_NULL);

	if (dft->brute_force)
		perform_dft_brute_force (dft, output, input);
	else
		perform_fft_radix2_dit (dft, output, input);

	visual_math_vectorized_complex_to_norm_scale (output, dft->real, dft->imag,
			dft->spectrum_size / 2,
			1.0 / dft->spectrum_size);

	return VISUAL_OK;
}

/**
 * Function to scale an ampltitude spectrum logarithmically.
 *
 * \note Scaled values are guaranteed to be in [0.0, 1.0].
 *
 * @param output Array of output samples
 * @param input  Array of input samples with values in [0.0, 1.0]
 * @param size Array size.
 *
 * @Return VISUAL_OK on success, VISUAL_ERROR_NULL on failure.
 */
int visual_dft_log_scale (float *output, float *input, int size)
{
	LogScaleCacheEntry *lcache;
	unsigned int i, j;
	float amp;

	visual_log_return_val_if_fail (output != NULL, -VISUAL_ERROR_NULL);
	visual_log_return_val_if_fail (input != NULL, -VISUAL_ERROR_NULL);

	return visual_dft_log_scale_standard (output, input, size);

	lcache = log_scale_cache_get (size);
	visual_object_ref (VISUAL_OBJECT (lcache));

	for (i = 0; i < size; i++) {
		amp = 0.0f;

		for (j = lcache->range[i]; j < lcache->range[i+1]; j++) {
			if (amp < input[j])
				amp = input[j];
		}

		if (amp > AMP_LOG_SCALE_THRESHOLD0)
			amp = 1.0f + log (input[i]) / AMP_LOG_SCALE_DIVISOR;
		else
			amp = 0.0f;

		output[i] = amp;
	}

	visual_object_unref (VISUAL_OBJECT (lcache));

	return VISUAL_OK;
}

int visual_dft_log_scale_standard (float *output, float *input, int size)
{
	int i;

	visual_log_return_val_if_fail (output != NULL, -VISUAL_ERROR_NULL);
	visual_log_return_val_if_fail (input != NULL, -VISUAL_ERROR_NULL);

	return visual_dft_log_scale_custom (output, input, size, AMP_LOG_SCALE_DIVISOR);
}

int visual_dft_log_scale_custom (float *output, float *input, int size, float log_scale_divisor)
{
	int i;

	visual_log_return_val_if_fail (output != NULL, -VISUAL_ERROR_NULL);
	visual_log_return_val_if_fail (input != NULL, -VISUAL_ERROR_NULL);

	for (i = 0; i < size; i++) {
		if (input[i] > AMP_LOG_SCALE_THRESHOLD0)
			output[i] = 1.0f + log (input[i]) / log_scale_divisor;
		else
			output[i] = 0.0f;
	}

	return VISUAL_OK;
}

/**
 * @}
 */