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;

    /**
     * Internal parameters of the biquad filter, which are actually the
     * final parameters of the filter's transfer function. This strucure is
     * only needed when these parameters should stored outside the
     * BiquadFilter class, e.g. to save calculation time by sharing them
     * between multiple filters.
     */
    struct biquad_param_t {
        bq_t a1;
        bq_t a2;
        bq_t b0;
        bq_t b1;
        bq_t b2;
    };

    /**
     * Bi-quadratic filter
     * (adapted from lisp code by Eli Brandt, http://www.cs.cmu.edu/~eli/)
     */
    class BiquadFilter {
        protected:
            // following five variables are only used if no external biquad_param_t reference is used
            bq_t a1;
            bq_t a2;
            bq_t b0;
            bq_t b1;
            bq_t b2;
            // following four variables are used to buffer the feedback
            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(bq_t& 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 Apply(biquad_param_t* param, const bq_t x) {
                bq_t y;

                y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
                    param->a1 * this->y1 + param->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;
            }

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

                x += this->y1 * fb * 0.98;
                y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
                    param->a1 * this->y1 + param->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);
            }

            inline void SetParameters(biquad_param_t* param, 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);
                param->b0 = a0r * (1.0 - cs) * 0.5;
                param->b1 = a0r * (1.0 - cs);
                param->b2 = a0r * (1.0 - cs) * 0.5;
                param->a1 = a0r * (2.0 * cs);
                param->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);
            }

            inline void SetParameters(biquad_param_t* param, 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);
                param->b0 = a0r * alpha;
                param->b1 = 0.0;
                param->b2 = a0r * -alpha;
                param->a1 = a0r * (2.0 * cs);
                param->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);
            }

            inline void SetParameters(biquad_param_t* param, 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);
                param->b0 = a0r * (1.0 + cs) * 0.5;
                param->b1 = a0r * -(1.0 + cs);
                param->b2 = a0r * (1.0 + cs) * 0.5;
                param->a1 = a0r * (2.0 * cs);
                param->a2 = a0r * (alpha - 1.0);
            }
    };

} // namespace LinuxSampler

