VCVRack
/
Befaco
mirror of https://github.com/VCVRack/Befaco.git
1
0
Fork 0
Code Issues 0 Releases 17 Wiki Activity
Browse Source

Switch to biquad for DC 

* only up/downsample if needed
tags/v1.1.0^2
hemmer 4 years ago
parent
a92f2df465
commit
b4241e035d
2 changed files with 150 additions and 382 deletions
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +30
    -23
      src/ChoppingKinky.cpp
  2. +120
    -359
      src/VariableOversampling.hpp

+ 30
- 23
src/ChoppingKinky.cpp View File

@@ -67,19 +67,14 @@ struct ChoppingKinky : Module {
dsp::SchmittTrigger trigger;
bool output_a_to_chopp;
float previous_a = 0.0;
VariableOversampling<> oversampler[NUM_CHANNELS];
int oversamplingIndex = 2;
dsp::RCFilter choppDC;
dsp::BiquadFilter blockDCFilter;
bool blockDC = false;
unsigned int currentFrame = 0;
ChoppingKinky() {
config(NUM_PARAMS, NUM_INPUTS, NUM_OUTPUTS, NUM_LIGHTS);
configParam(FOLD_A_PARAM, 0.f, 2.f, 0.f, "");
@@ -89,14 +84,13 @@ struct ChoppingKinky : Module {
// calculate up/downsampling rates
onSampleRateChange();
}
void onSampleRateChange() override {
// SampleRateConverter needs integer value
int sampleRate = APP->engine->getSampleRate();
choppDC.setCutoff(10.3f / sampleRate);
blockDCFilter.setParameters(dsp::BiquadFilter::HIGHPASS, 10.3f / sampleRate, M_SQRT1_2, 1.0f);
for (int channel_idx = 0; channel_idx < NUM_CHANNELS; channel_idx++) {
oversampler[channel_idx].setOversamplingIndex(oversamplingIndex);
@@ -150,33 +144,46 @@ struct ChoppingKinky : Module {
}
}
bool chopp_is_required = outputs[OUT_CHOPP_OUTPUT].isConnected();
bool a_is_required = outputs[OUT_A_OUTPUT].isConnected() || chopp_is_required;
bool b_is_required = outputs[OUT_B_OUTPUT].isConnected() || chopp_is_required;
in_a = in_a * gain_a;
in_b = in_b * gain_b;
float in_chopp = output_a_to_chopp ? 1.f : 0.f;
oversampler[CHANNEL_A].upsample(in_a);
oversampler[CHANNEL_B].upsample(in_b);
oversampler[CHANNEL_CHOPP].upsample(in_chopp);
if (a_is_required) {
oversampler[CHANNEL_A].upsample(in_a);
}
if (b_is_required) {
oversampler[CHANNEL_B].upsample(in_b);
}
if (chopp_is_required) {
oversampler[CHANNEL_CHOPP].upsample(in_chopp);
}
float* osBufferA = oversampler[CHANNEL_A].getOSBuffer();
float* osBufferB = oversampler[CHANNEL_B].getOSBuffer();
float* osBufferChopp = oversampler[CHANNEL_CHOPP].getOSBuffer();
for (int i = 0; i < oversampler[0].getOversamplingRatio(); i++) {
// TOO SLOW
//break;
osBufferA[i] = foldResponse(4.0f * osBufferA[i], -0.2);
osBufferB[i] = foldResponseSin(4.0f * osBufferB[i]);
osBufferChopp[i] = osBufferChopp[i] * osBufferA[i] + (1.f - osBufferChopp[i]) * osBufferB[i];
if (a_is_required) {
osBufferA[i] = foldResponse(4.0f * osBufferA[i], -0.2);
}
if (b_is_required) {
osBufferB[i] = foldResponseSin(4.0f * osBufferB[i]);
}
if (chopp_is_required) {
osBufferChopp[i] = osBufferChopp[i] * osBufferA[i] + (1.f - osBufferChopp[i]) * osBufferB[i];
}
}
float out_a = oversampler[CHANNEL_A].downsample();
float out_b = oversampler[CHANNEL_B].downsample();
float out_chopp = oversampler[CHANNEL_CHOPP].downsample();
float out_a = a_is_required ? oversampler[CHANNEL_A].downsample() : 0.f;
float out_b = b_is_required ? oversampler[CHANNEL_B].downsample() : 0.f;
float out_chopp = chopp_is_required ? oversampler[CHANNEL_CHOPP].downsample() : 0.f;
if (blockDC) {
choppDC.process(out_chopp);
out_chopp -= choppDC.lowpass();
out_chopp = blockDCFilter.process(out_chopp);
}
previous_a = in_a;


+ 120
- 359
src/VariableOversampling.hpp View File

@@ -2,336 +2,86 @@

#include <rack.hpp>
#include <dsp/common.hpp>
#include <type_traits>

/** Digital IIR filter processor. Using TDF-II structure:
* https://ccrma.stanford.edu/~jos/filters/Transposed_Direct_Forms.html
*/
template <int ORDER, typename T = float>
struct IIRFilter {
/** transfer function numerator coefficients: b_0, b_1, etc.*/
T b[ORDER] = {};

/** transfer function denominator coefficients: a_0, a_1, etc.*/
T a[ORDER] = {};

/** filter state */
T z[ORDER];
/**
High-order filter to be used for anti-aliasing or anti-imaging.
The template parameter N should be 1/2 the desired filter order.

IIRFilter() {
reset();
}
Currently uses an 2*N-th order Butterworth filter.
@TODO: implement Chebyshev, Elliptic filter options.
*/
template<int N>
class AAFilter {
public:
AAFilter() = default;

void reset() {
std::fill(z, &z[ORDER], 0.0f);
}
/** Calculate Q values for a Butterworth filter of a given order */
static std::vector<float> calculateButterQs(int order) {
const int lim = int (order / 2);
std::vector<float> Qs;

void setCoefficients(const T* b, const T* a) {
for (int i = 0; i < ORDER; i++) {
this->b[i] = b[i];
for (int k = 1; k <= lim; ++k) {
auto b = -2.0f * std::cos((2.0f * k + order - 1) * 3.14159 / (2.0f * order));
Qs.push_back(1.0f / b);
}
for (int i = 1; i < ORDER; i++) {
this->a[i] = a[i];
}
}

template <int N = ORDER>
inline typename std::enable_if <N == 2, T>::type
process(T x) noexcept {
T y = z[1] + x * b[0];
z[1] = x * b[1] - y * a[1];
return y;
}

template <int N = ORDER>
inline typename std::enable_if <N == 3, T>::type
process(T x) noexcept {
T y = z[1] + x * b[0];
z[1] = z[2] + x * b[1] - y * a[1];
z[2] = x * b[2] - y * a[2];
return y;
}

template <int N = ORDER>
inline typename std::enable_if <(N > 3), T>::type
process(T x) noexcept {
T y = z[1] + x * b[0];

for (int i = 1; i < ORDER-1; ++i)
z[i] = z[i+1] + x * b[i] - y * a[i];

z[ORDER-1] = x * b[ORDER-1] - y * a[ORDER-1];

return y;
}

/** Computes the complex transfer function $H(s)$ at a particular frequency
s: normalized angular frequency equal to $2 \pi f / f_{sr}$ ($\pi$ is the Nyquist frequency)
*/
std::complex<T> getTransferFunction(T s) {
// Compute sum(a_k z^-k) / sum(b_k z^-k) where z = e^(i s)
std::complex<T> bSum(b[0], 0);
std::complex<T> aSum(1, 0);
for (int i = 1; i < ORDER; i++) {
T p = -i * s;
std::complex<T> z(simd::cos(p), simd::sin(p));
bSum += b[i] * z;
aSum += a[i - 1] * z;
}
return bSum / aSum;
std::reverse(Qs.begin(), Qs.end());
return Qs;
}

T getFrequencyResponse(T f) {
return simd::abs(getTransferFunction(2 * M_PI * f));
/**
* Resets the filter to process at a new sample rate.
*
* @param sampleRate: The base (i.e. pre-oversampling) sample rate of the audio being processed
* @param osRatio: The oversampling ratio at which the filter is being used
*/
void reset(float sampleRate, int osRatio) {
float fc = 0.98f * (sampleRate / 2.0f);
auto Qs = calculateButterQs(2 * N);

for (int i = 0; i < N; ++i)
filters[i].setParameters(dsp::BiquadFilter::Type::LOWPASS, fc / (osRatio * sampleRate), Qs[i], 1.0f);
}

T getFrequencyPhase(T f) {
return simd::arg(getTransferFunction(2 * M_PI * f));
}
};

template <typename T = float>
struct TOnePole : IIRFilter<2, T> {
enum Type {
LOWPASS,
HIGHPASS,
NUM_TYPES
};

TOnePole() {
setParameters(LOWPASS, 0.0f);
}

/** Calculates and sets the one-pole transfer function coefficients.
f: normalized frequency (cutoff frequency / sample rate), must be less than 0.5
*/
void setParameters(Type type, float f) {
switch (type) {
case LOWPASS: {
this->a[1] = -std::exp(-2.f * M_PI * f);
this->b[0] = 1.f + this->a[0];
} break;

case HIGHPASS: {
this->a[1] = std::exp(-2.f * M_PI * (0.5f - f));
this->b[0] = 1.f - this->a[0];
} break;

default: break;
}
}
};
inline float process(float x) noexcept {
for (int i = 0; i < N; ++i)
x = filters[i].process(x);

typedef TOnePole<> OnePole;

template <typename T = float>
struct TBiquadFilter : IIRFilter<3, T> {
enum Type {
LOWPASS,
HIGHPASS,
LOWSHELF,
HIGHSHELF,
BANDPASS,
PEAK,
NOTCH,
NUM_TYPES
};

TBiquadFilter() {
setParameters(LOWPASS, 0.f, 0.f, 1.f);
return x;
}

/** Calculates and sets the biquad transfer function coefficients.
f: normalized frequency (cutoff frequency / sample rate), must be less than 0.5
Q: quality factor
V: gain
*/
void setParameters(Type type, float f, float Q, float V) {
float K = std::tan(M_PI * f);
switch (type) {
case LOWPASS: {
float norm = 1.f / (1.f + K / Q + K * K);
this->b[0] = K * K * norm;
this->b[1] = 2.f * this->b[0];
this->b[2] = this->b[0];
this->a[1] = 2.f * (K * K - 1.f) * norm;
this->a[2] = (1.f - K / Q + K * K) * norm;
} break;

case HIGHPASS: {
float norm = 1.f / (1.f + K / Q + K * K);
this->b[0] = norm;
this->b[1] = -2.f * this->b[0];
this->b[2] = this->b[0];
this->a[1] = 2.f * (K * K - 1.f) * norm;
this->a[2] = (1.f - K / Q + K * K) * norm;

} break;

case LOWSHELF: {
float sqrtV = std::sqrt(V);
if (V >= 1.f) {
float norm = 1.f / (1.f + M_SQRT2 * K + K * K);
this->b[0] = (1.f + M_SQRT2 * sqrtV * K + V * K * K) * norm;
this->b[1] = 2.f * (V * K * K - 1.f) * norm;
this->b[2] = (1.f - M_SQRT2 * sqrtV * K + V * K * K) * norm;
this->a[1] = 2.f * (K * K - 1.f) * norm;
this->a[2] = (1.f - M_SQRT2 * K + K * K) * norm;
}
else {
float norm = 1.f / (1.f + M_SQRT2 / sqrtV * K + K * K / V);
this->b[0] = (1.f + M_SQRT2 * K + K * K) * norm;
this->b[1] = 2.f * (K * K - 1) * norm;
this->b[2] = (1.f - M_SQRT2 * K + K * K) * norm;
this->a[1] = 2.f * (K * K / V - 1.f) * norm;
this->a[2] = (1.f - M_SQRT2 / sqrtV * K + K * K / V) * norm;
}
} break;

case HIGHSHELF: {
float sqrtV = std::sqrt(V);
if (V >= 1.f) {
float norm = 1.f / (1.f + M_SQRT2 * K + K * K);
this->b[0] = (V + M_SQRT2 * sqrtV * K + K * K) * norm;
this->b[1] = 2.f * (K * K - V) * norm;
this->b[2] = (V - M_SQRT2 * sqrtV * K + K * K) * norm;
this->a[1] = 2.f * (K * K - 1.f) * norm;
this->a[2] = (1.f - M_SQRT2 * K + K * K) * norm;
}
else {
float norm = 1.f / (1.f / V + M_SQRT2 / sqrtV * K + K * K);
this->b[0] = (1.f + M_SQRT2 * K + K * K) * norm;
this->b[1] = 2.f * (K * K - 1.f) * norm;
this->b[2] = (1.f - M_SQRT2 * K + K * K) * norm;
this->a[1] = 2.f * (K * K - 1.f / V) * norm;
this->a[2] = (1.f / V - M_SQRT2 / sqrtV * K + K * K) * norm;
}
} break;

case BANDPASS: {
float norm = 1.f / (1.f + K / Q + K * K);
this->b[0] = K / Q * norm;
this->b[1] = 0.f;
this->b[2] = -this->b[0];
this->a[1] = 2.f * (K * K - 1.f) * norm;
this->a[2] = (1.f - K / Q + K * K) * norm;
} break;

case PEAK: {
float c = 1.0f / K;
float phi = c * c;
float Knum = c / Q;
float Kdenom = Knum;

if(V > 1.0f)
Knum *= V;
else
Kdenom /= V;

float norm = phi + Kdenom + 1.0;
this->b[0] = (phi + Knum + 1.0f) / norm;
this->b[1] = 2.0f * (1.0f - phi) / norm;
this->b[2] = (phi - Knum + 1.0f) / norm;
this->a[1] = 2.0f * (1.0f - phi) / norm;
this->a[2] = (phi - Kdenom + 1.0f) / norm;
} break;

case NOTCH: {
float norm = 1.f / (1.f + K / Q + K * K);
this->b[0] = (1.f + K * K) * norm;
this->b[1] = 2.f * (K * K - 1.f) * norm;
this->b[2] = this->b[0];
this->a[1] = this->b[1];
this->a[2] = (1.f - K / Q + K * K) * norm;
} break;

default: break;
}
}
};

typedef TBiquadFilter<> BiquadFilter;


/**
High-order filter to be used for anti-aliasing or anti-imaging.
The template parameter N should be 1/2 the desired filter order.

Currently uses an 2*N-th order Butterworth filter.
@TODO: implement Chebyshev, Elliptic filter options.
*/
template<int N>
class AAFilter
{
public:
AAFilter() = default;

/** Calculate Q values for a Butterworth filter of a given order */
static std::vector<float> calculateButterQs(int order) {
const int lim = int (order / 2);
std::vector<float> Qs;

for(int k = 1; k <= lim; ++k) {
auto b = -2.0f * std::cos((2.0f * k + order - 1) * 3.14159 / (2.0f * order));
Qs.push_back(1.0f / b);
}

std::reverse(Qs.begin(), Qs.end());
return Qs;
}

/**
* Resets the filter to process at a new sample rate.
*
* @param sampleRate: The base (i.e. pre-oversampling) sample rate of the audio being processed
* @param osRatio: The oversampling ratio at which the filter is being used
*/
void reset(float sampleRate, int osRatio) {
float fc = 0.98f * (sampleRate / 2.0f);
auto Qs = calculateButterQs(2*N);

for(int i = 0; i < N; ++i)
filters[i].setParameters(BiquadFilter::Type::LOWPASS, fc / (osRatio * sampleRate), Qs[i], 1.0f);
}
inline float process(float x) noexcept {
for(int i = 0; i < N; ++i)
x = filters[i].process(x);
return x;
}

private:
BiquadFilter filters[N];
dsp::BiquadFilter filters[N];
};



/**
* Base class for oversampling of any order
*/
*/
class BaseOversampling {
public:
BaseOversampling() = default;
virtual ~BaseOversampling() {}
BaseOversampling() = default;
virtual ~BaseOversampling() {}

/** Resets the oversampler for processing at some base sample rate */
virtual void reset(float /*baseSampleRate*/) = 0;
/** Resets the oversampler for processing at some base sample rate */
virtual void reset(float /*baseSampleRate*/) = 0;

/** Upsample a single input sample and update the oversampled buffer */
virtual void upsample(float) noexcept = 0;
/** Upsample a single input sample and update the oversampled buffer */
virtual void upsample(float) noexcept = 0;

/** Output a downsampled output sample from the current oversampled buffer */
virtual float downsample() noexcept = 0;
/** Output a downsampled output sample from the current oversampled buffer */
virtual float downsample() noexcept = 0;

/** Returns a pointer to the oversampled buffer */
virtual float* getOSBuffer() noexcept = 0;
/** Returns a pointer to the oversampled buffer */
virtual float* getOSBuffer() noexcept = 0;
};

/**
/**
Class to implement an oversampled process.
To use, create an object and prepare using `reset()`.

Then use the following code to process samples:
@code
oversample.upsample(x);
@@ -341,51 +91,50 @@ public:
@endcode
*/
template<int ratio, int filtN = 4>
class Oversampling : public BaseOversampling
{
class Oversampling : public BaseOversampling {
public:
Oversampling() = default;
virtual ~Oversampling() {}
Oversampling() = default;
virtual ~Oversampling() {}

void reset(float baseSampleRate) override {
aaFilter.reset(baseSampleRate, ratio);
aiFilter.reset(baseSampleRate, ratio);
std::fill(osBuffer, &osBuffer[ratio], 0.0f);
}
void reset(float baseSampleRate) override {
aaFilter.reset(baseSampleRate, ratio);
aiFilter.reset(baseSampleRate, ratio);
std::fill(osBuffer, &osBuffer[ratio], 0.0f);
}

inline void upsample(float x) noexcept override {
osBuffer[0] = ratio * x;
std::fill(&osBuffer[1], &osBuffer[ratio], 0.0f);
inline void upsample(float x) noexcept override {
osBuffer[0] = ratio * x;
std::fill(&osBuffer[1], &osBuffer[ratio], 0.0f);

for(int k = 0; k < ratio; k++)
osBuffer[k] = aiFilter.process(osBuffer[k]);
}
for (int k = 0; k < ratio; k++)
osBuffer[k] = aiFilter.process(osBuffer[k]);
}

inline float downsample() noexcept override {
float y = 0.0f;
for(int k = 0; k < ratio; k++)
y = aaFilter.process(osBuffer[k]);
inline float downsample() noexcept override {
float y = 0.0f;
for (int k = 0; k < ratio; k++)
y = aaFilter.process(osBuffer[k]);

return y;
}
return y;
}

inline float* getOSBuffer() noexcept override {
return osBuffer;
}
inline float* getOSBuffer() noexcept override {
return osBuffer;
}

float osBuffer[ratio];
float osBuffer[ratio];

private:
AAFilter<filtN> aaFilter; // anti-aliasing filter
AAFilter<filtN> aiFilter; // anti-imaging filter
AAFilter<filtN> aaFilter; // anti-aliasing filter
AAFilter<filtN> aiFilter; // anti-imaging filter
};


/**
/**
Class to implement an oversampled process, with variable
oversampling factor. To use, create an object, set the oversampling
factor using `setOversamplingindex()` and prepare using `reset()`.
Then use the following code to process samples:
@code
oversample.upsample(x);
@@ -398,45 +147,57 @@ private:
template<int filtN = 4>
class VariableOversampling {
public:
VariableOversampling() = default;
VariableOversampling() = default;

/** Prepare the oversampler to process audio at a given sample rate */
void reset(float sampleRate) {
for(auto* os : oss)
os->reset(sampleRate);
}
/** Prepare the oversampler to process audio at a given sample rate */
void reset(float sampleRate) {
for (auto* os : oss)
os->reset(sampleRate);
}

/** Sets the oversampling factor as 2^idx */
void setOversamplingIndex (int newIdx) { osIdx = newIdx; }
/** Sets the oversampling factor as 2^idx */
void setOversamplingIndex(int newIdx) {
osIdx = newIdx;
}

/** Returns the oversampling index */
int getOversamplingIndex() const noexcept {
return osIdx;
}

/** Returns the oversampling index */
int getOversamplingIndex() const noexcept { return osIdx; }
/** Upsample a single input sample and update the oversampled buffer */
inline void upsample(float x) noexcept {
oss[osIdx]->upsample(x);
}

/** Upsample a single input sample and update the oversampled buffer */
inline void upsample(float x) noexcept { oss[osIdx]->upsample(x); }
/** Output a downsampled output sample from the current oversampled buffer */
inline float downsample() noexcept {
return oss[osIdx]->downsample();
}

/** Output a downsampled output sample from the current oversampled buffer */
inline float downsample() noexcept { return oss[osIdx]->downsample(); }
/** Returns a pointer to the oversampled buffer */
inline float* getOSBuffer() noexcept {
return oss[osIdx]->getOSBuffer();
}

/** Returns a pointer to the oversampled buffer */
inline float* getOSBuffer() noexcept { return oss[osIdx]->getOSBuffer(); }
/** Returns the current oversampling factor */
int getOversamplingRatio() const noexcept {
return 1 << osIdx;
}

/** Returns the current oversampling factor */
int getOversamplingRatio() const noexcept { return 1 << osIdx; }

private:
enum {
NumOS = 5, // number of oversampling options
};
int osIdx = 0;
Oversampling<1 << 0, filtN> os0; // 1x
Oversampling<1 << 1, filtN> os1; // 2x
Oversampling<1 << 2, filtN> os2; // 4x
Oversampling<1 << 3, filtN> os3; // 8x
Oversampling<1 << 4, filtN> os4; // 16x
BaseOversampling* oss[NumOS] = { &os0, &os1, &os2, &os3, &os4 };
enum {
NumOS = 5, // number of oversampling options
};
int osIdx = 0;
Oversampling < 1 << 0, filtN > os0; // 1x
Oversampling < 1 << 1, filtN > os1; // 2x
Oversampling < 1 << 2, filtN > os2; // 4x
Oversampling < 1 << 3, filtN > os3; // 8x
Oversampling < 1 << 4, filtN > os4; // 16x
BaseOversampling* oss[NumOS] = { &os0, &os1, &os2, &os3, &os4 };

};

Write Preview
Loading…
Cancel
Save