/* BSE Resampling Datahandles
 * Copyright (C) 2001-2003 Tim Janik
 * Copyright (C) 2006 Stefan Westerfeld
 *
 * This library 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 of the License, or (at your option) any later version.
 *
 * This library 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 library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/* 
 * halfband FIR filter for factor 2 resampling, created with octave
 *
 * design method: windowed sinc,  using ultraspherical window
 *
 *   coefficients = 32
 *             x0 = 1.01065
 *          alpha = 0.75
 *
 * design criteria (44100 Hz => 88200 Hz):
 * 
 *       passband = [     0, 18000 ]  1 - 2^-16 <= H(z) <= 1+2^-16
 *     transition = [ 18000, 26100 ]
 *       stopband = [ 26100, 44100 ]  | H(z) | <= -96 dB
 *
 * and for 48 kHz => 96 kHz:
 *
 *       passband = [     0, 19589 ]  1 - 2^-16 <= H(z) <= 1+2^-16
 *     transition = [ 19588, 29386 ]
 *       stopband = [ 29386, 48000 ]  | H(z) | <= -96 dB
 *
 * in order to keep the coefficient number down to 32, the filter
 * does only "almost" fulfill the spec, but its really really close
 * (no stopband ripple > -95 dB)
 */

const double halfband_fir_upsample2_96db_coeffs[32] =
{
  -3.48616530828033e-05,
   0.000112877490936198,
  -0.000278961878372482,
   0.000590495306376081,
  -0.00112566995029848,
   0.00198635062559427,
  -0.00330178798332932,
   0.00523534239035401,
  -0.00799905465189065,
   0.0118867161189188,
  -0.0173508611368417,
   0.0251928452706978,
  -0.0370909694665106,
   0.057408291607388,
  -0.102239638342325,
   0.317002929635456,
   /* here, a 0.5 coefficient will be used */
   0.317002929635456,
  -0.102239638342325,
   0.0574082916073878,
  -0.0370909694665105,
   0.0251928452706976,
  -0.0173508611368415,
   0.0118867161189186,
  -0.00799905465189052,
   0.0052353423903539,
  -0.00330178798332923,
   0.00198635062559419,
  -0.00112566995029842,
   0.000590495306376034,
  -0.00027896187837245,
   0.000112877490936177,
  -3.48616530827983e-05
};

/*
 *   coefficients = 16
 *             x0 = 1.013
 *          alpha = 0.2
 */
const double halfband_fir_upsample2_48db_coeffs[16] =
{
  -0.00270578824181636,
   0.00566964586625895,
  -0.0106460585587187,
   0.0185209590435965,
  -0.0310433957594089,
   0.0525722488176905,
  -0.0991138314110143,
   0.315921760444802,
   /* here, a 0.5 coefficient will be used */
   0.315921760444802,
  -0.0991138314110145,
   0.0525722488176907,
  -0.031043395759409,
   0.0185209590435966,
  -0.0106460585587187,
   0.00566964586625899,
  -0.00270578824181638
};

#include "bsedatahandle-resample.h"
#include "bseresampler.hh"
#include "gsldatahandle.h"
#include <sfi/sficxx.hh>
#include <vector>

namespace Bse
{

namespace Resampler
{

class Resampler2Factory {
public:
  enum Mode {
    MODE_UPSAMPLE,
    MODE_DOWNSAMPLE
  };
  enum Precision {
    PREC_48DB = 8,
    PREC_96DB = 16
  };
  /**
   * creates a resampler
   */
  static Resampler2* create (Mode mode, Precision precision);

