engines/common/BiquadFilter.h

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

// Note: the assembly code is currently disabled, as it doesn't fit into
// the new synthesis core introduced by LS 0.4.0

#ifndef __LS_BIQUADFILTER_H__
#define __LS_BIQUADFILTER_H__

#include <math.h>

#include "../../common/global_private.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 b0;
        bq_t b1;
        bq_t b2;
        bq_t a1;
        bq_t a2;
    };

    /**
     * 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 b0;
            bq_t b1;
            bq_t b2;
            bq_t a1;
            bq_t a2;
            // following four variables are used to buffer the feedback
            bq_t x1;
            bq_t x2;
            bq_t y1;
            bq_t y2;

#if __GNUC__ >= 4
            float fbc;
#else
            const static float fbc = 0.98;
#endif

            /**
             * 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:
            BiquadFilter() {
                Reset();
#if __GNUC__ >= 4
                fbc = 0.98f;
#endif
            }

            void Reset() {
                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(const biquad_param_t* __restrict 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;
            }

#if 0 // CONFIG_ASM && ARCH_X86
            // expects to find input in xmm0 (xmm0 stays unmodified) and finally leaves output in xmm6
            inline void Apply4StepsSSE(biquad_param_t* param) {
                __asm__ __volatile__ (
                    "movss (%2),%%xmm4                # b0\n\t"
                    "shufps   $0x00,%%xmm4,%%xmm4     # copy b0 to other cells\n\t"
                    "mulps  %%xmm0,%%xmm4             # xmm4 = x*b0\n\t"
                    "movups (%0),%%xmm2               # load b1,b2,a1,a2\n\t"
                    "movups (%1),%%xmm5               # load x1,x2,y1,y2\n\t"
                    /* sample 0 */
                    "movaps %%xmm5,%%xmm3\n\t"
                    "mulps  %%xmm2,%%xmm5             # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
                    "shufps $0x0a,%%xmm3,%%xmm3       # x2 = x1, y2 = y1\n\t"
                    "movss  %%xmm4,%%xmm6\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6             # xmm6 = b0*x + b1*x1 + b2*x2 + a1*y1 + a2*y2\n\t"
                    /* sample 1 */
                    "shufps $0x39,%%xmm4,%%xmm4       # rotate xmm4 down 1 cell\n\t"
                    "movss  %%xmm6,%%xmm3             # y1 = y\n\t"
                    "shufps $0x4e,%%xmm3,%%xmm3       # rotate 2 cells\n\t"
                    "movss  %%xmm0,%%xmm3             # x1 = x\n\t"
                    "shufps $0x93,%%xmm6,%%xmm6       # rotate output up 1 cell\n\t"
                    "movaps %%xmm3,%%xmm5\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0       # rotate input down 1 cell\n\t"
                    "mulps  %%xmm2,%%xmm5             # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
                    "movss  %%xmm5,%%xmm6\n\t"
                    "addss  %%xmm4,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6             # xmm6 = b0*x + b1*x1 + b2*x2 + a1*y1 + a2*y2\n\t"
                    /* sample 2 */
                    "shufps $0x0a,%%xmm3,%%xmm3       # x2 = x1, y2 = y1\n\t"
                    "shufps $0x39,%%xmm4,%%xmm4       # rotate xmm4 down 1 cell\n\t"
                    "movss  %%xmm6,%%xmm3             # y1 = y\n\t"
                    "shufps $0x4e,%%xmm3,%%xmm3       # rotate 2 cells\n\t"
                    "movss  %%xmm0,%%xmm3             # x1 = x\n\t"
                    "shufps $0x93,%%xmm6,%%xmm6       # rotate output up 1 cell\n\t"
                    "movaps %%xmm3,%%xmm5\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0       # rotate input down 1 cell\n\t"
                    "mulps  %%xmm2,%%xmm5             # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
                    "movss  %%xmm5,%%xmm6\n\t"
                    "addss  %%xmm4,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6\n\t"
                    "shufps $0x39,%%xmm5,%%xmm5\n\t"
                    "addss  %%xmm5,%%xmm6             # xmm6 = b0*x + b1*x1 + b2*x2 + a1*y1 + a2*y2\n\t"
                    /* sample 3 */
                    "shufps $0x0a,%%xmm3,%%xmm3       # x2 = x1, y2 = y1\n\t"
                    "shufps $0x39,%%xmm4,%%xmm4       # rotate xmm4 down 1 cell\n\t"
                    "movss  %%xmm6,%%xmm3             # y1 = y\n\t"
                    "shufps $0x4e,%%xmm3,%%xmm3       # rotate 2 cells\n\t"
                    "movss  %%xmm0,%%xmm3             # x1 = x\n\t"
                    "shufps $0x93,%%xmm6,%%xmm6       # rotate output up 1 cell\n\t"
                    "mulps  %%xmm3,%%xmm2             # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0       # rotate input down 1 cell\n\t"
                    "movss  %%xmm2,%%xmm6\n\t"
                    "shufps $0x39,%%xmm2,%%xmm2\n\t"
                    "addss  %%xmm2,%%xmm6\n\t"
                    "shufps $0x39,%%xmm2,%%xmm2\n\t"
                    "addss  %%xmm2,%%xmm6\n\t"
                    "shufps $0x39,%%xmm2,%%xmm2\n\t"
                    "addss  %%xmm2,%%xmm6\n\t"
                    "addss  %%xmm4,%%xmm6             # xmm6 = b0*x + b1*x1 + b2*x2 + a1*y1 + a2*y2\n\t"
                    /* done */
                    "shufps $0x0a,%%xmm3,%%xmm3       # x2 = x1, y2 = y1\n\t"
                    "movss  %%xmm6,%%xmm3             # y1 = y\n\t"
                    "shufps $0x4e,%%xmm3,%%xmm3       # rotate 2 cells\n\t"
                    "movss  %%xmm0,%%xmm3             # x1 = x\n\t"
                    "shufps $0x1b,%%xmm6,%%xmm6       # swap output to correct order\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0       # rotate input down 1 cell, to restore original input\n\t"
                    "movups %%xmm3,(%1)               # store x1,x2,y1,y2\n\t"
                    : /* no output */
                    : "r" (&param->b1), /* %0 - [b1,b2,a1,a2] */
                      "r" (&x1),        /* %1 - [x1,x2,y1,y2] */
                      "r" (&param->b0)  /* %2 */
                );
            }
