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-2017 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"

#if AC_APPLE_UNIVERSAL_BUILD
# include <libgig/gig.h>
#else
# include <gig.h>
#endif

#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 {
            if (fc < 1.f) fc = 1.f; // this lowpass turbo filter cannot cope with cutoff being zero (would cause click sounds)
            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 {
                // GigaStudio original filter types
                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,
                // own filter types
                vcf_type_1p_lowpass = ::gig::vcf_type_lowpass_1p,
                vcf_type_1p_highpass = ::gig::vcf_type_highpass_1p,
                vcf_type_2p_lowpass = ::gig::vcf_type_lowpass_2p,
                vcf_type_2p_highpass = ::gig::vcf_type_highpass_2p,
                vcf_type_2p_bandpass = ::gig::vcf_type_bandpass_2p,
                vcf_type_2p_bandreject = ::gig::vcf_type_bandreject_2p,
                vcf_type_4p_lowpass = ::gig::vcf_type_lowpass_4p,
                vcf_type_4p_highpass = ::gig::vcf_type_highpass_4p,
                vcf_type_6p_lowpass = ::gig::vcf_type_lowpass_6p,
                vcf_type_6p_highpass = ::gig::vcf_type_highpass_6p,
            };

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