413 lines
11 KiB
C++
413 lines
11 KiB
C++
/*
|
|
WDL - besselfilter.h
|
|
(c) Theo Niessink 2011
|
|
<http://www.taletn.com/>
|
|
|
|
This software is provided 'as-is', without any express or implied
|
|
warranty. In no event will the authors be held liable for any damages
|
|
arising from the use of this software.
|
|
|
|
Permission is granted to anyone to use this software for any purpose,
|
|
including commercial applications, and to alter it and redistribute it
|
|
freely, subject to the following restrictions:
|
|
|
|
1. The origin of this software must not be misrepresented; you must not
|
|
claim that you wrote the original software. If you use this software
|
|
in a product, an acknowledgment in the product documentation would be
|
|
appreciated but is not required.
|
|
2. Altered source versions must be plainly marked as such, and must not be
|
|
misrepresented as being the original software.
|
|
3. This notice may not be removed or altered from any source distribution.
|
|
|
|
|
|
This file provides classes for a low-pass Bessel filter design using the
|
|
matched Z-transform method.
|
|
|
|
This Bessel filter implementation was originally extracted from from the
|
|
source code of mkfilter, written by A.J. Fisher.
|
|
<http://www-users.cs.york.ac.uk/~fisher/mkfilter>
|
|
|
|
|
|
Example #1:
|
|
|
|
// 8th order anti-alias filter
|
|
#define WDL_BESSEL_FILTER_ORDER 8
|
|
#include "besselfilter.h"
|
|
|
|
int oversampling = 8;
|
|
|
|
WDL_BesselFilterCoeffs bessel;
|
|
WDL_BesselFilterStage filter;
|
|
|
|
bessel.Calc(0.5 / (double)oversampling);
|
|
filter.Reset();
|
|
|
|
for (int i = 0; i < nFrames; ++i)
|
|
{
|
|
filter.Process(inputs[0][i], bessel.Coeffs());
|
|
outputs[0][i] = filter.Output();
|
|
}
|
|
|
|
Example #2:
|
|
|
|
#include "besselfilter.h"
|
|
|
|
int order = 4;
|
|
int oversampling = 8;
|
|
|
|
// 2 cascaded filters
|
|
WDL_BesselFilterStage filter[2];
|
|
filter[0].Reset();
|
|
filter[1].Reset();
|
|
|
|
WDL_BesselFilterCoeffs coeffs;
|
|
coeffs.Calc(0.5 / (double)oversampling, order);
|
|
|
|
for (int i = 0; i < nFrames; ++i)
|
|
{
|
|
filter[0].Process(inputs[0][i], &coeffs);
|
|
filter[1].Process(filter[0].Output(), &coeffs);
|
|
outputs[0][i] = filter[1].Output();
|
|
}
|
|
|
|
Example #3:
|
|
|
|
#define WDL_BESSEL_DENORMAL_AGGRESSIVE
|
|
#include "besselfilter.h"
|
|
|
|
int order = 8;
|
|
int oversampling = 8;
|
|
|
|
WDL_BesselFilter bessel;
|
|
bessel.Calc(0.5 / (double)oversampling, order);
|
|
bessel.Reset();
|
|
|
|
for (int i = 0; i < nFrames; ++i)
|
|
{
|
|
bessel.Process(inputs[0][i]);
|
|
outputs[0][i] = bessel.Output();
|
|
}
|
|
|
|
*/
|
|
|
|
|
|
#ifndef _BESSELFILTER_H_
|
|
#define _BESSELFILTER_H_
|
|
|
|
|
|
#include <math.h>
|
|
#ifdef _MSC_VER
|
|
#pragma warning(disable:4996) // hypot
|
|
#endif
|
|
|
|
#include <string.h>
|
|
#include <assert.h>
|
|
|
|
#include "wdltypes.h"
|
|
|
|
// By default denormals are zeroed to prevent exessive CPU use. Defining
|
|
// WDL_BESSEL_DENORMAL_IGNORE will disable denormal filtering. Defining
|
|
// WDL_BESSEL_DENORMAL_AGGRESSIVE will filter out denormals more
|
|
// aggressively by zeroing anything below 5.6e-017.
|
|
#ifndef WDL_BESSEL_DENORMAL_IGNORE
|
|
#include "denormal.h"
|
|
#if defined(WDL_BESSEL_DENORMAL_AGGRESSIVE)
|
|
#define WDL_BESSEL_FIX_DENORMAL(a) (denormal_fix_double_aggressive(a))
|
|
#else
|
|
#define WDL_BESSEL_FIX_DENORMAL(a) (denormal_fix_double(a))
|
|
#endif
|
|
#else
|
|
#define WDL_BESSEL_FIX_DENORMAL(a) ((void)0)
|
|
#endif
|
|
|
|
|
|
// Defining WDL_BESSEL_FILTER_ORDER will make the Bessel filter order fixed,
|
|
// which increases efficiency.
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
#if !(WDL_BESSEL_FILTER_ORDER >= 1 && WDL_BESSEL_FILTER_ORDER <= 10)
|
|
#error WDL_BESSEL_FILTER_ORDER should be in 1..10 range
|
|
#endif
|
|
#define WDL_BESSEL_FILTER_MAX WDL_BESSEL_FILTER_ORDER
|
|
|
|
// Defining WDL_BESSEL_FILTER_MAX limits the maximum Bessel filter order,
|
|
// which reduces buffer sizes.
|
|
#else
|
|
#if defined(WDL_BESSEL_FILTER_MAX) && WDL_BESSEL_FILTER_MAX < 1
|
|
#define WDL_BESSEL_FILTER_MAX 1
|
|
#elif !defined(WDL_BESSEL_FILTER_MAX) || WDL_BESSEL_FILTER_MAX > 10
|
|
#define WDL_BESSEL_FILTER_MAX 10
|
|
#endif
|
|
#endif
|
|
|
|
|
|
class WDL_BesselFilterCoeffs
|
|
{
|
|
friend class WDL_BesselFilterStage;
|
|
|
|
public:
|
|
inline WDL_BesselFilterCoeffs() {}
|
|
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
// alpha = cornerFreq / (oversampling * sampleRate)
|
|
inline WDL_BesselFilterCoeffs(const double alpha) { Calc(alpha); }
|
|
|
|
void Calc(double alpha)
|
|
{
|
|
const int order = WDL_BESSEL_FILTER_ORDER;
|
|
|
|
#else
|
|
inline WDL_BesselFilterCoeffs(const double alpha, const int order) { Calc(alpha, order); }
|
|
|
|
void Calc(double alpha, const int order)
|
|
{
|
|
assert(order >= 1 && order <= WDL_BESSEL_FILTER_MAX);
|
|
mOrder = order;
|
|
#endif
|
|
|
|
assert(alpha >= 1e-37 && alpha < 0.5);
|
|
alpha *= 6.283185307179586476; // 2.*M_PI
|
|
|
|
// compute S-plane poles for prototype LP filter
|
|
// transform prototype into appropriate filter type (lp)
|
|
// given S-plane poles & zeros, compute Z-plane poles & zeros, by matched z-transform
|
|
complex zplane[WDL_BESSEL_FILTER_MAX];
|
|
int p = (order*order) / 4;
|
|
int n = 0;
|
|
if (order & 1) zplane[n++] = exp(multiply(alpha, mPoles[p++]));
|
|
for (int i = 0; i < order / 2; ++i)
|
|
{
|
|
zplane[n++] = exp(multiply(alpha, mPoles[p]));
|
|
zplane[n++] = exp(multiply(alpha, conjugate(mPoles[p++])));
|
|
}
|
|
|
|
// compute product of poles or zeros as a polynomial of z
|
|
complex coeffs[WDL_BESSEL_FILTER_MAX + 1];
|
|
coeffs[0] = 1.;
|
|
for (int i = 1; i <= order; ++i) coeffs[i] = 0.;
|
|
|
|
for (int i = 0; i < order; ++i)
|
|
{
|
|
// multiply factor (z-w) into coeffs
|
|
complex w = minus(zplane[i]);
|
|
for (int i = order; i >= 1; --i) coeffs[i] = add(multiply(w, coeffs[i]), coeffs[i - 1]);
|
|
coeffs[0] = multiply(w, coeffs[0]);
|
|
}
|
|
#ifndef NDEBUG
|
|
// check computed coeffs of z^k are all real
|
|
for (int n = 0; n <= order; ++n)
|
|
{
|
|
// mkfilter: coeff of z^n is not real; poles/zeros are not complex conjugates
|
|
assert(fabs(coeffs[n].im) <= 1e-10);
|
|
}
|
|
#endif
|
|
|
|
// given Z-plane poles [& zeros], compute [top &] bot polynomials in Z, and then recurrence relation
|
|
complex gain = 0.;
|
|
for (int i = order; i >= 0; --i) gain = add(gain, coeffs[i]);
|
|
gain = inverse(gain);
|
|
mCoeffs[0] = 1./hypot(gain.im, gain.re);
|
|
for (int i = 1, j = order - 1; i <= order; ++i, --j) mCoeffs[i] = -(coeffs[j].re / coeffs[order].re);
|
|
}
|
|
|
|
inline int Order() const
|
|
{
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
return WDL_BESSEL_FILTER_ORDER;
|
|
#else
|
|
return mOrder;
|
|
#endif
|
|
}
|
|
|
|
inline const double* Coeffs() const { return mCoeffs; }
|
|
inline double Gain() const { return mCoeffs[0]; }
|
|
|
|
protected:
|
|
double mCoeffs[WDL_BESSEL_FILTER_MAX + 1];
|
|
|
|
#ifndef WDL_BESSEL_FILTER_ORDER
|
|
int mOrder;
|
|
#endif
|
|
|
|
// Minimalistic complex number implementation
|
|
|
|
struct complex
|
|
{
|
|
double re, im;
|
|
complex() {}
|
|
complex(const double r, const double j = 0.): re(r), im(j) {}
|
|
};
|
|
|
|
// z = z1 + z2
|
|
inline complex add(const complex z1, const complex z2) const
|
|
{
|
|
return complex(z1.re + z2.re, z1.im + z2.im);
|
|
}
|
|
|
|
// z = r * z
|
|
inline complex multiply(const double r, const complex z) const
|
|
{
|
|
return complex(r * z.re, r * z.im);
|
|
}
|
|
|
|
// z = z1 * z2
|
|
inline complex multiply(const complex z1, const complex z2) const
|
|
{
|
|
return complex(z1.re * z2.re - z1.im * z2.im, z1.re * z2.im + z1.im * z2.re);
|
|
}
|
|
|
|
// z = -z
|
|
inline complex minus(const complex z) const
|
|
{
|
|
return complex(-z.re, -z.im);
|
|
}
|
|
|
|
// z = conjugate(z)
|
|
inline complex conjugate(const complex z) const
|
|
{
|
|
return complex(z.re, -z.im);
|
|
}
|
|
|
|
// z = 1/z
|
|
inline complex inverse(const complex z) const
|
|
{
|
|
const double r = z.re*z.re + z.im*z.im;
|
|
return complex(z.re / r, -z.im / r);
|
|
}
|
|
|
|
// z = exp(z)
|
|
inline complex exp(const complex z) const
|
|
{
|
|
const double r = ::exp(z.re);
|
|
return complex(r * cos(z.im), r * sin(z.im));
|
|
}
|
|
|
|
// Precalculated Bessel poles
|
|
static const complex mPoles[10 * 3];
|
|
} WDL_FIXALIGN;
|
|
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
#define WDL_BESSEL_FILTER_OUTPUT(n) (coeffs[WDL_BESSEL_FILTER_ORDER - n + 1] * mOutput[WDL_BESSEL_FILTER_ORDER - n + 1])
|
|
#endif
|
|
|
|
class WDL_BesselFilterStage
|
|
{
|
|
public:
|
|
inline WDL_BesselFilterStage() {}
|
|
inline WDL_BesselFilterStage(const double value) { Reset(value); }
|
|
|
|
inline void Reset() { memset(mOutput, 0, sizeof(mOutput)); }
|
|
|
|
inline void Reset(const double value)
|
|
{
|
|
for (int i = 0; i < WDL_BESSEL_FILTER_MAX; ++i) mOutput[i] = value;
|
|
}
|
|
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
|
|
inline void Process(const double input, const double* const coeffs)
|
|
{
|
|
#if WDL_BESSEL_FILTER_ORDER >= 10
|
|
mOutput[10] = mOutput[9];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 9
|
|
mOutput[9] = mOutput[8];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 8
|
|
mOutput[8] = mOutput[7];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 7
|
|
mOutput[7] = mOutput[6];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 6
|
|
mOutput[6] = mOutput[5];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 5
|
|
mOutput[5] = mOutput[4];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 4
|
|
mOutput[4] = mOutput[3];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 3
|
|
mOutput[3] = mOutput[2];
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 2
|
|
mOutput[2] = mOutput[1];
|
|
#endif
|
|
mOutput[1] = mOutput[0];
|
|
|
|
mOutput[0] = coeffs[0] * input
|
|
#if WDL_BESSEL_FILTER_ORDER >= 10
|
|
+ WDL_BESSEL_FILTER_OUTPUT(10)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 9
|
|
+ WDL_BESSEL_FILTER_OUTPUT(9)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 8
|
|
+ WDL_BESSEL_FILTER_OUTPUT(8)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 7
|
|
+ WDL_BESSEL_FILTER_OUTPUT(7)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 6
|
|
+ WDL_BESSEL_FILTER_OUTPUT(6)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 5
|
|
+ WDL_BESSEL_FILTER_OUTPUT(5)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 4
|
|
+ WDL_BESSEL_FILTER_OUTPUT(4)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 3
|
|
+ WDL_BESSEL_FILTER_OUTPUT(3)
|
|
#endif
|
|
#if WDL_BESSEL_FILTER_ORDER >= 2
|
|
+ WDL_BESSEL_FILTER_OUTPUT(2)
|
|
#endif
|
|
+ WDL_BESSEL_FILTER_OUTPUT(1);
|
|
WDL_BESSEL_FIX_DENORMAL(&mOutput[0]);
|
|
}
|
|
|
|
inline void Process(const double input, const WDL_BesselFilterCoeffs* const bessel) { Process(input, bessel->mCoeffs); }
|
|
|
|
#else // #elif !defined(WDL_BESSEL_FILTER_ORDER)
|
|
|
|
inline void Process(const double input, const double* const coeffs, const int order)
|
|
{
|
|
double output = coeffs[0] * input;
|
|
for (int i = order; i > 0; --i)
|
|
{
|
|
mOutput[i] = mOutput[i - 1];
|
|
output += coeffs[i] * mOutput[i];
|
|
WDL_BESSEL_FIX_DENORMAL(&output);
|
|
}
|
|
mOutput[0] = output;
|
|
}
|
|
|
|
inline void Process(const double input, const WDL_BesselFilterCoeffs* const bessel) { Process(input, bessel->mCoeffs, bessel->mOrder); }
|
|
|
|
#endif
|
|
|
|
inline double Output() const { return mOutput[0]; }
|
|
|
|
protected:
|
|
double mOutput[WDL_BESSEL_FILTER_MAX + 1];
|
|
} WDL_FIXALIGN;
|
|
|
|
|
|
class WDL_BesselFilter: public WDL_BesselFilterCoeffs, public WDL_BesselFilterStage
|
|
{
|
|
public:
|
|
inline WDL_BesselFilter() {}
|
|
|
|
#ifdef WDL_BESSEL_FILTER_ORDER
|
|
inline WDL_BesselFilter(const double alpha) { Calc(alpha); }
|
|
inline void Process(const double input) { WDL_BesselFilterStage::Process(input, mCoeffs); }
|
|
#else
|
|
inline WDL_BesselFilter(const double alpha, const int order) { Calc(alpha, order); }
|
|
inline void Process(const double input) { WDL_BesselFilterStage::Process(input, mCoeffs, mOrder); }
|
|
#endif
|
|
} WDL_FIXALIGN;
|
|
|
|
|
|
#endif // _BESSELFILTER_H_
|