Polyphase Resampling (rational resampling)

[xerpi]

Since CMSIS-DSP already provides functions for FIR decimation and FIR interpolation, it would be beneficial to include efficient rational resampling functionality using the polyphase resampling algorithm.

For reference, see SciPy's scipy.signal.resample_poly, GNU Radio's Rational Resampler, FutureSDR's PolyphaseResamplingFirKernel, SampleRateConverter, and liquid-dsp's rresamp.

[christophe0606]

Good suggestion. Tagged as "enhancement".

[dafx]

Been wanting this forever! Thanks.

[xerpi]

Some PoC code that seems to work well.
We want a polyphase filterbank with N sub-filters. Each filter has a delay of M samples.
Each sub-filter length will be S = 2 * M. The total filter length will be S * N = 2 * M * N.
In liquid-dsp, for the rational resampler, a filter with an extra coefficient (2 * M * N + 1) is designed, not sure why.

Then we can reorder the coefficients so that each sub-filter is contiguous in memory.
We can even do that in-place using an In-place matrix transposition algorithm:

/* Swaps two q15_t values */
static inline __attribute__((always_inline)) void swap_q15(q15_t *a, q15_t *b)
{
    q15_t tmp = *a;
    *a = *b;
    *b = tmp;
}

/* Rotates elements in range [first, last) by one position to the left */
static void rotate_left(q15_t *arr, size_t first, size_t last)
{
    size_t i;
    q15_t temp = arr[first];

    for (i = first; i < last - 1; i++)
        arr[i] = arr[i + 1];

    arr[last - 1] = temp;
}

/*
 * Source: https://stackoverflow.com/a/8733699
 */
void matrix_transpose_q15(q15_t *matrix, size_t width, size_t height)
{
    size_t x, y, step, count_adjustment, last, first;
    const size_t count = width * height;

    for (x = 0; x < width; ++x) {
        count_adjustment = width - x - 1;
        for (y = 0, step = 1; y < height; ++y, step += count_adjustment) {
            last = count - (y + x * height);
            first = last - step;
            rotate_left(matrix, first, last);
        }
    }
}

/* Reverses the order of elements in each row of the matrix */
void matrix_reverse_rows_q15(q15_t *matrix, size_t width, size_t height)
{
    size_t i, j;

    for (i = 0; i < height; ++i) {
        for (j = 0; j < width / 2; ++j)
            swap_q15(&matrix[i * width + j], &matrix[i * width + (width - 1 - j)]);
    }
}

/*
 * Reorders the matrix for polyphase filtering.
 * Each sub-filter is inverted so that it can be applied using a dot product (arm_dot_prod_q15).
 */
void matrix_polyphase_reorder_q15(q15_t *matrix, size_t width, size_t height)
{
    matrix_transpose_q15(matrix, width, height);
    matrix_reverse_rows_q15(matrix, height, width);
}

With this, we can implement a FIR polyphase filter bank (PFB). When we want to use the sub-filter at index P we just have to run an FIR (or arm_dot_prod_q15) starting at index P * S and of length S:

/* firpfb->coeffs contains the coefficients reordered by matrix_polyphase_reorder_q15 */
q63_t firpfb_execute(firpfb_t *firpfb, uint16_t index)
{
    q15_t *coeffs = &firpfb->coeffs[index * firpfb->subfilter_len];
    q15_t *data = window_read(&firpfb->win);
    q63_t result;

    assert(index < firpfb->num_subfilters);

    arm_dot_prod_q15(data, coeffs, firpfb->subfilter_len, &result);
    return result;
}

To implement a rational resampler with interpolation factor pand decimation factor q, we have to use a FIRPFB with p subfilters:

/* Input: array of q samples. Output: array of p samples. */
void rresamp_execute(rresamp_t *rresamp, const q15_t *input, q15_t *output)
{
    uint16_t i, n = 0, index = 0;

    for (i = 0; i < rresamp->q; i++) {
        firpfb_push(&rresamp->firpfb, input[i]);

        while (index < rresamp->p) {
            output[n++] = firpfb_execute(&rresamp->firpfb, index);
            index += rresamp->q;
        }

        index -= rresamp->p;
    }
}