/*
 *  Copyright (c) 2007 - 2008 by Damien Di Fede <ddf@compartmental.net>
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU Library 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 Library General Public License for more details.
 *
 *   You should have received a copy of the GNU Library General Public
 *   License along with this program; if not, write to the Free Software
 *   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

// This Chebyshev Filter implementation has been ported from the BASIC 
// implementation outlined in Chapter 20 of The Scientist and Engineer's
// Guide to Signal Processing, which can be found at:
//
//     http://www.dspguide.com/ch20.htm

package ddf.minim.effects;

import ddf.minim.Minim;

/**
 * A Chebyshev filter is an IIR filter that uses a particular method to
 * calculate the coefficients of the filter. It is defined by whether it is a
 * low pass filter or a high pass filter and the number of poles it has. You
 * needn't worry about what a pole is, exactly, just know that more poles
 * usually makes for a better filter. An additional limitation is that the
 * number of poles must be even. See {@link #setPoles(int)} for more information
 * about poles. Another characteristic of Chebyshev filters is how much "ripple"
 * they allow in the pass band. The pass band is the range of frequencies that
 * the filter lets through. The "ripple" in the pass band can be seen as wavy
 * line in the frequency response of the filter. Lots of ripple is bad, but more
 * ripple gives a faster rolloff from the pass band to the stop band (the range
 * of frequencies blocked by the filter). Faster rolloff is good because it
 * means the cutoff is sharper. Ripple is expressed as a percentage, such as
 * 0.5% ripple.
 * 
 * @author Damien Di Fede
 * @see <a href="http://www.dspguide.com/ch20.htm">Chebyshev Filters</a>
 * 
 */
public class ChebFilter extends IIRFilter
{
  /** A constant used to indicate a low pass filter. */
  public static final int LP = 1;
  /** A constant used to indicate a high pass filter. */
  public static final int HP = 2;
  private static final float PI = (float) Math.PI;
  private int type, poles;
  private float ripple;
  private boolean canCalc;

  /**
   * Constructs a Chebyshev filter with a cutoff of the given frequency, of the given
   * type, with the give amount of ripple in the pass band, and with the given
   * number of poles, that will be used to filter audio of that was recorded at
   * the given sample rate.
   * 
   * @param frequency
   *          the cutoff frequency of the filter
   * @param type
   *          the type of filter, either ChebFilter.LP or ChebFilter.HP
   * @param ripple
   *          the percentage of ripple, such as 0.005
   * @param poles
   *          the number of poles, must be even and in the range [2, 20]
   * @param sampleRate
   *          the sample rate of audio that will be filtered
   */
  public ChebFilter(float frequency, int type, float ripple, int poles,
      float sampleRate)
  {
    super(frequency, sampleRate);
    canCalc = false;
    setType(type);
    setRipple(ripple);
    setPoles(poles);
    canCalc = true;
    calcCoeff();
    initArrays();
  }

  /**
   * Sets the type of the filter. Either ChebFilter.LP or ChebFilter.HP
   * 
   * @param t
   *          the type of the filter
   */
  public void setType(int t)
  {
    if ( t != LP && t != HP )
    {
      Minim.error("Invalid filter type, defaulting to low pass.");
      t = LP;
    }
    type = t;
    calcCoeff();
  }

  /**
   * Returns the type of the filter.
   */
  public int getType()
  {
    return type;
  }

  /**
   * Sets the ripple percentage of the filter.
   * 
   * @param r
   *          the ripple percentage
   */
  public void setRipple(float r)
  {
    ripple = r;
    calcCoeff();
  }

  /**
   * Returns the ripple percentage of the filter.
   * 
   * @return the ripple percentage
   */
  public float getRipple()
  {
    return ripple;
  }

  /**
   * Sets the number of poles used in the filter. The number of poles must be
   * even and between 2 and 20. This function will report an error if either of
   * those conditions are not met. However, it should also be mentioned that
   * depending on the current cutoff frequency of the filter, the number of
   * poles that will result in a <i>stable</i> filter, can be a few as 4. The
   * function does not report an error in the case of the number of requested
   * poles resulting in an unstable filter. For reference, here is a table of
   * the maximum number of poles possible according to cutoff frequency:
   * <p>
   * <table border="1" cellpadding="5">
   * <tr>
   * <td>Cutoff Frequency<br />
   * (expressed as a fraction of the sampling rate)</td>
   * <td>0.02</td>
   * <td>0.05</td>
   * <td>0.10</td>
   * <td>0.25</td>
   * <td>0.40</td>
   * <td>0.45</td>
   * <td>0.48</td>
   * </tr>
   * <tr>
   * <td>Maximum poles</td>
   * <td>4</td>
   * <td>6</td>
   * <td>10</td>
   * <td>20</td>
   * <td>10</td>
   * <td>6</td>
   * <td>4</td>
   * </tr>
   * </table>
   * 
   * @param p -
   *          the number of poles
   */
  public void setPoles(int p)
  {
    if (p < 2)
    {
      Minim
          .error("ChebFilter.setPoles: The number of poles must be at least 2.");
      return;
    }
    if (p % 2 != 0)
    {
      Minim.error("ChebFilter.setPoles: The number of poles must be even.");
      return;
    }
    if (p > 20)
    {
      Minim.error("ChebFilter.setPoles: The maximum number of poles is 20.");
    }
    poles = p;
    calcCoeff();
    initArrays();
  }

