engines/gig/Filter.h

/***************************************************************************
 *                                                                         *
 *   LinuxSampler - modular, streaming capable sampler                     *
 *                                                                         *
 *   Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck   *
 *   Copyright (C) 2005 Christian Schoenebeck                              *
 *   Copyright (C) 2006-2011 Christian Schoenebeck and Andreas Persson     *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU 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 General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU 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                                                   *
 ***************************************************************************/

#ifndef __LS_GIG_FILTER_H__
#define __LS_GIG_FILTER_H__

#include "../../common/global.h"

#include <gig.h>

#include <cmath>

/* TODO: This file contains both generic filters (used by the sfz
   engine) and gig specific filters. It should probably be split up,
   and the generic parts should be moved out of the gig directory. */

/*
 * The formulas for the biquad coefficients come from Robert
 * Bristow-Johnson's Audio EQ Cookbook. The one pole filter formulas
 * come from a post on musicdsp.org. The one poles, biquads and
 * cascaded biquads are modeled after output from Dimension LE and SFZ
 * Player. The gig filters are modeled after output from GigaStudio.
 */
namespace LinuxSampler {

    /**
     * Filter state and parameters for a biquad filter.
     */
    class BiquadFilterData {
    public:
        float b0, b1, b2;
        float a1, a2;

        float x1, x2;
        float y1, y2;
    };

    /**
     * Filter state and parameters for cascaded biquad filters and gig
     * engine filters.
     */
    class FilterData : public BiquadFilterData
    {
    public:
        union {
            // gig filter parameters
            struct {
                float a3;
                float x3;
                float y3;

                float scale;
                float b20;
                float y21, y22, y23;
            };
            // cascaded biquad parameters
            struct {
                BiquadFilterData d2;
                BiquadFilterData d3;
            };
        };
    };

    /**
     * Abstract base class for all filter implementations.
     */
    class FilterBase {
    public:
        virtual float Apply(FilterData& d, float x) const = 0;
        virtual void SetParameters(FilterData& d, float fc, float r,
                                   float fs) const = 0;
        virtual void Reset(FilterData& d) const = 0;
    protected:
        void KillDenormal(float& f) const {
            f += 1e-18f;
            f -= 1e-18f;
        }
    };

    /**
     * One-pole lowpass filter.
     */
    class LowpassFilter1p : public FilterBase {
    public:
        float Apply(FilterData& d, float x) const {
            float y = x + d.a1 * (x - d.y1); // d.b0 * x - d.a1 * d.y1;
            KillDenormal(y);
            d.y1 = y;
            return y;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float c     = 2 - cos(omega);
            d.a1 = -(c - sqrt(c * c - 1));
            // d.b0 = 1 + d.a1;
        }

        void Reset(FilterData& d) const {
            d.y1 = 0;
        }
    };

    /**
     * One pole highpass filter.
     */
    class HighpassFilter1p : public FilterBase {
    public:
        float Apply(FilterData& d, float x) const {
            // d.b0 * x + d.b1 * d.x1 - d.a1 * d.y1;
            float y = d.a1 * (-x + d.x1 - d.y1);
            KillDenormal(y);
            d.x1 = x;
            d.y1 = y;
            return y;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float c     = 2 - cos(omega);
            d.a1 = -(c - sqrt(c * c - 1));
            // d.b0 = -d.a1
            // d.b1 = d.a1
        }

        void Reset(FilterData& d) const {
            d.x1 = 0;
            d.y1 = 0;
        }
    };

    /**
     * Base class for biquad filter implementations.
     */
    class BiquadFilter : public FilterBase {
    protected:
        float ApplyBQ(BiquadFilterData& d, float x) const {
            float y = d.b0 * x + d.b1 * d.x1 + d.b2 * d.x2 + 
                d.a1 * d.y1 + d.a2 * d.y2;
            KillDenormal(y);
            d.x2 = d.x1;
            d.x1 = x;
            d.y2 = d.y1;
            d.y1 = y;
            return y;
        }

    public:
        float Apply(FilterData& d, float x) const {
            return ApplyBQ(d, x);
        }

        void Reset(FilterData& d) const {
            d.x1 = d.x2 = 0;
            d.y1 = d.y2 = 0;
        }
    };

    /**
     * Base class for cascaded double biquad filter (four poles).
     */
    class DoubleBiquadFilter : public BiquadFilter {
    public:
        float Apply(FilterData& d, float x) const {
            return ApplyBQ(d.d2, BiquadFilter::Apply(d, x));
        }

        void Reset(FilterData& d) const {
            BiquadFilter::Reset(d);
            d.d2.x1 = d.d2.x2 = 0;
            d.d2.y1 = d.d2.y2 = 0;
        }
    };

    /**
     * Base class for cascaded triple biquad filter (six poles).
     */
    class TripleBiquadFilter : public DoubleBiquadFilter {
    public:
        float Apply(FilterData& d, float x) const {
            return ApplyBQ(d.d3, DoubleBiquadFilter::Apply(d, x));
        }

        void Reset(FilterData& d) const {
            DoubleBiquadFilter::Reset(d);
            d.d3.x1 = d.d3.x2 = 0;
            d.d3.y1 = d.d3.y2 = 0;
        }
    };