  /**
   * creates a filter with a specific coefficient sets
   */
  template<class Filter>
  static Resampler2* create (const double *d, guint order, double scaling);
};

template<class Filter> Resampler2*
Resampler2Factory::create (const double *d, guint order, double scaling)
{
  float taps[order];
  for (guint i = 0; i < order; i++)
    taps[i] = d[i] * scaling;

  Filter *filter = new Filter (taps);
  g_assert (order == filter->order());
  return filter;
}

Resampler2*
Resampler2Factory::create (Mode mode, Precision precision)
{
  if (mode == MODE_UPSAMPLE)
    {
      switch (precision)
	{
	case PREC_48DB: return create< Upsampler2<16,1> > (halfband_fir_upsample2_48db_coeffs, 16, 2.0);
	case PREC_96DB: return create< Upsampler2<32,1> > (halfband_fir_upsample2_96db_coeffs, 32, 2.0);
	}
    }
  return 0;
}

}

using std::vector;
using Resampler::Resampler2;
using Resampler::Resampler2Factory;

class DataHandleUpsample2 : public GslDataHandle,
			    public Sfi::GNewable /* 0 initialization */
{
  GslDataHandle	       *src_handle;
  int                   precision_bits;

public:
  gboolean		init_ok;
  vector<Resampler2 *>  upsamplers;
  int64			pcm_frame;
  vector<float>		pcm_data;
  int64			frame_size;
  int64                 filter_delay;
  int64                 filter_order;

  DataHandleUpsample2 (GslDataHandle *src_handle, int precision_bits)
    : src_handle (src_handle),
      precision_bits (precision_bits),
      init_ok (false)
  {
    g_return_if_fail (src_handle != NULL);
  
    init_ok = gsl_data_handle_common_init (this, NULL);
    if (init_ok)
      {
	name = g_strconcat (src_handle->name, "// #upsample2 /", NULL);
	src_handle = gsl_data_handle_ref (src_handle);
      }
  }

  ~DataHandleUpsample2()
  {
    gsl_data_handle_unref (src_handle);
    gsl_data_handle_common_free (this);
  }

  BseErrorType
  open (GslDataHandleSetup *setup)
  {
    BseErrorType error = gsl_data_handle_open (src_handle);
    if (error != BSE_ERROR_NONE)
      return error;

    *setup = src_handle->setup; /* copies setup.xinfos by pointer */
    setup->mix_freq *= 2.0;
    setup->n_values *= 2;

    frame_size = 1024 * setup->n_channels;
    pcm_frame = -2;
    pcm_data.resize (frame_size);

    for (guint i = 0; i < setup->n_channels; i++)
      {
	Resampler2Factory::Precision precision = static_cast<Resampler2Factory::Precision> (precision_bits);
	Resampler2 *resampler = Resampler2Factory::create (Resampler2Factory::MODE_UPSAMPLE, precision);
	g_assert (resampler); /* FIXME: better error handling */

	upsamplers.push_back (resampler);
	filter_order = resampler->order();

	g_assert (filter_order % 2 == 0);
	filter_delay = filter_order / 2 + 1;
      }
    return BSE_ERROR_NONE;
  }

  void
  close()
  {
    for (guint i = 0; i < setup.n_channels; i++)
      delete upsamplers[i];

    upsamplers.clear();
    pcm_data.clear();
    
    setup.xinfos = NULL;	/* cleanup pointer reference */
    gsl_data_handle_close (src_handle);
  }

  int64
  src_read (int64   voffset,
	    int64   n_values,
	    gfloat *values)
  {
    voffset += filter_delay * setup.n_channels;	/* compensate filter delay */

    int64 left = n_values;
    do
      {
	int64 l;
	if (voffset >= 0 && voffset < src_handle->setup.n_values)
	  {
	    l = gsl_data_handle_read (src_handle, voffset, std::min (left, src_handle->setup.n_values - voffset), values);
	    if (l < 0)
	      return l;	/* pass on errors */
	  }
	else
	  {
	    /*
	     * the resampler needs data before and after the data provided
	     * by src_handle; so we offer zero values here
	     */
	    *values = 0;
	    l = 1;
	  }

	voffset += l;
	left -= l;
	values += l;
      }
    while (left > 0);

    return n_values;
  }

  void
  deinterleave (float* src, float *dest, int64 n_values)
  {
    const int64 n_channels = setup.n_channels;

    for (int64 ch = 0; ch < n_channels; ch++)
      for (int64 v = ch; v < n_values; v += n_channels)
	*dest++ = src[v];
  }

  void
  interleave (float* src, float *dest, int64 n_values)
  {
    const int64 n_channels = setup.n_channels;

    for (int64 ch = 0; ch < n_channels; ch++)
      for (int64 v = ch; v < n_values; v += n_channels)
	dest[v] = *src++;
  }

  int64
  prepare_filter_history (int64 frame)
  {
    const int64 n_channels = setup.n_channels;
    const int64 n_input_samples = filter_order;

    float input_interleaved[n_input_samples * n_channels];
    float input[n_input_samples * n_channels];

    int64 l = src_read (frame * frame_size / 2 - n_input_samples * n_channels,
	                n_input_samples * n_channels, input_interleaved);
    if (l < 0)
      return l; /* pass on errors */

    deinterleave (input_interleaved, input, n_input_samples * setup.n_channels);

    for (guint ch = 0; ch < setup.n_channels; ch++)
      {
	/* we don't need the output, this is just for filling the filter history */
	float output[n_input_samples * 2];

	upsamplers[ch]->process_block (input + ch * n_input_samples, n_input_samples, output);
      }
    return 1;
  }

  int64
  read_frame (int64 frame)
  {
    /*
     * if we're seeking (not reading data linearily), we need to reinitialize
     * the filter history with new values
     */
    if (frame != pcm_frame + 1)
      {
	int64 l = prepare_filter_history (frame);
	if (l < 0)
	  return l; /* pass on errors */
      }

    float input_interleaved[frame_size / 2];
    float input[frame_size / 2];
    float output[frame_size];

    int64 l = src_read (frame * frame_size / 2, frame_size / 2, input_interleaved);
    if (l < 0)
      return l; /* pass on errors */

    deinterleave (input_interleaved, input, frame_size / 2);
    for (guint ch = 0; ch < setup.n_channels; ch++)
      {
	const int64 output_per_channel = frame_size / setup.n_channels;
	const int64 input_per_channel = output_per_channel / 2;

	upsamplers[ch]->process_block (input + ch * input_per_channel, input_per_channel, output + ch * output_per_channel);
      }
    interleave (output, &pcm_data[0], frame_size);

    pcm_frame = frame;
    return 1;
  }

  int64
  read (int64          voffset,
	int64          n_values,
	gfloat        *values)
  {
    int64 frame = voffset / pcm_data.size();
    if (frame != pcm_frame)
      {
	int64 l = read_frame (frame);
	if (l < 0)
	  return l;
      }
    g_assert (pcm_frame == frame);

    voffset -= pcm_frame * frame_size;
    g_assert (voffset >= 0);

    n_values = std::min (n_values, frame_size - voffset);
    for (int64 i = 0; i < n_values; i++)
      values[i] = pcm_data[voffset + i];

    return n_values;
  }
};

} // namespace Bse