#endif // __LS_BIQUADFILTER_H__
1	/***************************************************************************
2	* *
3	* LinuxSampler - modular, streaming capable sampler *
4	* *
5	* Copyright (C) 2003, 2004 by Benno Senoner and Christian Schoenebeck *
6	* *
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	* Internal parameters of the biquad filter, which are actually the
39	* final parameters of the filter's transfer function. This strucure is
40	* only needed when these parameters should stored outside the
41	* BiquadFilter class, e.g. to save calculation time by sharing them
42	* between multiple filters.
43	*/
44	struct biquad_param_t {
45	bq_t a1;
46	bq_t a2;
47	bq_t b0;
48	bq_t b1;
49	bq_t b2;
50	};
51
52	/**
53	* Bi-quadratic filter
54	* (adapted from lisp code by Eli Brandt, http://www.cs.cmu.edu/~eli/)
55	*/
56	class BiquadFilter {
57	protected:
58	// following five variables are only used if no external biquad_param_t reference is used
59	bq_t a1;
60	bq_t a2;
61	bq_t b0;
62	bq_t b1;
63	bq_t b2;
64	// following four variables are used to buffer the feedback
65	bq_t x1;
66	bq_t x2;
67	bq_t y1;
68	bq_t y2;
69
70	/**
71	* Prevent \a f from going into denormal mode which would slow down
72	* subsequent floating point calculations, we achieve that by setting
73	* \a f to zero when it falls under the denormal threshold value.
74	*/
75	inline void KillDenormal(bq_t& f) {
76	// TODO: this is a generic solution for 32bit floats, should be replaced by CPU specific asm code
77	f += 1e-18f;
78	f -= 1e-18f;
79	}
80	public:
81	inline BiquadFilter() {
82	x1 = 0.0f;
83	x2 = 0.0f;
84	y1 = 0.0f;
85	y2 = 0.0f;
86	}
87
88	inline bq_t Apply(const bq_t x) {
89	bq_t y;
90
91	y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
92	this->a1 * this->y1 + this->a2 * this->y2;
93	KillDenormal(y);
94	this->x2 = this->x1;
95	this->x1 = x;
96	this->y2 = this->y1;
97	this->y1 = y;
98
99	return y;
100	}
101
102	inline bq_t Apply(biquad_param_t* param, const bq_t x) {
103	bq_t y;
104
105	y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
106	param->a1 * this->y1 + param->a2 * this->y2;
107	KillDenormal(y);
108	this->x2 = this->x1;
109	this->x1 = x;
110	this->y2 = this->y1;
111	this->y1 = y;
112
113	return y;
114	}
115
116	inline bq_t ApplyFB(bq_t x, const bq_t fb) {
117	bq_t y;
118
119	x += this->y1 * fb * 0.98;
120	y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
121	this->a1 * this->y1 + this->a2 * this->y2;
122	KillDenormal(y);
123	this->x2 = this->x1;
124	this->x1 = x;
125	this->y2 = this->y1;
126	this->y1 = y;
127
128	return y;
129	}
130
131	inline bq_t ApplyFB(biquad_param_t* param, bq_t x, const bq_t fb) {
132	bq_t y;
133
134	x += this->y1 * fb * 0.98;
135	y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
136	param->a1 * this->y1 + param->a2 * this->y2;
137	KillDenormal(y);
138	this->x2 = this->x1;
139	this->x1 = x;
140	this->y2 = this->y1;
141	this->y1 = y;
142
143	return y;
144	}
145	};
146
147	class LowpassFilter : public BiquadFilter {
148	public:
149	inline LowpassFilter() : BiquadFilter() {}
150
151	inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
152	bq_t omega = 2.0 * M_PI * fc / fs;
153	bq_t sn = sin(omega);
154	bq_t cs = cos(omega);
155	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
156
157	const float a0r = 1.0 / (1.0 + alpha);
158	this->b0 = a0r * (1.0 - cs) * 0.5;
159	this->b1 = a0r * (1.0 - cs);
160	this->b2 = a0r * (1.0 - cs) * 0.5;
161	this->a1 = a0r * (2.0 * cs);
162	this->a2 = a0r * (alpha - 1.0);
163	}
164
165	inline void SetParameters(biquad_param_t* param, bq_t fc, bq_t bw, bq_t fs) {
166	bq_t omega = 2.0 * M_PI * fc / fs;
167	bq_t sn = sin(omega);
168	bq_t cs = cos(omega);
169	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
170
171	const float a0r = 1.0 / (1.0 + alpha);
172	param->b0 = a0r * (1.0 - cs) * 0.5;
173	param->b1 = a0r * (1.0 - cs);
174	param->b2 = a0r * (1.0 - cs) * 0.5;
175	param->a1 = a0r * (2.0 * cs);
176	param->a2 = a0r * (alpha - 1.0);
177	}
178	};
179
180	class BandpassFilter : public BiquadFilter {
181	public:
182	inline BandpassFilter() : BiquadFilter() {}
183
184	inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
185	bq_t omega = 2.0 * M_PI * fc / fs;
186	bq_t sn = sin(omega);
187	bq_t cs = cos(omega);
188	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
189
190	const float a0r = 1.0 / (1.0 + alpha);
191	this->b0 = a0r * alpha;
192	this->b1 = 0.0;
193	this->b2 = a0r * -alpha;
194	this->a1 = a0r * (2.0 * cs);
195	this->a2 = a0r * (alpha - 1.0);
196	}
197
198	inline void SetParameters(biquad_param_t* param, bq_t fc, bq_t bw, bq_t fs) {
199	bq_t omega = 2.0 * M_PI * fc / fs;
200	bq_t sn = sin(omega);
201	bq_t cs = cos(omega);
202	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
203
204	const float a0r = 1.0 / (1.0 + alpha);
205	param->b0 = a0r * alpha;
206	param->b1 = 0.0;
207	param->b2 = a0r * -alpha;
208	param->a1 = a0r * (2.0 * cs);
209	param->a2 = a0r * (alpha - 1.0);
210	}
211	};
212
213	class HighpassFilter : public BiquadFilter {
214	public:
215	inline HighpassFilter() : BiquadFilter() {}
216
217	inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
218	bq_t omega = 2.0 * M_PI * fc / fs;
219	bq_t sn = sin(omega);
220	bq_t cs = cos(omega);
221	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
222
223	const float a0r = 1.0 / (1.0 + alpha);
224	this->b0 = a0r * (1.0 + cs) * 0.5;
225	this->b1 = a0r * -(1.0 + cs);
226	this->b2 = a0r * (1.0 + cs) * 0.5;
227	this->a1 = a0r * (2.0 * cs);
228	this->a2 = a0r * (alpha - 1.0);
229	}
230
231	inline void SetParameters(biquad_param_t* param, bq_t fc, bq_t bw, bq_t fs) {
232	bq_t omega = 2.0 * M_PI * fc / fs;
233	bq_t sn = sin(omega);
234	bq_t cs = cos(omega);
235	bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
236
237	const float a0r = 1.0 / (1.0 + alpha);
238	param->b0 = a0r * (1.0 + cs) * 0.5;
239	param->b1 = a0r * -(1.0 + cs);
240	param->b2 = a0r * (1.0 + cs) * 0.5;
241	param->a1 = a0r * (2.0 * cs);
242	param->a2 = a0r * (alpha - 1.0);
243	}
244	};
245
246	} // namespace LinuxSampler
247
248	#endif // __LS_BIQUADFILTER_H__