    /** @brief Lowpass Filter
     *
     * Lowpass filter based on biquad filter implementation.
     */
    class LowpassFilter : public BiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
            float a0r = 1.0 / (1.0 + alpha);

            d.b0 = a0r * (1.0 - cs) * 0.5;
            d.b1 = a0r * (1.0 - cs);
            d.b2 = a0r * (1.0 - cs) * 0.5;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Four pole lowpass filter
     *
     * Lowpass filter based on two cascaded biquad filters.
     */
    class LowpassFilter4p : public DoubleBiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2;
            float a0r = 1.0 / (1.0 + alpha);

            d.b0 = a0r * (1.0 - cs) * 0.5;
            d.b1 = a0r * (1.0 - cs);
            d.b2 = a0r * (1.0 - cs) * 0.5;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);

            alpha *= exp(-M_LN10 / 20 * r);
            a0r = 1.0 / (1.0 + alpha);

            d.d2.b0 = a0r * (1.0 - cs) * 0.5;
            d.d2.b1 = a0r * (1.0 - cs);
            d.d2.b2 = a0r * (1.0 - cs) * 0.5;
            d.d2.a1 = a0r * (2.0 * cs);
            d.d2.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Six pole lowpass filter
     *
     * Lowpass filter based on three cascaded biquad filters.
     */
    class LowpassFilter6p : public TripleBiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2;
            float a0r = 1.0 / (1.0 + alpha);

            d.b0 = d.d2.b0 = a0r * (1.0 - cs) * 0.5;
            d.b1 = d.d2.b1 = a0r * (1.0 - cs);
            d.b2 = d.d2.b2 = a0r * (1.0 - cs) * 0.5;
            d.a1 = d.d2.a1 = a0r * (2.0 * cs);
            d.a2 = d.d2.a2 = a0r * (alpha - 1.0);

            alpha *= exp(-M_LN10 / 20 * r);
            a0r = 1.0 / (1.0 + alpha);

            d.d3.b0 = a0r * (1.0 - cs) * 0.5;
            d.d3.b1 = a0r * (1.0 - cs);
            d.d3.b2 = a0r * (1.0 - cs) * 0.5;
            d.d3.a1 = a0r * (2.0 * cs);
            d.d3.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Bandpass filter
     *
     * Bandpass filter based on biquad filter implementation.
     */
    class BandpassFilter : public BiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);

            float a0r = 1.0 / (1.0 + alpha);
            d.b0 = a0r * alpha;
            d.b1 = 0.0;
            d.b2 = a0r * -alpha;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Bandreject filter
     *
     * Bandreject filter based on biquad filter implementation.
     */
    class BandrejectFilter : public BiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);

            float a0r = 1.0 / (1.0 + alpha);
            d.b0 = a0r;
            d.b1 = a0r * (-2.0 * cs);
            d.b2 = a0r;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Highpass filter
     *
     * Highpass filter based on biquad filter implementation.
     */
    class HighpassFilter : public BiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);

            float a0r = 1.0 / (1.0 + alpha);
            d.b0 = a0r * (1.0 + cs) * 0.5;
            d.b1 = a0r * -(1.0 + cs);
            d.b2 = a0r * (1.0 + cs) * 0.5;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Four pole highpass filter
     *
     * Highpass filter based on three cascaded biquad filters.
     */
    class HighpassFilter4p : public DoubleBiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2;

            float a0r = 1.0 / (1.0 + alpha);
            d.b0 = a0r * (1.0 + cs) * 0.5;
            d.b1 = a0r * -(1.0 + cs);
            d.b2 = a0r * (1.0 + cs) * 0.5;
            d.a1 = a0r * (2.0 * cs);
            d.a2 = a0r * (alpha - 1.0);

            alpha *= exp(-M_LN10 / 20 * r);
            a0r = 1.0 / (1.0 + alpha);

            d.d2.b0 = a0r * (1.0 + cs) * 0.5;
            d.d2.b1 = a0r * -(1.0 + cs);
            d.d2.b2 = a0r * (1.0 + cs) * 0.5;
            d.d2.a1 = a0r * (2.0 * cs);
            d.d2.a2 = a0r * (alpha - 1.0);
        }
    };

    /** @brief Six pole highpass filter
     *
     * Highpass filter based on three cascaded biquad filters.
     */
    class HighpassFilter6p : public TripleBiquadFilter {
    public:
        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            float omega = 2.0 * M_PI * fc / fs;
            float sn    = sin(omega);
            float cs    = cos(omega);
            float alpha = sn * M_SQRT1_2;

            float a0r = 1.0 / (1.0 + alpha);
            d.b0 = d.d2.b0 = a0r * (1.0 + cs) * 0.5;
            d.b1 = d.d2.b1 = a0r * -(1.0 + cs);
            d.b2 = d.d2.b2 = a0r * (1.0 + cs) * 0.5;
            d.a1 = d.d2.a1 = a0r * (2.0 * cs);
            d.a2 = d.d2.a2 = a0r * (alpha - 1.0);

            alpha *= exp(-M_LN10 / 20 * r);
            a0r = 1.0 / (1.0 + alpha);

            d.d3.b0 = a0r * (1.0 + cs) * 0.5;
            d.d3.b1 = a0r * -(1.0 + cs);
            d.d3.b2 = a0r * (1.0 + cs) * 0.5;
            d.d3.a1 = a0r * (2.0 * cs);
            d.d3.a2 = a0r * (alpha - 1.0);
        }
    };