#endif // CONFIG_ASM && ARCH_X86

            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(const biquad_param_t* __restrict 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;
            }

#if 0 // CONFIG_ASM && ARCH_X86
            // expects to find input in xmm0 (xmm0 stays unmodified) and finally leaves output in xmm7
            inline void ApplyFB4StepsSSE(biquad_param_t* param, const bq_t &fb) {
                float xs, ys;
                float t0, t1, t2, t3, t4, t5, t6, t7, t8; // temporary stack space
                __asm__ __volatile__ (
                    /* prepare input */
                    "movss  %15,%%xmm5\n\t"
                    "movss  %%xmm0,(%14)\n\t"
                    /* sample 0 */
                    "movss   %0, %%xmm3\n\t"
                    "movss   %1, %%xmm4\n\t"
                    "mulss   %%xmm4, %%xmm5\n\t"
                    "movss   %%xmm3, %2\n\t"
                    "movss   %%xmm5, %16\n\t"
                    "mulss   %%xmm3, %%xmm5\n\t"
                    "movss   %19, %%xmm2\n\t"
                    "movss   %3, %%xmm6\n\t"
                    "movss   %21, %%xmm3\n\t"
                    "addss   %%xmm5, %%xmm6\n\t"
                    "movss  %%xmm2, %%xmm5\n\t"
                    "movss   %20, %%xmm4\n\t"
                    "movss   %%xmm6, %4\n\t"
                    "mulss   %%xmm6, %%xmm5\n\t"
                    "movss   %5, %%xmm6\n\t"
                    "movss   %%xmm2, %6\n\t"
                    "movss   %%xmm4, %7\n\t"
                    "movss   %%xmm3, %%xmm2\n\t"
                    "mulss   %%xmm6, %%xmm4\n\t"
                    "mulss   %8, %%xmm2\n\t"
                    "movss   %%xmm3, %9\n\t"
                    "addss   %%xmm4, %%xmm5\n\t"
                    "movss   %18, %%xmm3\n\t"
                    "movss   %17, %%xmm4\n\t"
                    "addss   %%xmm2, %%xmm5\n\t"
                    "movss   %%xmm4, %10\n\t"
                    "movss   %%xmm3, %%xmm2\n\t"
                    "mulss   %11, %%xmm4\n\t"
                    "mulss   %12, %%xmm2\n\t"
                    "movss   %%xmm3, %13\n\t"
                    "addss   %%xmm4, %%xmm5\n\t"
                    "movss   %11, %%xmm3\n\t"
                    "movss   %4, %%xmm4\n\t"
                    "addss   %%xmm2, %%xmm5\n\t"
                    :: "m" (y1),  /* %0 */
                       "m" (fbc), /* %1 */
                       "m" (t0),  /* %2 */
                       "m" (xs),  /* %3 */
                       "m" (t7),  /* %4 */
                       "m" (x1),  /* %5 */
                       "m" (t1),  /* %6 */
                       "m" (t2),  /* %7 */
                       "m" (x2),  /* %8 */
                       "m" (t3),  /* %9 */
                       "m" (t4),  /* %10 */
                       "m" (t0),  /* %11 */
                       "m" (y2),  /* %12 */
                       "m" (t5),  /* %13 */
                       "r" (&xs),  /* %14 */
                       "m" (fb),   /* %15 */
                       "m" (ys),   /* %16 */
                       "m" (param->a1), /* %17 */
                       "m" (param->a2), /* %18 */
                       "m" (param->b0), /* %19 */
                       "m" (param->b1), /* %20 */
                       "m" (param->b2) /* %21 */
                );
                __asm__ __volatile__ (
                    "shufps $0x39,%%xmm0,%%xmm0  # rotate down one cell\n\t"
                    "movss  %%xmm5,%%xmm7\n\t"
                    ::
                );
                /* sample 1 */
                __asm__ __volatile__ (
                    "movss   %0, %%xmm4\n\t"
                    "movss   %%xmm0, %%xmm3\n\t"
                    "mulss   %%xmm5, %%xmm4\n\t"
                    "mulss   %3, %%xmm6\n\t"
                    "movss   %5, %%xmm2\n\t"
                    "addss   %%xmm4, %%xmm3\n\t"
                    "mulss   %7, %%xmm2\n\t"
                    "movss   %6, %%xmm4\n\t"
                    "movss   %%xmm3, %8\n\t"
                    "mulss   %%xmm3, %%xmm4\n\t"
                    "addss   %%xmm2, %%xmm4\n\t"
                    "movss   %9, %%xmm3\n\t"
                    "mulss   %%xmm5, %%xmm3\n\t"
                    "movss   %10, %%xmm2\n\t"
                    "addss   %%xmm6, %%xmm4\n\t"
                    "mulss   %11, %%xmm2\n\t"
                    "addss   %%xmm3, %%xmm4\n\t"
                    "addss   %%xmm2, %%xmm4\n\t"
                    :: "m" (ys),  /* %0 */
                       "m" (fbc), /* %1 */
                       "m" (xs),  /* %2 */
                       "m" (t3),  /* %3 */
                       "m" (y2),  /* %4 */
                       "m" (t2),  /* %5 */
                       "m" (t1),  /* %6 */
                       "m" (t7),  /* %7 */
                       "m" (t8),  /* %8 */
                       "m" (t4),  /* %9 */
                       "m" (t5),  /* %10 */
                       "m" (t0),  /* %11 */
                       "m" (x2),  /* %12 */
                       "m" (x1),  /* %13 */
                       "m" (y1)   /* %14 */
                );
                __asm__ __volatile__ (
                    "shufps $0x93,%%xmm7,%%xmm7  # rotate up one cell\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0  # rotate down one cell\n\t"
                    "movss  %%xmm4,%%xmm7\n\t"
                    ::
                );
                /* sample 2 */
                __asm__ __volatile__ (
                    "movss   %2, %%xmm6\n\t"
                    "movss   %3, %%xmm3\n\t"
                    "mulss   %%xmm4, %%xmm6\n\t"
                    "movss   %4, %%xmm2\n\t"
                    "mulss   %9, %%xmm2\n\t"
                    "addss   %%xmm0, %%xmm6\n\t"
                    "mulss   %7, %%xmm5\n\t"
                    "mulss   %%xmm6, %%xmm3\n\t"
                    "addss   %%xmm2, %%xmm3\n\t"
                    "movss   %5, %%xmm2\n\t"
                    "mulss   %8, %%xmm2\n\t"
                    "addss   %%xmm2, %%xmm3\n\t"
                    "movss   %6, %%xmm2\n\t"
                    "mulss   %%xmm4, %%xmm2\n\t"
                    "addss   %%xmm5, %%xmm2\n\t"
                    "addss   %%xmm2, %%xmm3\n\t"
                    :: "m" (xs),  /* %0 */
                       "m" (fb),  /* %1 */
                       "m" (ys), /* %2 */
                       "m" (t1),  /* %3 */
                       "m" (t2),  /* %4 */
                       "m" (t3),  /* %5 */
                       "m" (t4),  /* %6 */
                       "m" (t5),  /* %7 */
                       "m" (t7),  /* %8 */
                       "m" (t8),  /* %9 */
                       "m" (x1),  /* %10 */
                       "m" (x2),  /* %11 */
                       "m" (y1),  /* %12 */
                       "m" (y2)   /* %13 */
                );
                __asm__ __volatile__ (
                    "shufps $0x39,%%xmm0,%%xmm0  # rotate down one cell\n\t"
                    "shufps $0x93,%%xmm7,%%xmm7  # rotate up one cell\n\t"
                    "movss  %%xmm3,%%xmm7\n\t"
                    ::
                );
                /* sample 3 */
                __asm__ __volatile__ (
                    "movss   %1, %%xmm2\n\t"
                    "mulss   %7, %%xmm4\n\t"
                    "mulss   %%xmm3, %%xmm2\n\t"
                    "movss   %3, %%xmm5\n\t"
                    "movss   %%xmm6, %11\n\t"
                    "addss   %%xmm0, %%xmm2\n\t"
                    "movss   %%xmm3, %13\n\t"
                    "mulss   %%xmm2, %%xmm5\n\t"
                    "mulss   %4, %%xmm6\n\t"
                    "movss   %%xmm2, %10\n\t"
                    "addss   %%xmm6, %%xmm5\n\t"
                    "movss   %5, %%xmm2\n\t"
                    "mulss   %9, %%xmm2\n\t"
                    "mulss   %6, %%xmm3\n\t"
                    "addss   %%xmm2, %%xmm5\n\t"
                    "addss   %%xmm3, %%xmm4\n\t"
                    "addss   %%xmm4, %%xmm5\n\t"
                    "movss   %%xmm5, %12\n\t"
                    :: "m" (xs),  /* %0 */
                       "m" (ys),  /* %1 */
                       "m" (fbc), /* %2 */
                       "m" (t1),  /* %3 */
                       "m" (t2),  /* %4 */
                       "m" (t3),  /* %5 */
                       "m" (t4),  /* %6 */
                       "m" (t5),  /* %7 */
                       "m" (t6),  /* %8 */
                       "m" (t8),  /* %9 */
                       "m" (x1),  /* %10 */
                       "m" (x2),  /* %11 */
                       "m" (y1),  /* %12 */
                       "m" (y2)   /* %13 */
                );
                __asm__ __volatile__ (
                    "shufps $0x93,%%xmm7,%%xmm7  # rotate up one cell\n\t"
                    "shufps $0x39,%%xmm0,%%xmm0  # rotate down one cell to restore original input\n\t"
                    "movss  %%xmm5,%%xmm7\n\t"
                    "shufps $0x1b,%%xmm7,%%xmm7  # swap output to correct order\n\t"
                    ::
                );
            }
