orastron/brickworks
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 9 Wiki Activity

improved bw_src and added proper doc + fixed bw_src_int coeffs

Browse Source
This commit is contained in:
Stefano D'Angelo 2023-04-04 15:13:26 +02:00
parent 67fcbf8b7d
commit c05b6d84db
2 changed files with 40 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
include/bw_src.h
View File

@ -21,7 +21,7 @@
/*!
* module_type {{{ dsp }}}
* version {{{ 0.4.0 }}}
* requires {{{ bw_config bw_common }}}
* requires {{{ bw_config bw_common bw_math }}}
* description {{{
* Aribtrary-ratio IIR sample rate converter.
* }}}
@ -62,20 +62,30 @@ typedef struct _bw_src_state bw_src_state;
* ```>>> */
static inline void bw_src_init(bw_src_coeffs *BW_RESTRICT coeffs, float ratio);
/*! <<<```
* Initializes input parameter values in `coeffs`.
* XXX
* Initializes `coeffs` using the given resampling `ratio`.
*
* `ratio` must be positive and determines the sample rate of the output
* signal, which will be equal to `ratio` times the sample rate of the input
* signal.
*
* #### bw_src_reset_state()
* ```>>> */
static inline void bw_src_reset_state(const bw_src_coeffs *BW_RESTRICT coeffs, bw_src_state *BW_RESTRICT state, float x0);
/*! <<<```
* XXX
* Resets the given `state` to its initial values using the given `coeffs`
* and the quiescent/initial input value `x0`.
*
* #### bw_src_process()
* ```>>> */
static inline void bw_src_process(const bw_src_coeffs *BW_RESTRICT coeffs, bw_src_state *BW_RESTRICT state, const float *x, float *y, int *n_in_samples, int *n_out_samples);
/*! <<<```
* XXX
* Processes at most the first `n_in_samples` of the input buffer `x` and
* fills the output buffer `y` with at most `n_out_samples` using `coeffs`,
* while using and updating `state`.
*
* After the call `n_in_samples` and `n_out_samples` will contain the actual
* number of consumed input samples and generated output samples,
* respectively.
* }}} */
/*** Implementation ***/
@ -83,6 +93,8 @@ static inline void bw_src_process(const bw_src_coeffs *BW_RESTRICT coeffs, bw_sr
/* WARNING: This part of the file is not part of the public API. Its content may
* change at any time in future versions. Please, do not use it directly. */
#include <bw_math.h>
struct _bw_src_coeffs {
float k;
float a1;
@ -92,8 +104,6 @@ struct _bw_src_coeffs {
float b0;
float b1;
float b2;
float b3;
float b4;
};
struct _bw_src_state {
@ -109,17 +119,18 @@ struct _bw_src_state {
static inline void bw_src_init(bw_src_coeffs *BW_RESTRICT coeffs, float ratio) {
coeffs->k = ratio >= 1.f ? 1.f / ratio : -1.f / ratio;
// TODO: better filter, perhaps optimzied coefficients
ratio = ratio >= 1.f ? ratio : -coeffs->k;
coeffs->a1 = (12.56637061435917f - 8.f * ratio) / (2.f * ratio + 3.141592653589793f);
coeffs->a2 = (ratio * (24.f * ratio - 75.39822368615503f) + 59.21762640653615f) / (ratio * (4.f * ratio + 12.56637061435917f) + 9.869604401089358f);
coeffs->a3 = (ratio * ((150.7964473723101f - 32.f * ratio) * ratio - 236.8705056261446f) + 124.0251067211993f) / (ratio * (ratio * (8.f * ratio + 37.69911184307752f) + 59.21762640653615f) + 31.00627668029982f);
coeffs->a4 = (ratio * (ratio * (ratio * (16.f * ratio - 100.5309649148734f) + 236.8705056261446f) - 248.0502134423985f) + 97.40909103400242f) / (ratio * (ratio * (ratio * (16.f * ratio + 100.5309649148734f) + 236.8705056261446f) + 248.0502134423985f) + 97.40909103400242f);
coeffs->b0 = 97.40909103400242f / (ratio * (ratio * (ratio * (16.f * ratio + 100.5309649148734f) + 236.8705056261446f) + 248.0502134423985f) + 97.40909103400242f);
// 4th-degree Butterworth with cutoff at ratio * Nyquist, using bilinear transform w/ prewarping
const float fc = bw_minf((float)(ratio >= 1.f ? 1.f / ratio : ratio), 0.9f);
const float T = bw_tanf_3(1.570796326794896f / (float)fc);
const float T2 = T * T;
const float k = 1.f / (T * (T * (T * (T + 2.613125929752753f) + 3.414213562373095f) + 2.613125929752753f) + 1.f);
coeffs->b0 = T2 * T2;
coeffs->b1 = 4.f * coeffs->b0;
coeffs->b2 = 6.f * coeffs->b0;
coeffs->b3 = coeffs->b1;
coeffs->b4 = coeffs->b0;
coeffs->a1 = k * (T * (T2 * (4.f * T + 5.226251859505504f) - 5.226251859505504f) - 4.f);
coeffs->a2 = k * (T2 * (6.f * T2 - 6.82842712474619f) + 6.f);
coeffs->a3 = k * (T * (T2 * (4.f * T - 5.226251859505504f) + 5.226251859505504f) - 4.f);
coeffs->a4 = k * (T * (T * (T * (T - 2.613125929752753f) + 3.414213562373095f) - 2.613125929752753f) + 1.f);
}
static inline void bw_src_reset_state(const bw_src_coeffs *BW_RESTRICT coeffs, bw_src_state *BW_RESTRICT state, float x0) {
@ -131,8 +142,8 @@ static inline void bw_src_reset_state(const bw_src_coeffs *BW_RESTRICT coeffs, b
state->z4 = state->z3;
} else {
// TDF-II
state->z4 = (coeffs->b4 - coeffs->a4) * x0;
state->z3 = (coeffs->b3 - coeffs->a3) * x0 + state->z4;
state->z4 = (coeffs->b0 - coeffs->a4) * x0;
state->z3 = (coeffs->b1 - coeffs->a3) * x0 + state->z4;
state->z2 = (coeffs->b2 - coeffs->a2) * x0 + state->z3;
state->z1 = (coeffs->b1 - coeffs->a1) * x0 + state->z2;
}
@ -149,7 +160,7 @@ static inline void bw_src_process(const bw_src_coeffs *BW_RESTRICT coeffs, bw_sr
while (i < *n_in_samples && j < *n_out_samples) {
// DF-II
const float z0 = x[i] - coeffs->a1 * state->z1 - coeffs->a2 * state->z2 - coeffs->a3 * state->z3 - coeffs->a4 * state->z4;
const float o = coeffs->b0 * z0 + coeffs->b1 * state->z1 + coeffs->b2 * state->z2 + coeffs->b3 * state->z3 + coeffs->b4 * state->z4;
const float o = coeffs->b0 * (z0 + state->z4) + coeffs->b1 * (state->z1 + state->z3) + coeffs->b2 * state->z2;
if (state->i >= 0.f) {
// 3rd degree Lagrange interpolation + Horner's rule
const float k1 = state->xz1 - state->xz2;
@ -186,11 +197,14 @@ static inline void bw_src_process(const bw_src_coeffs *BW_RESTRICT coeffs, bw_sr
const float c = 0.5f * (k1 + k2);
const float o = state->xz3 + state->i * (a + state->i * (b + state->i * c));
// TDF-II
y[j] = coeffs->b0 * o + state->z1;
state->z1 = coeffs->b1 * o - coeffs->a1 * y[j] + state->z2;
const float v0 = coeffs->b0 * o;
const float v1 = coeffs->b1 * o;
const float v2 = coeffs->b2 * o;
y[j] = v0 + state->z1;
state->z1 = v1 * o - coeffs->a1 * y[j] + state->z2;
state->z2 = coeffs->b2 * o - coeffs->a2 * y[j] + state->z3;
state->z3 = coeffs->b3 * o - coeffs->a3 * y[j] + state->z4;
state->z4 = coeffs->b4 * o - coeffs->a4 * y[j];
state->z3 = v1 * o - coeffs->a3 * y[j] + state->z4;
state->z4 = v0 - coeffs->a4 * y[j];
state->i += coeffs->k;
j++;
}

6
include/bw_src_int.h
View File

@ -73,8 +73,8 @@ static inline void bw_src_int_init(bw_src_int_coeffs *BW_RESTRICT coeffs, int ra
* If `ratio` is positive, then the sample rate of the output signal will be
* `ratio` times the sample rate of the input signal, otherwise, if it is
* negative, then the sample rate of the output signal will be equal to the
* sample rate of the input signal divided by `-ratio`. `ratio` must not be
* `0`.
* sample rate of the input signal divided by `-ratio`. `ratio` *MUST NOT* be
* in [`-1`, `1`].
*
* #### bw_src_int_reset_state()
* ```>>> */
@ -127,7 +127,7 @@ static inline void bw_src_int_init(bw_src_int_coeffs *BW_RESTRICT coeffs, int ra
coeffs->ratio = ratio;
// 4th-degree Butterworth with cutoff at ratio * Nyquist, using bilinear transform w/ prewarping
const float fc = (float)(ratio >= 0 ? ratio : -ratio);
const float T = bw_tanf_3(0.785398163397448f / (float)fc);
const float T = bw_tanf_3(1.570796326794896f / (float)fc);
const float T2 = T * T;
const float k = 1.f / (T * (T * (T * (T + 2.613125929752753f) + 3.414213562373095f) + 2.613125929752753f) + 1.f);
coeffs->b0 = T2 * T2;