namespace gig {

    /**
     * Base class for the gig engine filters.
     */
    class GigFilter : public FilterBase {
    public:
        void Reset(FilterData& d) const {
            d.x1 = d.x2 = d.x3 = 0;
            d.y1 = d.y2 = d.y3 = 0;
        }
    protected:
        float ApplyA(FilterData& d, float x) const {
            float y = x - d.a1 * d.y1 - d.a2 * d.y2 - d.a3 * d.y3;
            KillDenormal(y);
            d.y3 = d.y2;
            d.y2 = d.y1;
            d.y1 = y;
            return y;
        }
    };

#define GIG_PARAM_INIT                                                  \
    float f1 = fc * 0.0075279;                                          \
    float f2 = f1 - 1 + r * fc * (-5.5389e-5 + 1.1982e-7 * fc);         \
    float scale = r < 51 ? 1.0f : 1.3762f - 0.0075073f * r

    class LowpassFilter : public GigFilter {
    public:
        float Apply(FilterData& d, float x) const {
            return ApplyA(d, d.b0 * x);
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            GIG_PARAM_INIT;

            float f1_2 = f1 * f1;
            d.b0 = f1_2 * scale;
            d.a1 = f2;
            d.a2 = f1_2 - 1;
            d.a3 = -f2;
        }
    };

    class BandpassFilter : public GigFilter {
    public:
        float Apply(FilterData& d, float x) const {
            float y = ApplyA(d, d.b0 * x + d.b2 * d.x2);
            d.x2 = d.x1;
            d.x1 = x;
            return y;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            GIG_PARAM_INIT;

            d.b0 = f1 * scale;
            d.b2 = -d.b0;
            d.a1 = f2;
            d.a2 = f1 * f1 - 1;
            d.a3 = -f2;
        }
    };

    class HighpassFilter : public GigFilter {
        float Apply(FilterData& d, float x) const {
            float y = ApplyA(d, -x + d.x1 + d.x2 - d.x3);
            d.x3 = d.x2;
            d.x2 = d.x1;
            d.x1 = x;
            return y * d.scale;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            GIG_PARAM_INIT;

            d.a1 = f2;
            d.a2 = f1 * f1 - 1;
            d.a3 = -f2;
            d.scale = scale;
        }
    };

    class BandrejectFilter : public GigFilter {
        float Apply(FilterData& d, float x) const {
            float y = ApplyA(d, x - d.x1 + d.b2 * d.x2 + d.x3);
            d.x3 = d.x2;
            d.x2 = d.x1;
            d.x1 = x;
            return y * d.scale;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            GIG_PARAM_INIT;

            d.b2 = f1 * f1 - 1;
            d.a1 = f2;
            d.a2 = d.b2;
            d.a3 = -f2;
            d.scale = scale;
        }
    };

    class LowpassTurboFilter : public LowpassFilter {
    public:
        float Apply(FilterData& d, float x) const {
            float y = d.b20 * LowpassFilter::Apply(d, x)
                - d.a1 * d.y21 - d.a2 * d.y22 - d.a3 * d.y23;
            KillDenormal(y);
            d.y23 = d.y22;
            d.y22 = d.y21;
            d.y21 = y;
            return y;
        }