#endif // CONFIG_ASM && ARCH_X86
    };

    /** @brief Lowpass Filter
     *
     * Lowpass filter based on biquad filter implementation.
     */
    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* __restrict 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);
            }
    };

    /** @brief Bandpass Filter
     *
     * Bandpass filter based on biquad filter implementation.
     */
    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 * sn * 0.71;
                this->b1 = 0.0;
                this->b2 = a0r * -sn * 0.71;
                this->a1 = a0r * (2.0 * cs);
                this->a2 = a0r * (alpha - 1.0);
            }

            inline void SetParameters(biquad_param_t* __restrict 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 * sn * 0.71;
                param->b1 = 0.0;
                param->b2 = a0r * -sn * 0.71;
                param->a1 = a0r * (2.0 * cs);
                param->a2 = a0r * (alpha - 1.0);
            }
    };

    /** @brief Highpass Filter
     *
     * Highpass filter based on biquad filter implementation.
     */
    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* __restrict 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	2377	* Copyright (C) 2005 - 2012 Christian Schoenebeck *
7	schoenebeck	53	* *
8			* This program is free software; you can redistribute it and/or modify *
9			* it under the terms of the GNU General Public License as published by *
10			* the Free Software Foundation; either version 2 of the License, or *
11			* (at your option) any later version. *
12			* *
13			* This program is distributed in the hope that it will be useful, *
14			* but WITHOUT ANY WARRANTY; without even the implied warranty of *
15			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16			* GNU General Public License for more details. *
17			* *
18			* You should have received a copy of the GNU General Public License *
19			* along with this program; if not, write to the Free Software *
20			* Foundation, Inc., 59 Temple Place, Suite 330, Boston, *
21			* MA 02111-1307 USA *
22			***************************************************************************/
23
24	schoenebeck	1259	// Note: the assembly code is currently disabled, as it doesn't fit into
25			// the new synthesis core introduced by LS 0.4.0
26
27	schoenebeck	53	#ifndef __LS_BIQUADFILTER_H__
28			#define __LS_BIQUADFILTER_H__
29
30	capela	268	#include <math.h>
31
32	schoenebeck	1424	#include "../../common/global_private.h"
33	schoenebeck	328
34	schoenebeck	53	/// ln(2) / 2
35			#define LN_2_2 0.34657359f
36
37			#ifndef LIMIT
38			# define LIMIT(v,l,u) (v < l ? l : (v > u ? u : v))
39			#endif
40
41			namespace LinuxSampler {
42
43			typedef float bq_t;
44
45			/**
46	schoenebeck	80	* Internal parameters of the biquad filter, which are actually the
47			* final parameters of the filter's transfer function. This strucure is
48			* only needed when these parameters should stored outside the
49			* BiquadFilter class, e.g. to save calculation time by sharing them
50			* between multiple filters.
51			*/
52			struct biquad_param_t {
53			bq_t b0;
54			bq_t b1;
55			bq_t b2;
56	senkov	324	bq_t a1;
57			bq_t a2;
58	schoenebeck	80	};
59
60			/**
61			* Bi-quadratic filter
62			* (adapted from lisp code by Eli Brandt, http://www.cs.cmu.edu/~eli/)
63			*/
64	schoenebeck	53	class BiquadFilter {
65			protected:
66	schoenebeck	80	// following five variables are only used if no external biquad_param_t reference is used
67	schoenebeck	53	bq_t b0;
68			bq_t b1;
69			bq_t b2;
70	senkov	324	bq_t a1;
71			bq_t a2;
72	schoenebeck	80	// following four variables are used to buffer the feedback
73	schoenebeck	53	bq_t x1;
74			bq_t x2;
75			bq_t y1;
76			bq_t y2;
77
78	persson	685	#if __GNUC__ >= 4
79			float fbc;
80			#else
81	schoenebeck	319	const static float fbc = 0.98;
82	persson	685	#endif
83	schoenebeck	319
84	schoenebeck	53	/**
85	schoenebeck	80	* Prevent \a f from going into denormal mode which would slow down
86			* subsequent floating point calculations, we achieve that by setting
87			* \a f to zero when it falls under the denormal threshold value.
88			*/
89			inline void KillDenormal(bq_t& f) {
90	schoenebeck	53	// TODO: this is a generic solution for 32bit floats, should be replaced by CPU specific asm code
91			f += 1e-18f;
92			f -= 1e-18f;
93			}
94			public:
95	schoenebeck	319	BiquadFilter() {
96			Reset();
97	persson	685	#if __GNUC__ >= 4
98			fbc = 0.98f;
99			#endif
100	schoenebeck	319	}
101
102			void Reset() {
103	schoenebeck	53	x1 = 0.0f;
104			x2 = 0.0f;
105			y1 = 0.0f;
106			y2 = 0.0f;
107			}
108
109			inline bq_t Apply(const bq_t x) {
110			bq_t y;
111
112			y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
113			this->a1 * this->y1 + this->a2 * this->y2;
114			KillDenormal(y);
115			this->x2 = this->x1;
116			this->x1 = x;
117			this->y2 = this->y1;
118			this->y1 = y;
119
120			return y;
121			}
122
123	schoenebeck	2377	inline bq_t Apply(const biquad_param_t* __restrict param, const bq_t x) {
124	schoenebeck	80	bq_t y;
125
126			y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
127			param->a1 * this->y1 + param->a2 * this->y2;
128			KillDenormal(y);
129			this->x2 = this->x1;
130			this->x1 = x;
131			this->y2 = this->y1;
132			this->y1 = y;
133
134			return y;
135			}
136
137	schoenebeck	1259	#if 0 // CONFIG_ASM && ARCH_X86
138	schoenebeck	319	// expects to find input in xmm0 (xmm0 stays unmodified) and finally leaves output in xmm6
139			inline void Apply4StepsSSE(biquad_param_t* param) {
140			__asm__ __volatile__ (
141			"movss (%2),%%xmm4 # b0\n\t"
142			"shufps $0x00,%%xmm4,%%xmm4 # copy b0 to other cells\n\t"
143			"mulps %%xmm0,%%xmm4 # xmm4 = x*b0\n\t"
144			"movups (%0),%%xmm2 # load b1,b2,a1,a2\n\t"
145			"movups (%1),%%xmm5 # load x1,x2,y1,y2\n\t"
146			/* sample 0 */
147			"movaps %%xmm5,%%xmm3\n\t"
148			"mulps %%xmm2,%%xmm5 # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
149			"shufps $0x0a,%%xmm3,%%xmm3 # x2 = x1, y2 = y1\n\t"
150			"movss %%xmm4,%%xmm6\n\t"
151			"addss %%xmm5,%%xmm6\n\t"
152			"shufps $0x39,%%xmm5,%%xmm5\n\t"
153			"addss %%xmm5,%%xmm6\n\t"
154			"shufps $0x39,%%xmm5,%%xmm5\n\t"
155			"addss %%xmm5,%%xmm6\n\t"
156			"shufps $0x39,%%xmm5,%%xmm5\n\t"
157			"addss %%xmm5,%%xmm6 # xmm6 = b0x + b1x1 + b2x2 + a1y1 + a2*y2\n\t"
158			/* sample 1 */
159			"shufps $0x39,%%xmm4,%%xmm4 # rotate xmm4 down 1 cell\n\t"
160			"movss %%xmm6,%%xmm3 # y1 = y\n\t"
161			"shufps $0x4e,%%xmm3,%%xmm3 # rotate 2 cells\n\t"
162			"movss %%xmm0,%%xmm3 # x1 = x\n\t"
163			"shufps $0x93,%%xmm6,%%xmm6 # rotate output up 1 cell\n\t"
164			"movaps %%xmm3,%%xmm5\n\t"
165			"shufps $0x39,%%xmm0,%%xmm0 # rotate input down 1 cell\n\t"
166			"mulps %%xmm2,%%xmm5 # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
167			"movss %%xmm5,%%xmm6\n\t"
168			"addss %%xmm4,%%xmm6\n\t"
169			"shufps $0x39,%%xmm5,%%xmm5\n\t"
170			"addss %%xmm5,%%xmm6\n\t"
171			"shufps $0x39,%%xmm5,%%xmm5\n\t"
172			"addss %%xmm5,%%xmm6\n\t"
173			"shufps $0x39,%%xmm5,%%xmm5\n\t"
174			"addss %%xmm5,%%xmm6 # xmm6 = b0x + b1x1 + b2x2 + a1y1 + a2*y2\n\t"
175			/* sample 2 */
176			"shufps $0x0a,%%xmm3,%%xmm3 # x2 = x1, y2 = y1\n\t"
177			"shufps $0x39,%%xmm4,%%xmm4 # rotate xmm4 down 1 cell\n\t"
178			"movss %%xmm6,%%xmm3 # y1 = y\n\t"
179			"shufps $0x4e,%%xmm3,%%xmm3 # rotate 2 cells\n\t"
180			"movss %%xmm0,%%xmm3 # x1 = x\n\t"
181			"shufps $0x93,%%xmm6,%%xmm6 # rotate output up 1 cell\n\t"
182			"movaps %%xmm3,%%xmm5\n\t"
183			"shufps $0x39,%%xmm0,%%xmm0 # rotate input down 1 cell\n\t"
184			"mulps %%xmm2,%%xmm5 # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
185			"movss %%xmm5,%%xmm6\n\t"
186			"addss %%xmm4,%%xmm6\n\t"
187			"shufps $0x39,%%xmm5,%%xmm5\n\t"
188			"addss %%xmm5,%%xmm6\n\t"
189			"shufps $0x39,%%xmm5,%%xmm5\n\t"
190			"addss %%xmm5,%%xmm6\n\t"
191			"shufps $0x39,%%xmm5,%%xmm5\n\t"
192			"addss %%xmm5,%%xmm6 # xmm6 = b0x + b1x1 + b2x2 + a1y1 + a2*y2\n\t"
193			/* sample 3 */
194			"shufps $0x0a,%%xmm3,%%xmm3 # x2 = x1, y2 = y1\n\t"
195			"shufps $0x39,%%xmm4,%%xmm4 # rotate xmm4 down 1 cell\n\t"
196			"movss %%xmm6,%%xmm3 # y1 = y\n\t"
197			"shufps $0x4e,%%xmm3,%%xmm3 # rotate 2 cells\n\t"
198			"movss %%xmm0,%%xmm3 # x1 = x\n\t"
199			"shufps $0x93,%%xmm6,%%xmm6 # rotate output up 1 cell\n\t"
200			"mulps %%xmm3,%%xmm2 # xmm5 = [b1,b2,a1,a2] * [x1,x2,y1,y2]\n\t"
201			"shufps $0x39,%%xmm0,%%xmm0 # rotate input down 1 cell\n\t"
202			"movss %%xmm2,%%xmm6\n\t"
203			"shufps $0x39,%%xmm2,%%xmm2\n\t"
204			"addss %%xmm2,%%xmm6\n\t"
205			"shufps $0x39,%%xmm2,%%xmm2\n\t"
206			"addss %%xmm2,%%xmm6\n\t"
207			"shufps $0x39,%%xmm2,%%xmm2\n\t"
208			"addss %%xmm2,%%xmm6\n\t"
209			"addss %%xmm4,%%xmm6 # xmm6 = b0x + b1x1 + b2x2 + a1y1 + a2*y2\n\t"
210			/* done */
211			"shufps $0x0a,%%xmm3,%%xmm3 # x2 = x1, y2 = y1\n\t"
212			"movss %%xmm6,%%xmm3 # y1 = y\n\t"
213			"shufps $0x4e,%%xmm3,%%xmm3 # rotate 2 cells\n\t"
214			"movss %%xmm0,%%xmm3 # x1 = x\n\t"
215			"shufps $0x1b,%%xmm6,%%xmm6 # swap output to correct order\n\t"
216			"shufps $0x39,%%xmm0,%%xmm0 # rotate input down 1 cell, to restore original input\n\t"
217			"movups %%xmm3,(%1) # store x1,x2,y1,y2\n\t"
218			: /* no output */
219			: "r" (&param->b1), /* %0 - [b1,b2,a1,a2] */
220			"r" (&x1), /* %1 - [x1,x2,y1,y2] */
221			"r" (&param->b0) /* %2 */
222			);
223			}
224	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
225	schoenebeck	319
226	schoenebeck	53	inline bq_t ApplyFB(bq_t x, const bq_t fb) {
227			bq_t y;
228
229			x += this->y1 * fb * 0.98;
230			y = this->b0 * x + this->b1 * this->x1 + this->b2 * this->x2 +
231			this->a1 * this->y1 + this->a2 * this->y2;
232			KillDenormal(y);
233			this->x2 = this->x1;
234			this->x1 = x;
235			this->y2 = this->y1;
236			this->y1 = y;
237
238			return y;
239			}
240	schoenebeck	80
241	schoenebeck	2377	inline bq_t ApplyFB(const biquad_param_t* __restrict param, bq_t x, const bq_t fb) {
242	schoenebeck	80	bq_t y;
243
244			x += this->y1 * fb * 0.98;
245			y = param->b0 * x + param->b1 * this->x1 + param->b2 * this->x2 +
246			param->a1 * this->y1 + param->a2 * this->y2;
247			KillDenormal(y);
248			this->x2 = this->x1;
249			this->x1 = x;
250			this->y2 = this->y1;
251			this->y1 = y;
252
253			return y;
254			}
255	schoenebeck	319
256	schoenebeck	1259	#if 0 // CONFIG_ASM && ARCH_X86
257	schoenebeck	319	// expects to find input in xmm0 (xmm0 stays unmodified) and finally leaves output in xmm7
258			inline void ApplyFB4StepsSSE(biquad_param_t* param, const bq_t &fb) {
259			float xs, ys;
260			float t0, t1, t2, t3, t4, t5, t6, t7, t8; // temporary stack space
261			__asm__ __volatile__ (
262			/* prepare input */
263			"movss %15,%%xmm5\n\t"
264			"movss %%xmm0,(%14)\n\t"
265			/* sample 0 */
266			"movss %0, %%xmm3\n\t"
267			"movss %1, %%xmm4\n\t"
268			"mulss %%xmm4, %%xmm5\n\t"
269			"movss %%xmm3, %2\n\t"
270			"movss %%xmm5, %16\n\t"
271			"mulss %%xmm3, %%xmm5\n\t"
272			"movss %19, %%xmm2\n\t"
273			"movss %3, %%xmm6\n\t"
274			"movss %21, %%xmm3\n\t"
275			"addss %%xmm5, %%xmm6\n\t"
276			"movss %%xmm2, %%xmm5\n\t"
277			"movss %20, %%xmm4\n\t"
278			"movss %%xmm6, %4\n\t"
279			"mulss %%xmm6, %%xmm5\n\t"
280			"movss %5, %%xmm6\n\t"
281			"movss %%xmm2, %6\n\t"
282			"movss %%xmm4, %7\n\t"
283			"movss %%xmm3, %%xmm2\n\t"
284			"mulss %%xmm6, %%xmm4\n\t"
285			"mulss %8, %%xmm2\n\t"
286			"movss %%xmm3, %9\n\t"
287			"addss %%xmm4, %%xmm5\n\t"
288			"movss %18, %%xmm3\n\t"
289			"movss %17, %%xmm4\n\t"
290			"addss %%xmm2, %%xmm5\n\t"
291			"movss %%xmm4, %10\n\t"
292			"movss %%xmm3, %%xmm2\n\t"
293			"mulss %11, %%xmm4\n\t"
294			"mulss %12, %%xmm2\n\t"
295			"movss %%xmm3, %13\n\t"
296			"addss %%xmm4, %%xmm5\n\t"
297			"movss %11, %%xmm3\n\t"
298			"movss %4, %%xmm4\n\t"
299			"addss %%xmm2, %%xmm5\n\t"
300			:: "m" (y1), /* %0 */
301			"m" (fbc), /* %1 */
302			"m" (t0), /* %2 */
303			"m" (xs), /* %3 */
304			"m" (t7), /* %4 */
305			"m" (x1), /* %5 */
306			"m" (t1), /* %6 */
307			"m" (t2), /* %7 */
308			"m" (x2), /* %8 */
309			"m" (t3), /* %9 */
310			"m" (t4), /* %10 */
311			"m" (t0), /* %11 */
312			"m" (y2), /* %12 */
313			"m" (t5), /* %13 */
314			"r" (&xs), /* %14 */
315			"m" (fb), /* %15 */
316			"m" (ys), /* %16 */
317			"m" (param->a1), /* %17 */
318			"m" (param->a2), /* %18 */
319			"m" (param->b0), /* %19 */
320			"m" (param->b1), /* %20 */
321			"m" (param->b2) /* %21 */
322			);
323			__asm__ __volatile__ (
324			"shufps $0x39,%%xmm0,%%xmm0 # rotate down one cell\n\t"
325			"movss %%xmm5,%%xmm7\n\t"
326			::
327			);
328			/* sample 1 */
329			__asm__ __volatile__ (
330			"movss %0, %%xmm4\n\t"
331			"movss %%xmm0, %%xmm3\n\t"
332			"mulss %%xmm5, %%xmm4\n\t"
333			"mulss %3, %%xmm6\n\t"
334			"movss %5, %%xmm2\n\t"
335			"addss %%xmm4, %%xmm3\n\t"
336			"mulss %7, %%xmm2\n\t"
337			"movss %6, %%xmm4\n\t"
338			"movss %%xmm3, %8\n\t"
339			"mulss %%xmm3, %%xmm4\n\t"
340			"addss %%xmm2, %%xmm4\n\t"
341			"movss %9, %%xmm3\n\t"
342			"mulss %%xmm5, %%xmm3\n\t"
343			"movss %10, %%xmm2\n\t"
344			"addss %%xmm6, %%xmm4\n\t"
345			"mulss %11, %%xmm2\n\t"
346			"addss %%xmm3, %%xmm4\n\t"
347			"addss %%xmm2, %%xmm4\n\t"
348			:: "m" (ys), /* %0 */
349			"m" (fbc), /* %1 */
350			"m" (xs), /* %2 */
351			"m" (t3), /* %3 */
352			"m" (y2), /* %4 */
353			"m" (t2), /* %5 */
354			"m" (t1), /* %6 */
355			"m" (t7), /* %7 */
356			"m" (t8), /* %8 */
357			"m" (t4), /* %9 */
358			"m" (t5), /* %10 */
359			"m" (t0), /* %11 */
360			"m" (x2), /* %12 */
361			"m" (x1), /* %13 */
362			"m" (y1) /* %14 */
363			);
364			__asm__ __volatile__ (
365			"shufps $0x93,%%xmm7,%%xmm7 # rotate up one cell\n\t"
366			"shufps $0x39,%%xmm0,%%xmm0 # rotate down one cell\n\t"
367			"movss %%xmm4,%%xmm7\n\t"
368			::
369			);
370			/* sample 2 */
371			__asm__ __volatile__ (
372			"movss %2, %%xmm6\n\t"
373			"movss %3, %%xmm3\n\t"
374			"mulss %%xmm4, %%xmm6\n\t"
375			"movss %4, %%xmm2\n\t"
376			"mulss %9, %%xmm2\n\t"
377			"addss %%xmm0, %%xmm6\n\t"
378			"mulss %7, %%xmm5\n\t"
379			"mulss %%xmm6, %%xmm3\n\t"
380			"addss %%xmm2, %%xmm3\n\t"
381			"movss %5, %%xmm2\n\t"
382			"mulss %8, %%xmm2\n\t"
383			"addss %%xmm2, %%xmm3\n\t"
384			"movss %6, %%xmm2\n\t"
385			"mulss %%xmm4, %%xmm2\n\t"
386			"addss %%xmm5, %%xmm2\n\t"
387			"addss %%xmm2, %%xmm3\n\t"
388			:: "m" (xs), /* %0 */
389			"m" (fb), /* %1 */
390			"m" (ys), /* %2 */
391			"m" (t1), /* %3 */
392			"m" (t2), /* %4 */
393			"m" (t3), /* %5 */
394			"m" (t4), /* %6 */
395			"m" (t5), /* %7 */
396			"m" (t7), /* %8 */
397			"m" (t8), /* %9 */
398			"m" (x1), /* %10 */
399			"m" (x2), /* %11 */
400			"m" (y1), /* %12 */
401			"m" (y2) /* %13 */
402			);
403			__asm__ __volatile__ (
404			"shufps $0x39,%%xmm0,%%xmm0 # rotate down one cell\n\t"
405			"shufps $0x93,%%xmm7,%%xmm7 # rotate up one cell\n\t"
406			"movss %%xmm3,%%xmm7\n\t"
407			::
408			);
409			/* sample 3 */
410			__asm__ __volatile__ (
411			"movss %1, %%xmm2\n\t"
412			"mulss %7, %%xmm4\n\t"
413			"mulss %%xmm3, %%xmm2\n\t"
414			"movss %3, %%xmm5\n\t"
415			"movss %%xmm6, %11\n\t"
416			"addss %%xmm0, %%xmm2\n\t"
417			"movss %%xmm3, %13\n\t"
418			"mulss %%xmm2, %%xmm5\n\t"
419			"mulss %4, %%xmm6\n\t"
420			"movss %%xmm2, %10\n\t"
421			"addss %%xmm6, %%xmm5\n\t"
422			"movss %5, %%xmm2\n\t"
423			"mulss %9, %%xmm2\n\t"
424			"mulss %6, %%xmm3\n\t"
425			"addss %%xmm2, %%xmm5\n\t"
426			"addss %%xmm3, %%xmm4\n\t"
427			"addss %%xmm4, %%xmm5\n\t"
428			"movss %%xmm5, %12\n\t"
429			:: "m" (xs), /* %0 */
430			"m" (ys), /* %1 */
431			"m" (fbc), /* %2 */
432			"m" (t1), /* %3 */
433			"m" (t2), /* %4 */
434			"m" (t3), /* %5 */
435			"m" (t4), /* %6 */
436			"m" (t5), /* %7 */
437			"m" (t6), /* %8 */
438			"m" (t8), /* %9 */
439			"m" (x1), /* %10 */
440			"m" (x2), /* %11 */
441			"m" (y1), /* %12 */
442			"m" (y2) /* %13 */
443			);
444			__asm__ __volatile__ (
445			"shufps $0x93,%%xmm7,%%xmm7 # rotate up one cell\n\t"
446			"shufps $0x39,%%xmm0,%%xmm0 # rotate down one cell to restore original input\n\t"
447			"movss %%xmm5,%%xmm7\n\t"
448			"shufps $0x1b,%%xmm7,%%xmm7 # swap output to correct order\n\t"
449			::
450			);
451			}
452	schoenebeck	617	#endif // CONFIG_ASM && ARCH_X86
453	schoenebeck	53	};
454
455	schoenebeck	563	/** @brief Lowpass Filter
456			*
457			* Lowpass filter based on biquad filter implementation.
458			*/
459	schoenebeck	53	class LowpassFilter : public BiquadFilter {
460			public:
461			inline LowpassFilter() : BiquadFilter() {}
462
463			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
464			bq_t omega = 2.0 * M_PI * fc / fs;
465			bq_t sn = sin(omega);
466			bq_t cs = cos(omega);
467			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
468
469			const float a0r = 1.0 / (1.0 + alpha);
470			this->b0 = a0r * (1.0 - cs) * 0.5;
471			this->b1 = a0r * (1.0 - cs);
472			this->b2 = a0r * (1.0 - cs) * 0.5;
473			this->a1 = a0r * (2.0 * cs);
474			this->a2 = a0r * (alpha - 1.0);
475			}
476	schoenebeck	80
477	schoenebeck	2377	inline void SetParameters(biquad_param_t* __restrict param, bq_t fc, bq_t bw, bq_t fs) {
478	schoenebeck	80	bq_t omega = 2.0 * M_PI * fc / fs;
479			bq_t sn = sin(omega);
480			bq_t cs = cos(omega);
481			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
482
483			const float a0r = 1.0 / (1.0 + alpha);
484			param->b0 = a0r * (1.0 - cs) * 0.5;
485			param->b1 = a0r * (1.0 - cs);
486			param->b2 = a0r * (1.0 - cs) * 0.5;
487			param->a1 = a0r * (2.0 * cs);
488			param->a2 = a0r * (alpha - 1.0);
489			}
490	schoenebeck	53	};
491
492	schoenebeck	563	/** @brief Bandpass Filter
493			*
494			* Bandpass filter based on biquad filter implementation.
495			*/
496	schoenebeck	53	class BandpassFilter : public BiquadFilter {
497			public:
498			inline BandpassFilter() : BiquadFilter() {}
499
500			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
501			bq_t omega = 2.0 * M_PI * fc / fs;
502			bq_t sn = sin(omega);
503			bq_t cs = cos(omega);
504			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
505
506			const float a0r = 1.0 / (1.0 + alpha);
507	persson	791	this->b0 = a0r * sn * 0.71;
508	schoenebeck	53	this->b1 = 0.0;
509	persson	791	this->b2 = a0r * -sn * 0.71;
510	schoenebeck	53	this->a1 = a0r * (2.0 * cs);
511			this->a2 = a0r * (alpha - 1.0);
512			}
513	schoenebeck	80
514	schoenebeck	2377	inline void SetParameters(biquad_param_t* __restrict param, bq_t fc, bq_t bw, bq_t fs) {
515	schoenebeck	80	bq_t omega = 2.0 * M_PI * fc / fs;
516			bq_t sn = sin(omega);
517			bq_t cs = cos(omega);
518			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
519
520			const float a0r = 1.0 / (1.0 + alpha);
521	persson	791	param->b0 = a0r * sn * 0.71;
522	schoenebeck	80	param->b1 = 0.0;
523	persson	791	param->b2 = a0r * -sn * 0.71;
524	schoenebeck	80	param->a1 = a0r * (2.0 * cs);
525			param->a2 = a0r * (alpha - 1.0);
526			}
527	schoenebeck	53	};
528
529	schoenebeck	563	/** @brief Highpass Filter
530			*
531			* Highpass filter based on biquad filter implementation.
532			*/
533	schoenebeck	53	class HighpassFilter : public BiquadFilter {
534			public:
535			inline HighpassFilter() : BiquadFilter() {}
536
537			inline void SetParameters(bq_t fc, bq_t bw, bq_t fs) {
538			bq_t omega = 2.0 * M_PI * fc / fs;
539			bq_t sn = sin(omega);
540			bq_t cs = cos(omega);
541			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
542
543			const float a0r = 1.0 / (1.0 + alpha);
544			this->b0 = a0r * (1.0 + cs) * 0.5;
545			this->b1 = a0r * -(1.0 + cs);
546			this->b2 = a0r * (1.0 + cs) * 0.5;
547			this->a1 = a0r * (2.0 * cs);
548			this->a2 = a0r * (alpha - 1.0);
549			}
550	schoenebeck	80
551	schoenebeck	2377	inline void SetParameters(biquad_param_t* __restrict param, bq_t fc, bq_t bw, bq_t fs) {
552	schoenebeck	80	bq_t omega = 2.0 * M_PI * fc / fs;
553			bq_t sn = sin(omega);
554			bq_t cs = cos(omega);
555			bq_t alpha = sn * sinh(M_LN2 / 2.0 * bw * omega / sn);
556
557			const float a0r = 1.0 / (1.0 + alpha);
558			param->b0 = a0r * (1.0 + cs) * 0.5;
559			param->b1 = a0r * -(1.0 + cs);
560			param->b2 = a0r * (1.0 + cs) * 0.5;
561			param->a1 = a0r * (2.0 * cs);
562			param->a2 = a0r * (alpha - 1.0);
563			}
564	schoenebeck	53	};
565
566			} // namespace LinuxSampler
567
568			#endif // __LS_BIQUADFILTER_H__