engines/common/BiquadFilter.h

/***************************************************************************
 *                                                                         *
 *   LinuxSampler - modular, streaming capable sampler                     *
 *                                                                         *
 *   Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck   *
 *                                                                         *
 *   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_BIQUADFILTER_H__
#define __LS_BIQUADFILTER_H__

/// ln(2) / 2
#define LN_2_2                  0.34657359f

#ifndef LIMIT
# define LIMIT(v,l,u)           (v < l ? l : (v > u ? u : v))
#endif

namespace LinuxSampler {

    typedef float bq_t;

    /**
    * Bi-quadratic filter
    * (adapted from lisp code by Eli Brandt, http://www.cs.cmu.edu/~eli/)
    */
    class BiquadFilter {
        protected:
            bq_t a1;
            bq_t a2;
            bq_t b0;
            bq_t b1;
            bq_t b2;
            bq_t x1;
            bq_t x2;
            bq_t y1;
            bq_t y2;

            /**
            * Prevent \a f from going into denormal mode which would slow down
            * subsequent floating point calculations, we achieve that by setting
            * \a f to zero when it falls under the denormal threshold value.
            */
            inline void KillDenormal(float& f) {
                // TODO: this is a generic solution for 32bit floats, should be replaced by CPU specific asm code
                f += 1e-18f;
                f -= 1e-18f;
            }
        public:
            inline BiquadFilter() {
                x1 = 0.0f;
                x2 = 0.0f;
                y1 = 0.0f;
                y2 = 0.0f;
            }

            inline bq_t Apply(const bq_t x) {
                bq_t y;

                y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
                    this->a1 * this->y1 + this->a2 * this->y2;
                KillDenormal(y);
                this->x2 = this->x1;
                this->x1 = x;
                this->y2 = this->y1;
                this->y1 = y;

                return y;
            }

            inline bq_t ApplyFB(bq_t x, const bq_t fb) {
                bq_t y;

                x += this->y1 * fb * 0.98;
                y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
                    this->a1 * this->y1 + this->a2 * this->y2;
                KillDenormal(y);
                this->x2 = this->x1;
                this->x1 = x;
                this->y2 = this->y1;
                this->y1 = y;

                return y;
            }
    };

    class LowpassFilter : public BiquadFilter {
        public:
            inline LowpassFilter() : BiquadFilter() {}

            inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
                bq_t omega = 2.0 * M_PI * fc / fs;
                bq_t sn    = sin(omega);
                bq_t cs    = cos(omega);
                bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);

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

    class BandpassFilter : public BiquadFilter {
        public:
            inline BandpassFilter() : BiquadFilter() {}

            inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
                bq_t omega = 2.0 * M_PI * fc / fs;
                bq_t sn    = sin(omega);
                bq_t cs    = cos(omega);
                bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);

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

    class HighpassFilter : public BiquadFilter {
        public:
            inline HighpassFilter() : BiquadFilter() {}

            inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
                bq_t omega = 2.0 * M_PI * fc / fs;
                bq_t sn    = sin(omega);
                bq_t cs    = cos(omega);
                bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);

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

} // namespace LinuxSampler

