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:
        LowpassFilter1p() { }

        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:
        HighpassFilter1p() { }

        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:
        LowpassFilter() { }

        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:
        LowpassFilter4p() { }

        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:
        LowpassFilter6p() { }

        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:
        BandpassFilter() { }

        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:
        BandrejectFilter() { }

        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:
        HighpassFilter() { }

        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:
        HighpassFilter4p() { }

        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:
        HighpassFilter6p() { }

        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:
        LowpassFilter() { }

        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:
        BandpassFilter() { }

        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 {
    public:
        HighpassFilter() { }

        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 {
    public:
        BandrejectFilter() { }

        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:
        LowpassTurboFilter() { }

        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	LowpassFilter1p() { }
107
108	float Apply(FilterData& d, float x) const {
109	float y = x + d.a1 * (x - d.y1); // d.b0 * x - d.a1 * d.y1;
110	KillDenormal(y);
111	d.y1 = y;
112	return y;
113	}
114
115	void SetParameters(FilterData& d, float fc, float r, float fs) const {
116	float omega = 2.0 * M_PI * fc / fs;
117	float c = 2 - cos(omega);
118	d.a1 = -(c - sqrt(c * c - 1));
119	// d.b0 = 1 + d.a1;
120	}
121
122	void Reset(FilterData& d) const {
123	d.y1 = 0;
124	}
125	};
126
127	/**
128	* One pole highpass filter.
129	*/
130	class HighpassFilter1p : public FilterBase {
131	public:
132	HighpassFilter1p() { }
133
134	float Apply(FilterData& d, float x) const {
135	// d.b0 * x + d.b1 * d.x1 - d.a1 * d.y1;
136	float y = d.a1 * (-x + d.x1 - d.y1);
137	KillDenormal(y);
138	d.x1 = x;
139	d.y1 = y;
140	return y;
141	}
142
143	void SetParameters(FilterData& d, float fc, float r, float fs) const {
144	float omega = 2.0 * M_PI * fc / fs;
145	float c = 2 - cos(omega);
146	d.a1 = -(c - sqrt(c * c - 1));
147	// d.b0 = -d.a1
148	// d.b1 = d.a1
149	}
150
151	void Reset(FilterData& d) const {
152	d.x1 = 0;
153	d.y1 = 0;
154	}
155	};
156
157	/**
158	* Base class for biquad filter implementations.
159	*/
160	class BiquadFilter : public FilterBase {
161	protected:
162	float ApplyBQ(BiquadFilterData& d, float x) const {
163	float y = d.b0 * x + d.b1 * d.x1 + d.b2 * d.x2 +
164	d.a1 * d.y1 + d.a2 * d.y2;
165	KillDenormal(y);
166	d.x2 = d.x1;
167	d.x1 = x;
168	d.y2 = d.y1;
169	d.y1 = y;
170	return y;
171	}
172
173	public:
174	float Apply(FilterData& d, float x) const {
175	return ApplyBQ(d, x);
176	}
177
178	void Reset(FilterData& d) const {
179	d.x1 = d.x2 = 0;
180	d.y1 = d.y2 = 0;
181	}
182	};
183
184	/**
185	* Base class for cascaded double biquad filter (four poles).
186	*/
187	class DoubleBiquadFilter : public BiquadFilter {
188	public:
189	float Apply(FilterData& d, float x) const {
190	return ApplyBQ(d.d2, BiquadFilter::Apply(d, x));
191	}
192
193	void Reset(FilterData& d) const {
194	BiquadFilter::Reset(d);
195	d.d2.x1 = d.d2.x2 = 0;
196	d.d2.y1 = d.d2.y2 = 0;
197	}
198	};
199
200	/**
201	* Base class for cascaded triple biquad filter (six poles).
202	*/
203	class TripleBiquadFilter : public DoubleBiquadFilter {
204	public:
205	float Apply(FilterData& d, float x) const {
206	return ApplyBQ(d.d3, DoubleBiquadFilter::Apply(d, x));
207	}
208
209	void Reset(FilterData& d) const {
210	DoubleBiquadFilter::Reset(d);
211	d.d3.x1 = d.d3.x2 = 0;
212	d.d3.y1 = d.d3.y2 = 0;
213	}
214	};
215
216
217	/** @brief Lowpass Filter
218	*
219	* Lowpass filter based on biquad filter implementation.
220	*/
221	class LowpassFilter : public BiquadFilter {
222	public:
223	LowpassFilter() { }
224
225	void SetParameters(FilterData& d, float fc, float r, float fs) const {
226	float omega = 2.0 * M_PI * fc / fs;
227	float sn = sin(omega);
228	float cs = cos(omega);
229	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
230	float a0r = 1.0 / (1.0 + alpha);
231
232	d.b0 = a0r * (1.0 - cs) * 0.5;
233	d.b1 = a0r * (1.0 - cs);
234	d.b2 = a0r * (1.0 - cs) * 0.5;
235	d.a1 = a0r * (2.0 * cs);
236	d.a2 = a0r * (alpha - 1.0);
237	}
238	};
239
240	/** @brief Four pole lowpass filter
241	*
242	* Lowpass filter based on two cascaded biquad filters.
243	*/
244	class LowpassFilter4p : public DoubleBiquadFilter {
245	public:
246	LowpassFilter4p() { }
247
248	void SetParameters(FilterData& d, float fc, float r, float fs) const {
249	float omega = 2.0 * M_PI * fc / fs;
250	float sn = sin(omega);
251	float cs = cos(omega);
252	float alpha = sn * M_SQRT1_2;
253	float a0r = 1.0 / (1.0 + alpha);
254
255	d.b0 = a0r * (1.0 - cs) * 0.5;
256	d.b1 = a0r * (1.0 - cs);
257	d.b2 = a0r * (1.0 - cs) * 0.5;
258	d.a1 = a0r * (2.0 * cs);
259	d.a2 = a0r * (alpha - 1.0);
260
261	alpha = exp(-M_LN10 / 20 r);
262	a0r = 1.0 / (1.0 + alpha);
263
264	d.d2.b0 = a0r * (1.0 - cs) * 0.5;
265	d.d2.b1 = a0r * (1.0 - cs);
266	d.d2.b2 = a0r * (1.0 - cs) * 0.5;
267	d.d2.a1 = a0r * (2.0 * cs);
268	d.d2.a2 = a0r * (alpha - 1.0);
269	}
270	};
271
272	/** @brief Six pole lowpass filter
273	*
274	* Lowpass filter based on three cascaded biquad filters.
275	*/
276	class LowpassFilter6p : public TripleBiquadFilter {
277	public:
278	LowpassFilter6p() { }
279
280	void SetParameters(FilterData& d, float fc, float r, float fs) const {
281	float omega = 2.0 * M_PI * fc / fs;
282	float sn = sin(omega);
283	float cs = cos(omega);
284	float alpha = sn * M_SQRT1_2;
285	float a0r = 1.0 / (1.0 + alpha);
286
287	d.b0 = d.d2.b0 = a0r * (1.0 - cs) * 0.5;
288	d.b1 = d.d2.b1 = a0r * (1.0 - cs);
289	d.b2 = d.d2.b2 = a0r * (1.0 - cs) * 0.5;
290	d.a1 = d.d2.a1 = a0r * (2.0 * cs);
291	d.a2 = d.d2.a2 = a0r * (alpha - 1.0);
292
293	alpha = exp(-M_LN10 / 20 r);
294	a0r = 1.0 / (1.0 + alpha);
295
296	d.d3.b0 = a0r * (1.0 - cs) * 0.5;
297	d.d3.b1 = a0r * (1.0 - cs);
298	d.d3.b2 = a0r * (1.0 - cs) * 0.5;
299	d.d3.a1 = a0r * (2.0 * cs);
300	d.d3.a2 = a0r * (alpha - 1.0);
301	}
302	};
303
304	/** @brief Bandpass filter
305	*
306	* Bandpass filter based on biquad filter implementation.
307	*/
308	class BandpassFilter : public BiquadFilter {
309	public:
310	BandpassFilter() { }
311
312	void SetParameters(FilterData& d, float fc, float r, float fs) const {
313	float omega = 2.0 * M_PI * fc / fs;
314	float sn = sin(omega);
315	float cs = cos(omega);
316	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
317
318	float a0r = 1.0 / (1.0 + alpha);
319	d.b0 = a0r * alpha;
320	d.b1 = 0.0;
321	d.b2 = a0r * -alpha;
322	d.a1 = a0r * (2.0 * cs);
323	d.a2 = a0r * (alpha - 1.0);
324	}
325	};
326
327	/** @brief Bandreject filter
328	*
329	* Bandreject filter based on biquad filter implementation.
330	*/
331	class BandrejectFilter : public BiquadFilter {
332	public:
333	BandrejectFilter() { }
334
335	void SetParameters(FilterData& d, float fc, float r, float fs) const {
336	float omega = 2.0 * M_PI * fc / fs;
337	float sn = sin(omega);
338	float cs = cos(omega);
339	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
340
341	float a0r = 1.0 / (1.0 + alpha);
342	d.b0 = a0r;
343	d.b1 = a0r * (-2.0 * cs);
344	d.b2 = a0r;
345	d.a1 = a0r * (2.0 * cs);
346	d.a2 = a0r * (alpha - 1.0);
347	}
348	};
349
350	/** @brief Highpass filter
351	*
352	* Highpass filter based on biquad filter implementation.
353	*/
354	class HighpassFilter : public BiquadFilter {
355	public:
356	HighpassFilter() { }
357
358	void SetParameters(FilterData& d, float fc, float r, float fs) const {
359	float omega = 2.0 * M_PI * fc / fs;
360	float sn = sin(omega);
361	float cs = cos(omega);
362	float alpha = sn * M_SQRT1_2 * exp(-M_LN10 / 20 * r);
363
364	float a0r = 1.0 / (1.0 + alpha);
365	d.b0 = a0r * (1.0 + cs) * 0.5;
366	d.b1 = a0r * -(1.0 + cs);
367	d.b2 = a0r * (1.0 + cs) * 0.5;
368	d.a1 = a0r * (2.0 * cs);
369	d.a2 = a0r * (alpha - 1.0);
370	}
371	};
372
373	/** @brief Four pole highpass filter
374	*
375	* Highpass filter based on three cascaded biquad filters.
376	*/
377	class HighpassFilter4p : public DoubleBiquadFilter {
378	public:
379	HighpassFilter4p() { }
380
381	void SetParameters(FilterData& d, float fc, float r, float fs) const {
382	float omega = 2.0 * M_PI * fc / fs;
383	float sn = sin(omega);
384	float cs = cos(omega);
385	float alpha = sn * M_SQRT1_2;
386
387	float a0r = 1.0 / (1.0 + alpha);
388	d.b0 = a0r * (1.0 + cs) * 0.5;
389	d.b1 = a0r * -(1.0 + cs);
390	d.b2 = a0r * (1.0 + cs) * 0.5;
391	d.a1 = a0r * (2.0 * cs);
392	d.a2 = a0r * (alpha - 1.0);
393
394	alpha = exp(-M_LN10 / 20 r);
395	a0r = 1.0 / (1.0 + alpha);
396
397	d.d2.b0 = a0r * (1.0 + cs) * 0.5;
398	d.d2.b1 = a0r * -(1.0 + cs);
399	d.d2.b2 = a0r * (1.0 + cs) * 0.5;
400	d.d2.a1 = a0r * (2.0 * cs);
401	d.d2.a2 = a0r * (alpha - 1.0);
402	}
403	};
404
405	/** @brief Six pole highpass filter
406	*
407	* Highpass filter based on three cascaded biquad filters.
408	*/
409	class HighpassFilter6p : public TripleBiquadFilter {
410	public:
411	HighpassFilter6p() { }
412
413	void SetParameters(FilterData& d, float fc, float r, float fs) const {
414	float omega = 2.0 * M_PI * fc / fs;
415	float sn = sin(omega);
416	float cs = cos(omega);
417	float alpha = sn * M_SQRT1_2;
418
419	float a0r = 1.0 / (1.0 + alpha);
420	d.b0 = d.d2.b0 = a0r * (1.0 + cs) * 0.5;
421	d.b1 = d.d2.b1 = a0r * -(1.0 + cs);
422	d.b2 = d.d2.b2 = a0r * (1.0 + cs) * 0.5;
423	d.a1 = d.d2.a1 = a0r * (2.0 * cs);
424	d.a2 = d.d2.a2 = a0r * (alpha - 1.0);
425
426	alpha = exp(-M_LN10 / 20 r);
427	a0r = 1.0 / (1.0 + alpha);
428
429	d.d3.b0 = a0r * (1.0 + cs) * 0.5;
430	d.d3.b1 = a0r * -(1.0 + cs);
431	d.d3.b2 = a0r * (1.0 + cs) * 0.5;
432	d.d3.a1 = a0r * (2.0 * cs);
433	d.d3.a2 = a0r * (alpha - 1.0);
434	}
435	};
436
437	namespace gig {
438
439	/**
440	* Base class for the gig engine filters.
441	*/
442	class GigFilter : public FilterBase {
443	public:
444	void Reset(FilterData& d) const {
445	d.x1 = d.x2 = d.x3 = 0;
446	d.y1 = d.y2 = d.y3 = 0;
447	}
448	protected:
449	float ApplyA(FilterData& d, float x) const {
450	float y = x - d.a1 * d.y1 - d.a2 * d.y2 - d.a3 * d.y3;
451	KillDenormal(y);
452	d.y3 = d.y2;
453	d.y2 = d.y1;
454	d.y1 = y;
455	return y;
456	}
457	};
458
459	#define GIG_PARAM_INIT \
460	float f1 = fc * 0.0075279; \
461	float f2 = f1 - 1 + r * fc * (-5.5389e-5 + 1.1982e-7 * fc); \
462	float scale = r < 51 ? 1.0f : 1.3762f - 0.0075073f * r
463
464	class LowpassFilter : public GigFilter {
465	public:
466	LowpassFilter() { }
467
468	float Apply(FilterData& d, float x) const {
469	return ApplyA(d, d.b0 * x);
470	}
471
472	void SetParameters(FilterData& d, float fc, float r, float fs) const {
473	GIG_PARAM_INIT;
474
475	float f1_2 = f1 * f1;
476	d.b0 = f1_2 * scale;
477	d.a1 = f2;
478	d.a2 = f1_2 - 1;
479	d.a3 = -f2;
480	}
481	};
482
483	class BandpassFilter : public GigFilter {
484	public:
485	BandpassFilter() { }
486
487	float Apply(FilterData& d, float x) const {
488	float y = ApplyA(d, d.b0 * x + d.b2 * d.x2);
489	d.x2 = d.x1;
490	d.x1 = x;
491	return y;
492	}
493
494	void SetParameters(FilterData& d, float fc, float r, float fs) const {
495	GIG_PARAM_INIT;
496
497	d.b0 = f1 * scale;
498	d.b2 = -d.b0;
499	d.a1 = f2;
500	d.a2 = f1 * f1 - 1;
501	d.a3 = -f2;
502	}
503	};
504
505	class HighpassFilter : public GigFilter {
506	public:
507	HighpassFilter() { }
508
509	float Apply(FilterData& d, float x) const {
510	float y = ApplyA(d, -x + d.x1 + d.x2 - d.x3);
511	d.x3 = d.x2;
512	d.x2 = d.x1;
513	d.x1 = x;
514	return y * d.scale;
515	}
516
517	void SetParameters(FilterData& d, float fc, float r, float fs) const {
518	GIG_PARAM_INIT;
519
520	d.a1 = f2;
521	d.a2 = f1 * f1 - 1;
522	d.a3 = -f2;
523	d.scale = scale;
524	}
525	};
526
527	class BandrejectFilter : public GigFilter {
528	public:
529	BandrejectFilter() { }
530
531	float Apply(FilterData& d, float x) const {
532	float y = ApplyA(d, x - d.x1 + d.b2 * d.x2 + d.x3);
533	d.x3 = d.x2;
534	d.x2 = d.x1;
535	d.x1 = x;
536	return y * d.scale;
537	}
538
539	void SetParameters(FilterData& d, float fc, float r, float fs) const {
540	GIG_PARAM_INIT;
541
542	d.b2 = f1 * f1 - 1;
543	d.a1 = f2;
544	d.a2 = d.b2;
545	d.a3 = -f2;
546	d.scale = scale;
547	}
548	};
549
550	class LowpassTurboFilter : public LowpassFilter {
551	public:
552	LowpassTurboFilter() { }
553
554	float Apply(FilterData& d, float x) const {
555	float y = d.b20 * LowpassFilter::Apply(d, x)
556	- d.a1 * d.y21 - d.a2 * d.y22 - d.a3 * d.y23;
557	KillDenormal(y);
558	d.y23 = d.y22;
559	d.y22 = d.y21;
560	d.y21 = y;
561	return y;
562	}
563
564	void SetParameters(FilterData& d, float fc, float r, float fs) const {
565	LowpassFilter::SetParameters(d, fc, r, fs);
566	d.b20 = d.b0 * 0.5;
567	}
568	};
569	} //namespace gig
570
571
572	/**
573	* Main filter class.
574	*/
575	class Filter {
576	protected:
577	static const LowpassFilter1p lp1p;
578	static const LowpassFilter lp2p;
579	static const LowpassFilter4p lp4p;
580	static const LowpassFilter6p lp6p;
581	static const BandpassFilter bp2p;
582	static const BandrejectFilter br2p;
583	static const HighpassFilter1p hp1p;
584	static const HighpassFilter hp2p;
585	static const HighpassFilter4p hp4p;
586	static const HighpassFilter6p hp6p;
587	/**
588	* These are filters similar to the ones from Gigasampler.
589	*/
590	static const gig::HighpassFilter HPFilter;
591	static const gig::BandpassFilter BPFilter;
592	static const gig::LowpassFilter LPFilter;
593	static const gig::BandrejectFilter BRFilter;
594	static const gig::LowpassTurboFilter LPTFilter;
595
596	FilterData d;
597	const FilterBase* pFilter;
598
599	public:
600	Filter() {
601	// set filter type to 'lowpass' by default
602	pFilter = &LPFilter;
603	pFilter->Reset(d);
604	}
605
606	enum vcf_type_t {
607	vcf_type_gig_lowpass = ::gig::vcf_type_lowpass,
608	vcf_type_gig_lowpassturbo = ::gig::vcf_type_lowpassturbo,
609	vcf_type_gig_bandpass = ::gig::vcf_type_bandpass,
610	vcf_type_gig_highpass = ::gig::vcf_type_highpass,
611	vcf_type_gig_bandreject = ::gig::vcf_type_bandreject,
612	vcf_type_1p_lowpass,
613	vcf_type_1p_highpass,
614	vcf_type_2p_lowpass,
615	vcf_type_2p_highpass,
616	vcf_type_2p_bandpass,
617	vcf_type_2p_bandreject,
618	vcf_type_4p_lowpass,
619	vcf_type_4p_highpass,
620	vcf_type_6p_lowpass,
621	vcf_type_6p_highpass
622	};
623
624	void SetType(vcf_type_t FilterType) {
625	switch (FilterType) {
626	case vcf_type_gig_highpass:
627	pFilter = &HPFilter;
628	break;
629	case vcf_type_gig_bandreject:
630	pFilter = &BRFilter;
631	break;
632	case vcf_type_gig_bandpass:
633	pFilter = &BPFilter;
634	break;
635	case vcf_type_gig_lowpassturbo:
636	pFilter = &LPTFilter;
637	break;
638	case vcf_type_1p_lowpass:
639	pFilter = &lp1p;
640	break;
641	case vcf_type_1p_highpass:
642	pFilter = &hp1p;
643	break;
644	case vcf_type_2p_lowpass:
645	pFilter = &lp2p;
646	break;
647	case vcf_type_2p_highpass:
648	pFilter = &hp2p;
649	break;
650	case vcf_type_2p_bandpass:
651	pFilter = &bp2p;
652	break;
653	case vcf_type_2p_bandreject:
654	pFilter = &br2p;
655	break;
656	case vcf_type_4p_lowpass:
657	pFilter = &lp4p;
658	break;
659	case vcf_type_4p_highpass:
660	pFilter = &hp4p;
661	break;
662	case vcf_type_6p_lowpass:
663	pFilter = &lp6p;
664	break;
665	case vcf_type_6p_highpass:
666	pFilter = &hp6p;
667	break;
668	default:
669	pFilter = &LPFilter;
670	}
671	pFilter->Reset(d);
672	}
673
674	void SetParameters(float cutoff, float resonance, float fs) {
675	pFilter->SetParameters(d, cutoff, resonance, fs);
676	}
677
678	void Reset() {
679	return pFilter->Reset(d);
680	}
681
682	float Apply(float in) {
683	return pFilter->Apply(d, in);
684	}
685	};
686
687	} //namespace LinuxSampler
688
689	#endif // __LS_GIG_FILTER_H__