using namespace Bse;

static BseErrorType
dh_upsample2_open (GslDataHandle *dhandle, GslDataHandleSetup *setup)
{
  return reinterpret_cast<DataHandleUpsample2 *> (dhandle)->open (setup);
}

static void
dh_upsample2_close (GslDataHandle *dhandle)
{
  reinterpret_cast<DataHandleUpsample2 *> (dhandle)->close();
}


static void
dh_upsample2_destroy (GslDataHandle *dhandle)
{
  delete reinterpret_cast<DataHandleUpsample2 *> (dhandle);
}

static int64
dh_upsample2_read (GslDataHandle *dhandle,
		   int64          voffset,
		   int64          n_values,
		   gfloat        *values)
{
  return reinterpret_cast<DataHandleUpsample2 *> (dhandle)->read (voffset, n_values, values);
}

GslDataHandle*
bse_data_handle_new_upsample2 (GslDataHandle *src_handle, int precision_bits)
{
  static GslDataHandleFuncs dh_upsample2_vtable = {
    dh_upsample2_open,
    dh_upsample2_read,
    dh_upsample2_close,
    NULL,
    dh_upsample2_destroy,
  };
  DataHandleUpsample2 *dhandle = new DataHandleUpsample2 (src_handle, precision_bits);
  if (dhandle->init_ok)
    {
      dhandle->vtable = &dh_upsample2_vtable;
      return dhandle;
    }
  else
    {
      delete dhandle;
      return NULL;
    }
}