#endif // __LS_BIQUADFILTER_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	53	* *
7			* This program is free software; you can redistribute it and/or modify *
8			* it under the terms of the GNU General Public License as published by *
9			* the Free Software Foundation; either version 2 of the License, or *
10			* (at your option) any later version. *
11			* *
12			* This program is distributed in the hope that it will be useful, *
13			* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15			* GNU General Public License for more details. *
16			* *
17			* You should have received a copy of the GNU General Public License *
18			* along with this program; if not, write to the Free Software *
19			* Foundation, Inc., 59 Temple Place, Suite 330, Boston, *
20			* MA 02111-1307 USA *
21			***************************************************************************/
22
23			#ifndef __LS_BIQUADFILTER_H__
24			#define __LS_BIQUADFILTER_H__
25
26			/// ln(2) / 2
27			#define LN_2_2 0.34657359f
28
29			#ifndef LIMIT
30			# define LIMIT(v,l,u) (v < l ? l : (v > u ? u : v))
31			#endif
32
33			namespace LinuxSampler {
34
35			typedef float bq_t;
36
37			/**
38			* Bi-quadratic filter
39			* (adapted from lisp code by Eli Brandt, http://www.cs.cmu.edu/~eli/)
40			*/
41			class BiquadFilter {
42			protected:
43			bq_t a1;
44			bq_t a2;
45			bq_t b0;
46			bq_t b1;
47			bq_t b2;
48			bq_t x1;
49			bq_t x2;
50			bq_t y1;
51			bq_t y2;
52
53			/**
54			* Prevent \a f from going into denormal mode which would slow down
55			* subsequent floating point calculations, we achieve that by setting
56			* \a f to zero when it falls under the denormal threshold value.
57			*/
58			inline void KillDenormal(float& f) {
59			// TODO: this is a generic solution for 32bit floats, should be replaced by CPU specific asm code
60			f += 1e-18f;
61			f -= 1e-18f;
62			}
63			public:
64			inline BiquadFilter() {
65			x1 = 0.0f;
66			x2 = 0.0f;
67			y1 = 0.0f;
68			y2 = 0.0f;
69			}
70
71			inline bq_t Apply(const bq_t x) {
72			bq_t y;
73
74			y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
75			this->a1 * this->y1 + this->a2 * this->y2;
76			KillDenormal(y);
77			this->x2 = this->x1;
78			this->x1 = x;
79			this->y2 = this->y1;
80			this->y1 = y;
81
82			return y;
83			}
84
85			inline bq_t ApplyFB(bq_t x, const bq_t fb) {
86			bq_t y;
87
88			x += this->y1 * fb * 0.98;
89			y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
90			this->a1 * this->y1 + this->a2 * this->y2;
91			KillDenormal(y);
92			this->x2 = this->x1;
93			this->x1 = x;
94			this->y2 = this->y1;
95			this->y1 = y;
96
97			return y;
98			}
99			};
100
101			class LowpassFilter : public BiquadFilter {
102			public:
103			inline LowpassFilter() : BiquadFilter() {}
104
105			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
106			bq_t omega = 2.0 * M_PI * fc / fs;
107			bq_t sn = sin(omega);
108			bq_t cs = cos(omega);
109			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
110
111			const float a0r = 1.0 / (1.0 + alpha);
112			this->b0 = a0r * (1.0 - cs) * 0.5;
113			this->b1 = a0r * (1.0 - cs);
114			this->b2 = a0r * (1.0 - cs) * 0.5;
115			this->a1 = a0r * (2.0 * cs);
116			this->a2 = a0r * (alpha - 1.0);
117			}
118			};
119
120			class BandpassFilter : public BiquadFilter {
121			public:
122			inline BandpassFilter() : BiquadFilter() {}
123
124			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
125			bq_t omega = 2.0 * M_PI * fc / fs;
126			bq_t sn = sin(omega);
127			bq_t cs = cos(omega);
128			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
129
130			const float a0r = 1.0 / (1.0 + alpha);
131			this->b0 = a0r * alpha;
132			this->b1 = 0.0;
133			this->b2 = a0r * -alpha;
134			this->a1 = a0r * (2.0 * cs);
135			this->a2 = a0r * (alpha - 1.0);
136			}
137			};
138
139			class HighpassFilter : public BiquadFilter {
140			public:
141			inline HighpassFilter() : BiquadFilter() {}
142
143			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
144			bq_t omega = 2.0 * M_PI * fc / fs;
145			bq_t sn = sin(omega);
146			bq_t cs = cos(omega);
147			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
148
149			const float a0r = 1.0 / (1.0 + alpha);
150			this->b0 = a0r * (1.0 + cs) * 0.5;
151			this->b1 = a0r * -(1.0 + cs);
152			this->b2 = a0r * (1.0 + cs) * 0.5;
153			this->a1 = a0r * (2.0 * cs);
154			this->a2 = a0r * (alpha - 1.0);
155			}
156			};
157
158			} // namespace LinuxSampler
159
160			#endif // __LS_BIQUADFILTER_H__