  /**
   * Returns the number of poles in the filter.
   * 
   * @return the number of poles
   */
  public int getPoles()
  {
    return poles;
  }

  protected void calcCoeff()
  {
    a = new float[0];
    b = new float[0];
    if (!canCalc) return;
    // where the poles will wind up
    float[] ca = new float[23];
    float[] cb = new float[23];
    // temporary arrays for working with ca and cb
    float[] ta = new float[23];
    float[] tb = new float[23];
    // arrays to hold the two-pole coefficients
    // used during the aggregation process
    float[] pa = new float[3];
    float[] pb = new float[2];
    // I don't know why this must be done
    ca[2] = 1;
    cb[2] = 1;
    // calculate two poles at a time
    for (int p = 1; p <= poles / 2; p++)
    {
      // calc pair p, put the results in pa and pb
      calcTwoPole(p, pa, pb);
      // copy ca and cb into ta and tb
      System.arraycopy(ca, 0, ta, 0, ta.length);
      System.arraycopy(cb, 0, tb, 0, tb.length);
      // add coefficients to the cascade
      for (int i = 2; i < 23; i++)
      {
        ca[i] = pa[0] * ta[i] + pa[1] * ta[i - 1] + pa[2] * ta[i - 2];
        cb[i] = tb[i] - pb[0] * tb[i - 1] - pb[1] * tb[i - 2];
      }
    }
    // final stage of combining coefficients
    cb[2] = 0;
    for (int i = 0; i < 21; i++)
    {
      ca[i] = ca[i + 2];
      cb[i] = -cb[i + 2];
    }
    // normalize the gain
    float sa = 0;
    float sb = 0;
    for (int i = 0; i < 21; i++)
    {
      if (type == LP)
      {
        sa += ca[i];
        sb += cb[i];
      }
      else
      {
        sa += ca[i] * (float) Math.pow(-1, i);
        sb += cb[i] * (float) Math.pow(-1, i);
      }
    }
    float gain = sa / (1 - sb);
    for (int i = 0; i < 21; i++)
    {
      ca[i] /= gain;
    }
    // initialize the coefficient arrays used by process()
    a = new float[poles + 1];
    b = new float[poles];
    // copy the values from ca and cb into a and b
    // in this implementation cb[0] = 0 and cb[1] is where
    // the b coefficients begin, so they are numbered the way
    // one normally numbers coefficients when talking about IIR filters
    // however, process() expects b[0] to be the coefficient B1
    // so we copy cb over to b starting at index 1
    System.arraycopy(ca, 0, a, 0, a.length);
    System.arraycopy(cb, 1, b, 0, b.length);
  }

  private void calcTwoPole(int p, float[] pa, float[] pb)
  {
    float np = (float) poles;
    float angle = PI / (np * 2) + (p - 1) * PI / np;
    float rp = -(float) Math.cos(angle);
    float ip = (float) Math.sin(angle);
    // warp from a circle to an ellipse
    if (ripple > 0)
    {
      float es = 1 / (float) Math
          .sqrt(Math.pow(100.0f / (100.0f - ripple), 2) - 1);
      float vx = (1 / np) * (float) Math.log(es + Math.sqrt(es * es + 1));
      float kx = (1 / np) * (float) Math.log(es + Math.sqrt(es * es - 1));
      kx = (float) (Math.exp(kx) + Math.exp(-kx)) / 2;
      rp *= ((Math.exp(vx) - Math.exp(-vx)) / 2) / kx;
      ip *= ((Math.exp(vx) + Math.exp(-vx)) / 2) / kx;
    }
    // s-domain to z-domain conversion
    float t = 2 * (float) Math.tan(0.5);
    float w = 2 * PI * (frequency()/sampleRate());
    float m = rp * rp + ip * ip;
    float d = 4 - 4 * rp * t + m * t * t;
    float x0 = (t * t) / d;
    float x1 = (2 * t * t) / d;
    float x2 = (t * t) / d;
    float y1 = (8 - 2 * m * t * t) / d;
    float y2 = (-4 - 4 * rp * t - m * t * t) / d;
    // LP to LP, or LP to HP transform
    float k;
    if (type == HP)
      k = -(float) Math.cos(w / 2 + 0.5f) / (float) Math.cos(w / 2 - 0.5f);
    else
      k = (float) Math.sin(0.5f - w / 2) / (float) Math.sin(0.5 + w / 2);
    d = 1 + y1 * k - y2 * k * k;
    pa[0] = (x0 - x1 * k + x2 * k * k) / d;
    pa[1] = (-2 * x0 * k + x1 + x1 * k * k - 2 * x2 * k) / d;
    pa[2] = (x0 * k * k - x1 * k + x2) / d;
    pb[0] = (2 * k + y1 + y1 * k * k - 2 * y2 * k) / d;
    pb[1] = (-(k * k) - y1 * k + y2) / d;
    if (type == HP)
    {
      pa[1] = -pa[1];
      pb[0] = -pb[0];
    }
  }
}
