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

improved bw_src_int and added proper doc

Browse Source
This commit is contained in:
Stefano D'Angelo 2023-04-04 14:37:38 +02:00
parent 389c6266b5
commit 67fcbf8b7d
1 changed files with 35 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
include/bw_src_int.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 {{{
* Integer-ratio IIR sample rate converter.
*
@ -68,20 +68,33 @@ typedef struct _bw_src_int_state bw_src_int_state;
* ```>>> */
static inline void bw_src_int_init(bw_src_int_coeffs *BW_RESTRICT coeffs, int ratio);
/*! <<<```
* Initializes input parameter values in `coeffs`.
* XXX
* Initializes `coeffs` using the given resampling `ratio`.
*
* 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`.
*
* #### bw_src_int_reset_state()
* ```>>> */
static inline void bw_src_int_reset_state(const bw_src_int_coeffs *BW_RESTRICT coeffs, bw_src_int_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_int_process()
* ```>>> */
static inline int bw_src_int_process(const bw_src_int_coeffs *BW_RESTRICT coeffs, bw_src_int_state *BW_RESTRICT state, const float *x, float *y, int n_in_samples);
/*! <<<```
* XXX
* Processes the first `n_in_samples` of the input buffer `x` and fills the
* output buffer `y` using `coeffs`, while using and updating `state`.
*
* The number of generated output samples will be `ratio` times
* `n_in_samples` if `ratio` is positive, otherwise at most `n_in_samples`
* divided by `-ratio` and then rounded towards positive infinity.
*
* Returns the number of generated output samples.
* }}} */
/*** Implementation ***/
@ -89,6 +102,8 @@ static inline int bw_src_int_process(const bw_src_int_coeffs *BW_RESTRICT coeffs
/* 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_int_coeffs {
int ratio;
float a1;
@ -98,8 +113,6 @@ struct _bw_src_int_coeffs {
float b0;
float b1;
float b2;
float b3;
float b4;
};
struct _bw_src_int_state {
@ -112,17 +125,18 @@ struct _bw_src_int_state {
static inline void bw_src_int_init(bw_src_int_coeffs *BW_RESTRICT coeffs, int ratio) {
coeffs->ratio = ratio;
// TODO: better filter, perhaps optimzied coefficients
ratio = ratio >= 0 ? ratio : -ratio;
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 = (float)(ratio >= 0 ? ratio : -ratio);
const float T = bw_tanf_3(0.785398163397448f / (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_int_reset_state(const bw_src_int_coeffs *BW_RESTRICT coeffs, bw_src_int_state *BW_RESTRICT state, float x0) {
@ -135,8 +149,8 @@ static inline void bw_src_int_reset_state(const bw_src_int_coeffs *BW_RESTRICT c
state->i = 0;
} 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;
}
@ -150,7 +164,7 @@ static inline int bw_src_int_process(const bw_src_int_coeffs *BW_RESTRICT coeffs
const float z0 = x[i] - coeffs->a1 * state->z1 - coeffs->a2 * state->z2 - coeffs->a3 * state->z3 - coeffs->a4 * state->z4;
if (!state->i) {
state->i = -coeffs->ratio;
y[n] = coeffs->b0 * z0 + coeffs->b1 * state->z1 + coeffs->b2 * state->z2 + coeffs->b3 * state->z3 + coeffs->b4 * state->z4;
y[n] = coeffs->b0 * (z0 + state->z4) + coeffs->b1 * (state->z1 + state->z3) + coeffs->b2 * state->z2;
n++;
}
state->i--;
@ -165,14 +179,12 @@ static inline int bw_src_int_process(const bw_src_int_coeffs *BW_RESTRICT coeffs
const float v0 = coeffs->b0 * x[i];
const float v1 = coeffs->b1 * x[i];
const float v2 = coeffs->b2 * x[i];
const float v3 = coeffs->b3 * x[i];
const float v4 = coeffs->b4 * x[i];
for (int j = 0; j < coeffs->ratio; j++) {
y[n] = v0 + state->z1;
state->z1 = v1 - coeffs->a1 * y[n] + state->z2;
state->z2 = v2 - coeffs->a2 * y[n] + state->z3;
state->z3 = v3 - coeffs->a3 * y[n] + state->z4;
state->z4 = v4 - coeffs->a4 * y[n];
state->z3 = v1 - coeffs->a3 * y[n] + state->z4;
state->z4 = v0 - coeffs->a4 * y[n];
n++;
}
}