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__

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