--- linuxsampler/trunk/src/engines/gig/Filter.h 2004/12/25 21:58:58 328 +++ linuxsampler/trunk/src/engines/gig/Filter.h 2011/04/25 08:12:36 2175 @@ -3,6 +3,8 @@ * 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 * @@ -25,153 +27,631 @@ #include "../../common/global.h" -#if DEBUG_HEADERS -# warning Filter.h included -#endif // DEBUG_HEADERS +#include -#include "../common/BiquadFilter.h" +#include -// TODO: Gigasampler's "Turbo Lowpass" and "Bandreject" filters not implemented yet +/* 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 { -#include "../../lib/fileloader/libgig/gig.h" + /** + * Filter state and parameters for a biquad filter. + */ + class BiquadFilterData { + public: + float b0, b1, b2; + float a1, a2; -#define LSF_BW 0.9 -#define LSF_FB 0.9f + float x1, x2; + float y1, y2; + }; -namespace LinuxSampler { namespace gig { + /** + * 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; + }; + }; + }; /** - * These are filters similar to the ones from Gigasampler. + * 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: - BandpassFilter BasicBPFilter; - HighpassFilter HPFilter; - BandpassFilter BPFilter; - LowpassFilter LPFilter; - BiquadFilter* pFilter; - bq_t scale; - bq_t resonance; - bq_t cutoff; - ::gig::vcf_type_t Type; - static const float fFB = LSF_FB; - public: + 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; - Type = ::gig::vcf_type_lowpass; + pFilter->Reset(d); } - inline bq_t Cutoff() { return cutoff; } - - inline bq_t Resonance() { return resonance; } + 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 + }; - inline void SetType(::gig::vcf_type_t FilterType) { + void SetType(vcf_type_t FilterType) { switch (FilterType) { - case ::gig::vcf_type_highpass: + case vcf_type_gig_highpass: pFilter = &HPFilter; break; - case ::gig::vcf_type_bandreject: //TODO: not implemented yet - Type = ::gig::vcf_type_bandpass; - case ::gig::vcf_type_bandpass: + case vcf_type_gig_bandreject: + pFilter = &BRFilter; + break; + case vcf_type_gig_bandpass: pFilter = &BPFilter; break; - case ::gig::vcf_type_lowpassturbo: //TODO: not implemented yet - default: - Type = ::gig::vcf_type_lowpass; - case ::gig::vcf_type_lowpass: - pFilter = &LPFilter; - - } - Type = FilterType; - } - - inline void SetParameters(bq_t cutoff, bq_t resonance, bq_t fs) { - BasicBPFilter.SetParameters(cutoff, 0.7, fs); - switch (Type) { - case ::gig::vcf_type_highpass: - HPFilter.SetParameters(cutoff, 1.0 - resonance * LSF_BW, fs); + case vcf_type_gig_lowpassturbo: + pFilter = &LPTFilter; break; - case ::gig::vcf_type_bandpass: - BPFilter.SetParameters(cutoff, 1.0 - resonance * LSF_BW, fs); + case vcf_type_1p_lowpass: + pFilter = &lp1p; break; - case ::gig::vcf_type_lowpass: - LPFilter.SetParameters(cutoff, 1.0 - resonance * LSF_BW, fs); + case vcf_type_1p_highpass: + pFilter = &hp1p; break; - } - this->scale = 1.0f - resonance * 0.7f; - this->resonance = resonance; - this->cutoff = cutoff; - } - - inline void SetParameters(biquad_param_t* base, biquad_param_t* main, bq_t cutoff, bq_t resonance, bq_t fs) { - BasicBPFilter.SetParameters(base, cutoff, 0.7, fs); - switch (Type) { - case ::gig::vcf_type_highpass: - HPFilter.SetParameters(main, cutoff, 1.0 - resonance * LSF_BW, fs); + case vcf_type_2p_lowpass: + pFilter = &lp2p; break; - case ::gig::vcf_type_bandpass: - BPFilter.SetParameters(main, cutoff, 1.0 - resonance * LSF_BW, fs); + case vcf_type_2p_highpass: + pFilter = &hp2p; break; - case ::gig::vcf_type_lowpass: - LPFilter.SetParameters(main, cutoff, 1.0 - resonance * LSF_BW, fs); + 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; } - this->scale = 1.0f - resonance * 0.7f; - this->resonance = resonance; - this->cutoff = cutoff; - } - - void Reset() { - BasicBPFilter.Reset(); - HPFilter.Reset(); - BPFilter.Reset(); - LPFilter.Reset(); + pFilter->Reset(d); } - inline bq_t Apply(const bq_t in) { - return pFilter->Apply(in) * this->scale + - BasicBPFilter.ApplyFB(in, this->resonance * LSF_FB) * this->resonance; + void SetParameters(float cutoff, float resonance, float fs) { + pFilter->SetParameters(d, cutoff, resonance, fs); } - inline bq_t Apply(biquad_param_t* base, biquad_param_t* main, const bq_t in) { - return pFilter->Apply(main, in) * this->scale + - BasicBPFilter.ApplyFB(base, in, this->resonance * LSF_FB) * this->resonance; + void Reset() { + return pFilter->Reset(d); } -#if ARCH_X86 - // expects to find input in xmm0 and leaves output in xmm7 - inline void Apply4StepsSSE(biquad_param_t* base, biquad_param_t* main) { - float fb; - __asm__ __volatile__ ( - "movss %0, %%xmm4\n\t" - "mulss %1, %%xmm4 # this->resonance * LSF_FB\n\t" - "movss %%xmm4, %2\n\t" - :: "m" (fFB), /* %0 */ - "m" (resonance), /* %1 */ - "m" (fb) /* %2 */ - ); - BasicBPFilter.ApplyFB4StepsSSE(base, fb); // leaves output in xmm7 - __asm__ __volatile__ ( - "movss %0, %%xmm4\n\t" - "shufps $0, %%xmm4, %%xmm4 # copy to other 3 cells\n\t" - "mulps %%xmm4, %%xmm7 # ApplyFB() * this->resonance\n\t" - :: "m" (resonance) /* %0 */ - ); - pFilter->Apply4StepsSSE(main); // leaves output in xmm6 - __asm__ __volatile__ ( - "movss %0, %%xmm5\n\t" - "shufps $0, %%xmm5, %%xmm5 # copy to other 3 cells\n\t" - "mulps %%xmm5, %%xmm6 # Apply() * this->scale\n\t" - "addps %%xmm6, %%xmm7 # xmm7 = result\n\t" - :: "m" (scale) /* %0 */ - ); + float Apply(float in) { + return pFilter->Apply(d, in); } -#endif // ARCH_X86 - }; -}} //namespace LinuxSampler::gig +} //namespace LinuxSampler #endif // __LS_GIG_FILTER_H__