From 67fcbf8b7daf23a909e6338fbdc164546eb8fce9 Mon Sep 17 00:00:00 2001 From: Stefano D'Angelo Date: Tue, 4 Apr 2023 14:37:38 +0200 Subject: [PATCH] improved bw_src_int and added proper doc --- include/bw_src_int.h | 58 ++++++++++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 23 deletions(-) diff --git a/include/bw_src_int.h b/include/bw_src_int.h index c9c8a4f..3c2a16c 100644 --- a/include/bw_src_int.h +++ b/include/bw_src_int.h @@ -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 + 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++; } }