        void SetParameters(FilterData& d, float fc, float r, float fs) const {
            LowpassFilter::SetParameters(d, fc, r, fs);
            d.b20 = d.b0 * 0.5;
        }
    };
} //namespace gig


    /**
     * Main filter class.
     */
    class Filter {
        protected:
            static const LowpassFilter1p  lp1p;
            static const LowpassFilter    lp2p;
            static const LowpassFilter4p  lp4p;
            static const LowpassFilter6p  lp6p;
            static const BandpassFilter   bp2p;
            static const BandrejectFilter br2p;
            static const HighpassFilter1p hp1p;
            static const HighpassFilter   hp2p;
            static const HighpassFilter4p hp4p;
            static const HighpassFilter6p hp6p;
            /**
             * These are filters similar to the ones from Gigasampler.
             */
            static const gig::HighpassFilter     HPFilter;
            static const gig::BandpassFilter     BPFilter;
            static const gig::LowpassFilter      LPFilter;
            static const gig::BandrejectFilter   BRFilter;
            static const gig::LowpassTurboFilter LPTFilter;

            FilterData d;
            const FilterBase* pFilter;

        public:
            Filter() {
                // set filter type to 'lowpass' by default
                pFilter = &LPFilter;
                pFilter->Reset(d);
            }

            enum vcf_type_t {
                vcf_type_gig_lowpass = ::gig::vcf_type_lowpass,
                vcf_type_gig_lowpassturbo = ::gig::vcf_type_lowpassturbo,
                vcf_type_gig_bandpass = ::gig::vcf_type_bandpass,
                vcf_type_gig_highpass = ::gig::vcf_type_highpass,
                vcf_type_gig_bandreject = ::gig::vcf_type_bandreject,
                vcf_type_1p_lowpass,
                vcf_type_1p_highpass,
                vcf_type_2p_lowpass,
                vcf_type_2p_highpass,
                vcf_type_2p_bandpass,
                vcf_type_2p_bandreject,
                vcf_type_4p_lowpass,
                vcf_type_4p_highpass,
                vcf_type_6p_lowpass,
                vcf_type_6p_highpass
            };

            void SetType(vcf_type_t FilterType) {
                switch (FilterType) {
                    case vcf_type_gig_highpass:
                        pFilter = &HPFilter;
                        break;
                    case vcf_type_gig_bandreject:
                        pFilter = &BRFilter;
                        break;
                    case vcf_type_gig_bandpass:
                        pFilter = &BPFilter;
                        break;
                    case vcf_type_gig_lowpassturbo:
                        pFilter = &LPTFilter;
                        break;
                    case vcf_type_1p_lowpass:
                        pFilter = &lp1p;
                        break;
                    case vcf_type_1p_highpass:
                        pFilter = &hp1p;
                        break;
                    case vcf_type_2p_lowpass:
                        pFilter = &lp2p;
                        break;
                    case vcf_type_2p_highpass:
                        pFilter = &hp2p;
                        break;
                    case vcf_type_2p_bandpass:
                        pFilter = &bp2p;
                        break;
                    case vcf_type_2p_bandreject:
                        pFilter = &br2p;
                        break;
                    case vcf_type_4p_lowpass:
                        pFilter = &lp4p;
                        break;
                    case vcf_type_4p_highpass:
                        pFilter = &hp4p;
                        break;
                    case vcf_type_6p_lowpass:
                        pFilter = &lp6p;
                        break;
                    case vcf_type_6p_highpass:
                        pFilter = &hp6p;
                        break;
                    default:
                        pFilter = &LPFilter;
                }
                pFilter->Reset(d);
            }

            void SetParameters(float cutoff, float resonance, float fs) {
                pFilter->SetParameters(d, cutoff, resonance, fs);
            }

            void Reset() {
                return pFilter->Reset(d);
            }

            float Apply(float in) {
                return pFilter->Apply(d, in);
            }
    };

} //namespace LinuxSampler

