add oversampler
This commit is contained in:
412
oversampling/WDL/besselfilter.h
Normal file
412
oversampling/WDL/besselfilter.h
Normal file
@@ -0,0 +1,412 @@
|
||||
/*
|
||||
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_
|
||||
Reference in New Issue
Block a user