commit 87f9c013eae4c5993b7f40b5ae98e70c6bc29a89
parent 526c302105feaa6da2c9c2a4d3539364ce7bbdae
Author: Matt Demanett <matt@demanett.net>
Date: Tue, 20 Mar 2018 17:09:53 -0400
Experiemnt with multi-pole Butterworth/Chebyshev LP/HP filter.
Diffstat:
6 files changed, 249 insertions(+), 42 deletions(-)
diff --git a/src/Test.hpp b/src/Test.hpp
@@ -4,7 +4,7 @@
extern Model* modelTest;
-// #define LPF 1
+#define LPF 1
// #define LPFNOISE 1
// #define SINE 1
// #define SQUARE 1
@@ -15,7 +15,7 @@ extern Model* modelTest;
// #define OVERSAMPLING 1
// #define OVERSAMPLED_BL 1
// #define FM 1
-#define PM 1
+// #define PM 1
// #define FEEDBACK_PM 1
// #define EG 1
// #define TABLES 1
diff --git a/src/Test2.cpp b/src/Test2.cpp
@@ -26,6 +26,32 @@ void Test2::step() {
in = inputs[IN_INPUT].value;
}
outputs[OUT_OUTPUT].value = _complexBiquad.next(in);
+
+#elif MULTIPOLE
+ ++_steps;
+ if (_steps >= maxSteps) {
+ _steps = 0;
+
+ _filter.setParams(
+ params[PARAM2A_PARAM].value <= 0.5f ? MultipoleFilter::LP_TYPE : MultipoleFilter::HP_TYPE,
+ 2 * clamp((int)(params[PARAM1B_PARAM].value * (MultipoleFilter::maxPoles / 2)), 1, MultipoleFilter::maxPoles / 2),
+ engineGetSampleRate(),
+ params[PARAM1A_PARAM].value * engineGetSampleRate() / 2.0f,
+ params[PARAM2B_PARAM].value * MultipoleFilter::maxRipple
+ );
+ // _filter.setParams(
+ // MultipoleFilter::HP_TYPE,
+ // 4,
+ // engineGetSampleRate(),
+ // 0.1f * engineGetSampleRate(),
+ // 0.1f
+ // );
+ }
+ float in = 0.0f;
+ if (inputs[IN_INPUT].active) {
+ in = inputs[IN_INPUT].value;
+ }
+ outputs[OUT_OUTPUT].value = _filter.next(in);
#endif
}
diff --git a/src/Test2.hpp b/src/Test2.hpp
@@ -4,10 +4,13 @@
extern Model* modelTest2;
-#define COMPLEX_BIQUAD 1
+// #define COMPLEX_BIQUAD 1
+#define MULTIPOLE 1
#ifdef COMPLEX_BIQUAD
#include "dsp/filter.hpp"
+#elif MULTIPOLE
+#include "dsp/filter.hpp"
#elif
#error what
#endif
@@ -49,6 +52,10 @@ struct Test2 : Module {
#ifdef COMPLEX_BIQUAD
ComplexBiquadFilter _complexBiquad;
+#elif MULTIPOLE
+ MultipoleFilter _filter;
+ const int maxSteps = 10;
+ int _steps = maxSteps;
#endif
Test2() : Module(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS, NUM_LIGHTS) {
diff --git a/src/dsp/buffer.hpp b/src/dsp/buffer.hpp
@@ -121,5 +121,40 @@ struct AveragingBuffer {
}
};
+template<typename T>
+struct HistoryBuffer {
+ int _size, _i;
+ T* _buf;
+
+ HistoryBuffer(int size, T initialValue)
+ : _size(size)
+ , _i(size)
+ {
+ assert(size > 0);
+ _buf = new T[size];
+ std::fill(_buf, _buf + size, initialValue);
+ }
+ ~HistoryBuffer() {
+ delete[] _buf;
+ }
+
+ inline void push(T s) {
+ ++_i;
+ if (_i >= _size) {
+ _i = 0;
+ }
+ _buf[_i] = s;
+ }
+
+ inline T value(int i) {
+ assert(i <= 0 && -i <= _size);
+ int j = _i - i;
+ if (j < 0) {
+ j += _size;
+ }
+ return _buf[j];
+ }
+};
+
} // namespace dsp
} // namespace bogaudio
diff --git a/src/dsp/filter.cpp b/src/dsp/filter.cpp
@@ -6,34 +6,6 @@
using namespace bogaudio::dsp;
-// See: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
-
-void BiquadFilter::setParams(float b0, float b1, float b2, float a0, float a1, float a2) {
- _b0 = b0 / a0;
- _b1 = b1 / a0;
- _b2 = b2 / a0;
- _a1 = a1 / a0;
- _a2 = a2 / a0;
- // printf("Biquad set params: b0=%f b1=%f b2=%f a1=%f a2=%f\n", _b0, _b1, _b2, _a1, _a2);
-}
-
-float BiquadFilter::next(float sample) {
- _x[2] = _x[1];
- _x[1] = _x[0];
- _x[0] = sample;
-
- _y[2] = _y[1];
- _y[1] = _y[0];
- _y[0] = _b0 * _x[0];
- _y[0] += _b1 * _x[1];
- _y[0] += _b2 * _x[2];
- _y[0] -= _a1 * _y[1];
- _y[0] -= _a2 * _y[2];
-
- return _y[0];
-}
-
-
void ComplexBiquadFilter::setComplexParams(
float gain,
float zeroRadius,
@@ -74,6 +46,7 @@ void ComplexBiquadFilter::updateParams() {
}
+// See: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
void LowPassFilter::setParams(float sampleRate, float cutoff, float q) {
if (_sampleRate == sampleRate && _cutoff == cutoff && _q == q) {
return;
@@ -96,3 +69,109 @@ void LowPassFilter::setParams(float sampleRate, float cutoff, float q) {
1.0 - alpha
);
}
+
+
+// Adapted from Smith 1997, "The Scientist and Engineer's Guide to DSP"
+void MultipoleFilter::setParams(
+ Type type,
+ int poles,
+ float sampleRate,
+ float cutoff,
+ float ripple
+) {
+ if (
+ _type == type &&
+ _poles == poles &&
+ _sampleRate == sampleRate &&
+ _cutoff == cutoff &&
+ _ripple == ripple
+ ) {
+ return;
+ }
+ assert(poles >= 1 && poles <= maxPoles);
+ assert(poles % 2 == 0); // relax this later?
+ assert(sampleRate >= 0.0f);
+ assert(cutoff >= 0.0f && cutoff <= sampleRate / 2.0f);
+ assert(ripple >= 0.0f && ripple <= maxRipple);
+ _type = type;
+ _poles = poles;
+ _sampleRate = sampleRate;
+ _cutoff = cutoff;
+ _ripple = ripple;
+
+ for (int p = 0, pn = _poles / 2; p < pn; ++p) {
+ double a0 = 0.0;
+ double a1 = 0.0;
+ double a2 = 0.0;
+ double b1 = 0.0;
+ double b2 = 0.0;
+ {
+ double angle = M_PI/(_poles*2.0) + p*M_PI/_poles;
+ double rp = -cos(angle);
+ double ip = sin(angle);
+
+ if (_ripple > 0.01f) {
+ double es = sqrt(pow(1.0 / (1.0 - _ripple), 2.0) - 1.0);
+ double esi = 1.0 / es;
+ double esis = esi * esi;
+ double polesi = 1.0 / (float)_poles;
+ double vx = polesi * log(esi + sqrt(esis + 1.0));
+ double kx = polesi * log(esi + sqrt(esis - 1.0));
+ kx = (exp(kx) + exp(-kx)) / 2.0;
+ rp *= ((exp(vx) - exp(-vx)) / 2.0) / kx;
+ ip *= ((exp(vx) + exp(-vx)) / 2.0) / kx;
+ // printf("\n\n\ncheb p=%d pn=%d rip=%f rp=%f ip=%f es=%f vx=%f kx=%f\n", p, pn, _ripple, rp, ip, es, vx, kx);
+ }
+
+ const double t = 2.0 * tan(0.5);
+ const double ts = t * t;
+ // printf("\n\nco=%f sr=%f fc=%f\n", _cutoff, _sampleRate, _cutoff / _sampleRate);
+ double m = rp*rp + ip*ip;
+ double mts = m * ts;
+ double d = 4.0 - 4.0*rp*t + mts;
+ double x0 = ts / d;
+ double x1 = 2.0 * x0;
+ double x2 = x0;
+ double y1 = (8.0 - 2.0*mts) / d;
+ double y2 = (-4.0 - 4.0*rp*t - mts) / d;
+ // printf("vs p=%d t=%f w=%f m=%f d=%f x0=%f x1=%f x2=%f y1=%f y2=%f\n", p, t, w, m, d, x0, x1, x2, y1, y2);
+
+ // FIXME: optimization: everything above here only depends on type, poles and ripple.
+
+ double w = 2.0 * M_PI * (_cutoff / _sampleRate);
+ double w2 = w / 2.0;
+ double k = 0.0;
+ switch (_type) {
+ case LP_TYPE: {
+ k = sin(0.5 - w2) / sin(0.5 + w2);
+ break;
+ }
+ case HP_TYPE: {
+ k = -cos(w2 + 0.5) / cos(w2 - 0.5);
+ break;
+ }
+ }
+ double ks = k * k;
+ d = 1.0 + y1*k - y2*ks;
+ a0 = (x0 - x1*k + x2*ks) / d;
+ a1 = (-2.0*x0*k + x1 + x1*ks - 2.0*x2*k) / d;
+ a2 = (x0*ks - x1*k + x2) / d;
+ b1 = (2.0*k + y1 + y1*ks - 2.0*y2*k) / d;
+ b2 = (-ks - y1*k + y2) / d;
+ if (_type == HP_TYPE) {
+ a1 = -a1;
+ b1 = -b1;
+ }
+
+ // printf("pole %d: rp=%f ip=%f a0=%f a1=%f a2=%f b1=%f b2=%f\n\n\n", p, rp, ip, a0, a1, a2, b1, b2);
+ _biquads[p].setParams(a0, a1, a2, 1.0, -b1, -b2);
+ }
+ }
+}
+
+float MultipoleFilter::next(float sample) {
+ for (int p = 0, pn = _poles / 2; p < pn; ++p) {
+ sample = _biquads[p].next(sample);
+ }
+ return sample;
+}
diff --git a/src/dsp/filter.hpp b/src/dsp/filter.hpp
@@ -2,6 +2,8 @@
#include <math.h>
+#include "buffer.hpp"
+
namespace bogaudio {
namespace dsp {
@@ -12,23 +14,54 @@ struct Filter {
virtual float next(float sample) = 0;
};
+template<typename T>
struct BiquadFilter : Filter {
- float _b0 = 0.0;
- float _b1 = 0.0;
- float _b2 = 0.0;
- float _a1 = 0.0;
- float _a2 = 0.0 ;
+ T _a0 = 0.0;
+ T _a1 = 0.0;
+ T _a2 = 0.0;
+ T _b1 = 0.0;
+ T _b2 = 0.0 ;
- float _x[3] {};
- float _y[3] {};
+ T _x[3] {};
+ T _y[3] {};
BiquadFilter() {}
- void setParams(float b0, float b1, float b2, float a0, float a1, float a2);
- virtual float next(float sample) override;
+ void setParams(T a0, T a1, T a2, T b0, T b1, T b2) {
+ if (b0 == 1.0) {
+ _a0 = a0;
+ _a1 = a1;
+ _a2 = a2;
+ _b1 = b1;
+ _b2 = b2;
+ }
+ else {
+ _a0 = a0 / b0;
+ _a1 = a1 / b0;
+ _a2 = a2 / b0;
+ _b1 = b1 / b0;
+ _b2 = b2 / b0;
+ }
+ }
+
+ virtual float next(float sample) override {
+ _x[2] = _x[1];
+ _x[1] = _x[0];
+ _x[0] = sample;
+
+ _y[2] = _y[1];
+ _y[1] = _y[0];
+ _y[0] = _a0 * _x[0];
+ _y[0] += _a1 * _x[1];
+ _y[0] += _a2 * _x[2];
+ _y[0] -= _b1 * _y[1];
+ _y[0] -= _b2 * _y[2];
+
+ return _y[0];
+ }
};
-struct ComplexBiquadFilter : BiquadFilter {
+struct ComplexBiquadFilter : BiquadFilter<float> {
float _gain = 1.0f;
float _zeroRadius = 1.0f;
float _zeroTheta = M_PI;
@@ -54,7 +87,7 @@ struct LowPassFilter : Filter {
float _cutoff;
float _q;
- BiquadFilter _biquad;
+ BiquadFilter<float> _biquad;
LowPassFilter(float sampleRate, float cutoff, float q)
: _sampleRate(sampleRate)
@@ -69,5 +102,32 @@ struct LowPassFilter : Filter {
}
};
+struct MultipoleFilter : Filter {
+ enum Type {
+ LP_TYPE,
+ HP_TYPE
+ };
+
+ static constexpr int maxPoles = 20;
+ static constexpr float maxRipple = 0.29f;
+ Type _type = LP_TYPE;
+ int _poles = 1;
+ float _sampleRate = 0.0f;
+ float _cutoff = 0.0f;
+ float _ripple = 0.0f;
+ BiquadFilter<double> _biquads[maxPoles / 2];
+
+ MultipoleFilter() {}
+
+ void setParams(
+ Type type,
+ int poles,
+ float sampleRate,
+ float cutoff,
+ float ripple // FIXME: using this with more than two poles creates large gain, need compensation.
+ );
+ virtual float next(float sample) override;
+};
+
} // namespace dsp
} // namespace bogaudio