#endif // __LS_GIG_FILTER_H__
1	/***************************************************************************
2	* *
3	* LinuxSampler - modular, streaming capable sampler *
4	* *
5	* Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck *
6	* Copyright (C) 2005 Christian Schoenebeck *
7	* Copyright (C) 2006-2011 Christian Schoenebeck and Andreas Persson *
8	* *
9	* This program is free software; you can redistribute it and/or modify *
10	* it under the terms of the GNU General Public License as published by *
11	* the Free Software Foundation; either version 2 of the License, or *
12	* (at your option) any later version. *
13	* *
14	* This program is distributed in the hope that it will be useful, *
15	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
16	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17	* GNU General Public License for more details. *
18	* *
19	* You should have received a copy of the GNU General Public License *
20	* along with this program; if not, write to the Free Software *
21	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, *
22	* MA 02111-1307 USA *
23	***************************************************************************/
24
25	#ifndef __LS_GIG_FILTER_H__
26	#define __LS_GIG_FILTER_H__
27
28	#include "../../common/global.h"
29
30	#include <gig.h>
31
32	#include <cmath>
33
34	/* TODO: This file contains both generic filters (used by the sfz
35	engine) and gig specific filters. It should probably be split up,
36	and the generic parts should be moved out of the gig directory. */
37
38	/*
39	* The formulas for the biquad coefficients come from Robert
40	* Bristow-Johnson's Audio EQ Cookbook. The one pole filter formulas
41	* come from a post on musicdsp.org. The one poles, biquads and
42	* cascaded biquads are modeled after output from Dimension LE and SFZ
43	* Player. The gig filters are modeled after output from GigaStudio.
44	*/
45	namespace LinuxSampler {
46
47	/**
48	* Filter state and parameters for a biquad filter.
49	*/
50	class BiquadFilterData {
51	public:
52	float b0, b1, b2;
53	float a1, a2;
54
55	float x1, x2;
56	float y1, y2;
57	};
58
59	/**
60	* Filter state and parameters for cascaded biquad filters and gig
61	* engine filters.
62	*/
63	class FilterData : public BiquadFilterData
64	{
65	public:
66	union {
67	// gig filter parameters
68	struct {
69	float a3;
70	float x3;
71	float y3;
72
73	float scale;
74	float b20;
75	float y21, y22, y23;
76	};
77	// cascaded biquad parameters
78	struct {
79	BiquadFilterData d2;
80	BiquadFilterData d3;
81	};
82	};
83	};
84
85	/**
86	* Abstract base class for all filter implementations.
87	*/
88	class FilterBase {
89	public:
90	virtual float Apply(FilterData& d, float x) const = 0;
91	virtual void SetParameters(FilterData& d, float fc, float r,
92	float fs) const = 0;
93	virtual void Reset(FilterData& d) const = 0;
94	protected:
95	void KillDenormal(float& f) const {
96	f += 1e-18f;
97	f -= 1e-18f;
98	}
99	};
100
101	/**
102	* One-pole lowpass filter.
103	*/
104	class LowpassFilter1p : public FilterBase {
105	public:
106	float Apply(FilterData& d, float x) const {
107	float y = x + d.a1 * (x - d.y1); // d.b0 * x - d.a1 * d.y1;
108	KillDenormal(y);
109	d.y1 = y;
110	return y;
111	}
112
113	void SetParameters(FilterData& d, float fc, float r, float fs) const {
114	float omega = 2.0 * M_PI * fc / fs;
115	float c = 2 - cos(omega);
116	d.a1 = -(c - sqrt(c * c - 1));
117	// d.b0 = 1 + d.a1;
118	}
119
120	void Reset(FilterData& d) const {
121	d.y1 = 0;
122	}
123	};
124
125	/**
126	* One pole highpass filter.
127	*/
128	class HighpassFilter1p : public FilterBase {
129	public:
130	float Apply(FilterData& d, float x) const {
131	// d.b0 * x + d.b1 * d.x1 - d.a1 * d.y1;
132	float y = d.a1 * (-x + d.x1 - d.y1);
133	KillDenormal(y);
134	d.x1 = x;
135	d.y1 = y;
136	return y;
137	}
138
139	void SetParameters(FilterData& d, float fc, float r, float fs) const {
140	float omega = 2.0 * M_PI * fc / fs;
141	float c = 2 - cos(omega);
142	d.a1 = -(c - sqrt(c * c - 1));
143	// d.b0 = -d.a1
144	// d.b1 = d.a1
145	}
146
147	void Reset(FilterData& d) const {
148	d.x1 = 0;
149	d.y1 = 0;
150	}
151	};
152
153	/**
154	* Base class for biquad filter implementations.
155	*/
156	class BiquadFilter : public FilterBase {
157	protected:
158	float ApplyBQ(BiquadFilterData& d, float x) const {
159	float y = d.b0 * x + d.b1 * d.x1 + d.b2 * d.x2 +
160	d.a1 * d.y1 + d.a2 * d.y2;
161	KillDenormal(y);
162	d.x2 = d.x1;
163	d.x1 = x;
164	d.y2 = d.y1;
165	d.y1 = y;
166	return y;
167	}
168
169	public:
170	float Apply(FilterData& d, float x) const {
171	return ApplyBQ(d, x);
172	}
173
174	void Reset(FilterData& d) const {
175	d.x1 = d.x2 = 0;
176	d.y1 = d.y2 = 0;
177	}
178	};
179
180	/**
181	* Base class for cascaded double biquad filter (four poles).
182	*/
183	class DoubleBiquadFilter : public BiquadFilter {
184	public:
185	float Apply(FilterData& d, float x) const {
186	return ApplyBQ(d.d2, BiquadFilter::Apply(d, x));
187	}
188
189	void Reset(FilterData& d) const {
190	BiquadFilter::Reset(d);
191	d.d2.x1 = d.d2.x2 = 0;
192	d.d2.y1 = d.d2.y2 = 0;
193	}
194	};
195
196	/**
197	* Base class for cascaded triple biquad filter (six poles).
198	*/
199	class TripleBiquadFilter : public DoubleBiquadFilter {
200	public:
201	float Apply(FilterData& d, float x) const {
202	return ApplyBQ(d.d3, DoubleBiquadFilter::Apply(d, x));
203	}
204
205	void Reset(FilterData& d) const {
206	DoubleBiquadFilter::Reset(d);
207	d.d3.x1 = d.d3.x2 = 0;
208	d.d3.y1 = d.d3.y2 = 0;
209	}
210	};
211
212
213	/** @brief Lowpass Filter
214	*
215	* Lowpass filter based on biquad filter implementation.
216	*/
217	class LowpassFilter : public BiquadFilter {
218	public:
219	void SetParameters(FilterData& d, float fc, float r, float fs) const {
220	float omega = 2.0 * M_PI * fc / fs;
221	float sn = sin(omega);
222	float cs = cos(omega);
223	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
224	float a0r = 1.0 / (1.0 + alpha);
225
226	d.b0 = a0r * (1.0 - cs) * 0.5;
227	d.b1 = a0r * (1.0 - cs);
228	d.b2 = a0r * (1.0 - cs) * 0.5;
229	d.a1 = a0r * (2.0 * cs);
230	d.a2 = a0r * (alpha - 1.0);
231	}
232	};
233
234	/** @brief Four pole lowpass filter
235	*
236	* Lowpass filter based on two cascaded biquad filters.
237	*/
238	class LowpassFilter4p : public DoubleBiquadFilter {
239	public:
240	void SetParameters(FilterData& d, float fc, float r, float fs) const {
241	float omega = 2.0 * M_PI * fc / fs;
242	float sn = sin(omega);
243	float cs = cos(omega);
244	float alpha = sn * M_SQRT1_2;
245	float a0r = 1.0 / (1.0 + alpha);
246
247	d.b0 = a0r * (1.0 - cs) * 0.5;
248	d.b1 = a0r * (1.0 - cs);
249	d.b2 = a0r * (1.0 - cs) * 0.5;
250	d.a1 = a0r * (2.0 * cs);
251	d.a2 = a0r * (alpha - 1.0);
252
253	alpha = exp(-M_LN10 / 20 r);
254	a0r = 1.0 / (1.0 + alpha);
255
256	d.d2.b0 = a0r * (1.0 - cs) * 0.5;
257	d.d2.b1 = a0r * (1.0 - cs);
258	d.d2.b2 = a0r * (1.0 - cs) * 0.5;
259	d.d2.a1 = a0r * (2.0 * cs);
260	d.d2.a2 = a0r * (alpha - 1.0);
261	}
262	};
263
264	/** @brief Six pole lowpass filter
265	*
266	* Lowpass filter based on three cascaded biquad filters.
267	*/
268	class LowpassFilter6p : public TripleBiquadFilter {
269	public:
270	void SetParameters(FilterData& d, float fc, float r, float fs) const {
271	float omega = 2.0 * M_PI * fc / fs;
272	float sn = sin(omega);
273	float cs = cos(omega);
274	float alpha = sn * M_SQRT1_2;
275	float a0r = 1.0 / (1.0 + alpha);
276
277	d.b0 = d.d2.b0 = a0r * (1.0 - cs) * 0.5;
278	d.b1 = d.d2.b1 = a0r * (1.0 - cs);
279	d.b2 = d.d2.b2 = a0r * (1.0 - cs) * 0.5;
280	d.a1 = d.d2.a1 = a0r * (2.0 * cs);
281	d.a2 = d.d2.a2 = a0r * (alpha - 1.0);
282
283	alpha = exp(-M_LN10 / 20 r);
284	a0r = 1.0 / (1.0 + alpha);
285
286	d.d3.b0 = a0r * (1.0 - cs) * 0.5;
287	d.d3.b1 = a0r * (1.0 - cs);
288	d.d3.b2 = a0r * (1.0 - cs) * 0.5;
289	d.d3.a1 = a0r * (2.0 * cs);
290	d.d3.a2 = a0r * (alpha - 1.0);
291	}
292	};
293
294	/** @brief Bandpass filter
295	*
296	* Bandpass filter based on biquad filter implementation.
297	*/
298	class BandpassFilter : public BiquadFilter {
299	public:
300	void SetParameters(FilterData& d, float fc, float r, float fs) const {
301	float omega = 2.0 * M_PI * fc / fs;
302	float sn = sin(omega);
303	float cs = cos(omega);
304	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
305
306	float a0r = 1.0 / (1.0 + alpha);
307	d.b0 = a0r * alpha;
308	d.b1 = 0.0;
309	d.b2 = a0r * -alpha;
310	d.a1 = a0r * (2.0 * cs);
311	d.a2 = a0r * (alpha - 1.0);
312	}
313	};
314
315	/** @brief Bandreject filter
316	*
317	* Bandreject filter based on biquad filter implementation.
318	*/
319	class BandrejectFilter : public BiquadFilter {
320	public:
321	void SetParameters(FilterData& d, float fc, float r, float fs) const {
322	float omega = 2.0 * M_PI * fc / fs;
323	float sn = sin(omega);
324	float cs = cos(omega);
325	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
326
327	float a0r = 1.0 / (1.0 + alpha);
328	d.b0 = a0r;
329	d.b1 = a0r * (-2.0 * cs);
330	d.b2 = a0r;
331	d.a1 = a0r * (2.0 * cs);
332	d.a2 = a0r * (alpha - 1.0);
333	}
334	};
335
336	/** @brief Highpass filter
337	*
338	* Highpass filter based on biquad filter implementation.
339	*/
340	class HighpassFilter : public BiquadFilter {
341	public:
342	void SetParameters(FilterData& d, float fc, float r, float fs) const {
343	float omega = 2.0 * M_PI * fc / fs;
344	float sn = sin(omega);
345	float cs = cos(omega);
346	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
347
348	float a0r = 1.0 / (1.0 + alpha);
349	d.b0 = a0r * (1.0 + cs) * 0.5;
350	d.b1 = a0r * -(1.0 + cs);
351	d.b2 = a0r * (1.0 + cs) * 0.5;
352	d.a1 = a0r * (2.0 * cs);
353	d.a2 = a0r * (alpha - 1.0);
354	}
355	};
356
357	/** @brief Four pole highpass filter
358	*
359	* Highpass filter based on three cascaded biquad filters.
360	*/
361	class HighpassFilter4p : public DoubleBiquadFilter {
362	public:
363	void SetParameters(FilterData& d, float fc, float r, float fs) const {
364	float omega = 2.0 * M_PI * fc / fs;
365	float sn = sin(omega);
366	float cs = cos(omega);
367	float alpha = sn * M_SQRT1_2;
368
369	float a0r = 1.0 / (1.0 + alpha);
370	d.b0 = a0r * (1.0 + cs) * 0.5;
371	d.b1 = a0r * -(1.0 + cs);
372	d.b2 = a0r * (1.0 + cs) * 0.5;
373	d.a1 = a0r * (2.0 * cs);
374	d.a2 = a0r * (alpha - 1.0);
375
376	alpha = exp(-M_LN10 / 20 r);
377	a0r = 1.0 / (1.0 + alpha);
378
379	d.d2.b0 = a0r * (1.0 + cs) * 0.5;
380	d.d2.b1 = a0r * -(1.0 + cs);
381	d.d2.b2 = a0r * (1.0 + cs) * 0.5;
382	d.d2.a1 = a0r * (2.0 * cs);
383	d.d2.a2 = a0r * (alpha - 1.0);
384	}
385	};
386
387	/** @brief Six pole highpass filter
388	*
389	* Highpass filter based on three cascaded biquad filters.
390	*/
391	class HighpassFilter6p : public TripleBiquadFilter {
392	public:
393	void SetParameters(FilterData& d, float fc, float r, float fs) const {
394	float omega = 2.0 * M_PI * fc / fs;
395	float sn = sin(omega);
396	float cs = cos(omega);
397	float alpha = sn * M_SQRT1_2;
398
399	float a0r = 1.0 / (1.0 + alpha);
400	d.b0 = d.d2.b0 = a0r * (1.0 + cs) * 0.5;
401	d.b1 = d.d2.b1 = a0r * -(1.0 + cs);
402	d.b2 = d.d2.b2 = a0r * (1.0 + cs) * 0.5;
403	d.a1 = d.d2.a1 = a0r * (2.0 * cs);
404	d.a2 = d.d2.a2 = a0r * (alpha - 1.0);
405
406	alpha = exp(-M_LN10 / 20 r);
407	a0r = 1.0 / (1.0 + alpha);
408
409	d.d3.b0 = a0r * (1.0 + cs) * 0.5;
410	d.d3.b1 = a0r * -(1.0 + cs);
411	d.d3.b2 = a0r * (1.0 + cs) * 0.5;
412	d.d3.a1 = a0r * (2.0 * cs);
413	d.d3.a2 = a0r * (alpha - 1.0);
414	}
415	};
416
417	namespace gig {
418
419	/**
420	* Base class for the gig engine filters.
421	*/
422	class GigFilter : public FilterBase {
423	public:
424	void Reset(FilterData& d) const {
425	d.x1 = d.x2 = d.x3 = 0;
426	d.y1 = d.y2 = d.y3 = 0;
427	}
428	protected:
429	float ApplyA(FilterData& d, float x) const {
430	float y = x - d.a1 * d.y1 - d.a2 * d.y2 - d.a3 * d.y3;
431	KillDenormal(y);
432	d.y3 = d.y2;
433	d.y2 = d.y1;
434	d.y1 = y;
435	return y;
436	}
437	};
438
439	#define GIG_PARAM_INIT \
440	float f1 = fc * 0.0075279; \
441	float f2 = f1 - 1 + r * fc * (-5.5389e-5 + 1.1982e-7 * fc); \
442	float scale = r < 51 ? 1.0f : 1.3762f - 0.0075073f * r
443
444	class LowpassFilter : public GigFilter {
445	public:
446	float Apply(FilterData& d, float x) const {
447	return ApplyA(d, d.b0 * x);
448	}
449
450	void SetParameters(FilterData& d, float fc, float r, float fs) const {
451	GIG_PARAM_INIT;
452
453	float f1_2 = f1 * f1;
454	d.b0 = f1_2 * scale;
455	d.a1 = f2;
456	d.a2 = f1_2 - 1;
457	d.a3 = -f2;
458	}
459	};
460
461	class BandpassFilter : public GigFilter {
462	public:
463	float Apply(FilterData& d, float x) const {
464	float y = ApplyA(d, d.b0 * x + d.b2 * d.x2);
465	d.x2 = d.x1;
466	d.x1 = x;
467	return y;
468	}
469
470	void SetParameters(FilterData& d, float fc, float r, float fs) const {
471	GIG_PARAM_INIT;
472
473	d.b0 = f1 * scale;
474	d.b2 = -d.b0;
475	d.a1 = f2;
476	d.a2 = f1 * f1 - 1;
477	d.a3 = -f2;
478	}
479	};
480
481	class HighpassFilter : public GigFilter {
482	float Apply(FilterData& d, float x) const {
483	float y = ApplyA(d, -x + d.x1 + d.x2 - d.x3);
484	d.x3 = d.x2;
485	d.x2 = d.x1;
486	d.x1 = x;
487	return y * d.scale;
488	}
489
490	void SetParameters(FilterData& d, float fc, float r, float fs) const {
491	GIG_PARAM_INIT;
492
493	d.a1 = f2;
494	d.a2 = f1 * f1 - 1;
495	d.a3 = -f2;
496	d.scale = scale;
497	}
498	};
499
500	class BandrejectFilter : public GigFilter {
501	float Apply(FilterData& d, float x) const {
502	float y = ApplyA(d, x - d.x1 + d.b2 * d.x2 + d.x3);
503	d.x3 = d.x2;
504	d.x2 = d.x1;
505	d.x1 = x;
506	return y * d.scale;
507	}
508
509	void SetParameters(FilterData& d, float fc, float r, float fs) const {
510	GIG_PARAM_INIT;
511
512	d.b2 = f1 * f1 - 1;
513	d.a1 = f2;
514	d.a2 = d.b2;
515	d.a3 = -f2;
516	d.scale = scale;
517	}
518	};
519
520	class LowpassTurboFilter : public LowpassFilter {
521	public:
522	float Apply(FilterData& d, float x) const {
523	float y = d.b20 * LowpassFilter::Apply(d, x)
524	- d.a1 * d.y21 - d.a2 * d.y22 - d.a3 * d.y23;
525	KillDenormal(y);
526	d.y23 = d.y22;
527	d.y22 = d.y21;
528	d.y21 = y;
529	return y;
530	}
531
532	void SetParameters(FilterData& d, float fc, float r, float fs) const {
533	LowpassFilter::SetParameters(d, fc, r, fs);
534	d.b20 = d.b0 * 0.5;
535	}
536	};
537	} //namespace gig
538
539
540	/**
541	* Main filter class.
542	*/
543	class Filter {
544	protected:
545	static const LowpassFilter1p lp1p;
546	static const LowpassFilter lp2p;
547	static const LowpassFilter4p lp4p;
548	static const LowpassFilter6p lp6p;
549	static const BandpassFilter bp2p;
550	static const BandrejectFilter br2p;
551	static const HighpassFilter1p hp1p;
552	static const HighpassFilter hp2p;
553	static const HighpassFilter4p hp4p;
554	static const HighpassFilter6p hp6p;
555	/**
556	* These are filters similar to the ones from Gigasampler.
557	*/
558	static const gig::HighpassFilter HPFilter;
559	static const gig::BandpassFilter BPFilter;
560	static const gig::LowpassFilter LPFilter;
561	static const gig::BandrejectFilter BRFilter;
562	static const gig::LowpassTurboFilter LPTFilter;
563
564	FilterData d;
565	const FilterBase* pFilter;
566
567	public:
568	Filter() {
569	// set filter type to 'lowpass' by default
570	pFilter = &LPFilter;
571	pFilter->Reset(d);
572	}
573
574	enum vcf_type_t {
575	vcf_type_gig_lowpass = ::gig::vcf_type_lowpass,
576	vcf_type_gig_lowpassturbo = ::gig::vcf_type_lowpassturbo,
577	vcf_type_gig_bandpass = ::gig::vcf_type_bandpass,
578	vcf_type_gig_highpass = ::gig::vcf_type_highpass,
579	vcf_type_gig_bandreject = ::gig::vcf_type_bandreject,
580	vcf_type_1p_lowpass,
581	vcf_type_1p_highpass,
582	vcf_type_2p_lowpass,
583	vcf_type_2p_highpass,
584	vcf_type_2p_bandpass,
585	vcf_type_2p_bandreject,
586	vcf_type_4p_lowpass,
587	vcf_type_4p_highpass,
588	vcf_type_6p_lowpass,
589	vcf_type_6p_highpass
590	};
591
592	void SetType(vcf_type_t FilterType) {
593	switch (FilterType) {
594	case vcf_type_gig_highpass:
595	pFilter = &HPFilter;
596	break;
597	case vcf_type_gig_bandreject:
598	pFilter = &BRFilter;
599	break;
600	case vcf_type_gig_bandpass:
601	pFilter = &BPFilter;
602	break;
603	case vcf_type_gig_lowpassturbo:
604	pFilter = &LPTFilter;
605	break;
606	case vcf_type_1p_lowpass:
607	pFilter = &lp1p;
608	break;
609	case vcf_type_1p_highpass:
610	pFilter = &hp1p;
611	break;
612	case vcf_type_2p_lowpass:
613	pFilter = &lp2p;
614	break;
615	case vcf_type_2p_highpass:
616	pFilter = &hp2p;
617	break;
618	case vcf_type_2p_bandpass:
619	pFilter = &bp2p;
620	break;
621	case vcf_type_2p_bandreject:
622	pFilter = &br2p;
623	break;
624	case vcf_type_4p_lowpass:
625	pFilter = &lp4p;
626	break;
627	case vcf_type_4p_highpass:
628	pFilter = &hp4p;
629	break;
630	case vcf_type_6p_lowpass:
631	pFilter = &lp6p;
632	break;
633	case vcf_type_6p_highpass:
634	pFilter = &hp6p;
635	break;
636	default:
637	pFilter = &LPFilter;
638	}
639	pFilter->Reset(d);
640	}
641
642	void SetParameters(float cutoff, float resonance, float fs) {
643	pFilter->SetParameters(d, cutoff, resonance, fs);
644	}
645
646	void Reset() {
647	return pFilter->Reset(d);
648	}
649
650	float Apply(float in) {
651	return pFilter->Apply(d, in);
652	}
653	};
654
655	} //namespace LinuxSampler
656
657	#endif // __LS_GIG_FILTER_H__