Multirate on Hardware

Decimation, CIC filters, and polyphase FIR on ESP32-S3, STM32F4, and an FPGA (VHDL)

Multirate processing is where embedded DSP shines. A sigma-delta ADC oversampling at 3 MHz must be decimated to 48 kHz before any meaningful processing can happen, and the decimation filter runs at the input rate, so efficiency matters enormously. The CIC filter exists precisely because of this constraint: it achieves high decimation ratios using only additions, no multiplications.

This page covers practical implementations of decimation, CIC, and polyphase FIR on microcontrollers. For the theory (spectral folding, polyphase decomposition, multi-stage design), see the main multirate page.

CIC decimator: the multiply-free filter

The CIC (Cascaded Integrator-Comb) filter is the quintessential “exists because of hardware constraints” algorithm (Hogenauer 1981). A $K$-th order CIC decimating by $M$:

\[H(z) = \left(\frac{1 - z^{-M}}{1 - z^{-1}}\right)^K\]

decomposes into integrators running at the high rate and comb sections at the low rate. The entire structure requires zero multiplications, only additions, subtractions, and delays.

Bare-metal CIC implementation

A first-order CIC ($K=1$) decimating by $M$:

#include <stdint.h>

typedef struct {
    int32_t integrator;   // runs at high rate
    int32_t comb_delay;   // runs at low rate
    int count;
    int M;                // decimation factor
} CicDecimator;

void cic_init(CicDecimator *cic, int M) {
    cic->integrator = 0;
    cic->comb_delay = 0;
    cic->count = 0;
    cic->M = M;
}

// Call at the HIGH sample rate. Returns 1 when a decimated output is ready.
int cic_process(CicDecimator *cic, int32_t in, int32_t *out) {
    // Integrator: accumulate every input sample
    cic->integrator += in;

    if (++cic->count >= cic->M) {
        cic->count = 0;
        // Comb: difference at the low rate
        *out = cic->integrator - cic->comb_delay;
        cic->comb_delay = cic->integrator;
        return 1;  // output ready
    }
    return 0;  // no output this sample
}

For a higher-order CIC ($K > 1$), cascade $K$ integrators at the high rate and $K$ combs at the low rate:

#define CIC_ORDER 3

typedef struct {
    int32_t integrators[CIC_ORDER];
    int32_t comb_delays[CIC_ORDER];
    int count;
    int M;
} CicDecimatorN;

void cic_n_init(CicDecimatorN *cic, int M) {
    for (int k = 0; k < CIC_ORDER; k++) {
        cic->integrators[k] = 0;
        cic->comb_delays[k] = 0;
    }
    cic->count = 0;
    cic->M = M;
}

int cic_n_process(CicDecimatorN *cic, int32_t in, int32_t *out) {
    // K integrators at the high rate
    int32_t val = in;
    for (int k = 0; k < CIC_ORDER; k++) {
        cic->integrators[k] += val;
        val = cic->integrators[k];
    }

    if (++cic->count >= cic->M) {
        cic->count = 0;
        // K combs at the low rate
        for (int k = 0; k < CIC_ORDER; k++) {
            int32_t tmp = val;
            val = val - cic->comb_delays[k];
            cic->comb_delays[k] = tmp;
        }
        *out = val;
        return 1;
    }
    return 0;
}

Integer overflow in CIC

The integrator word width must accommodate the worst-case accumulation. For a $K$-th order CIC decimating by $M$ with $B$-bit input, the required register width is:

\[B_{\text{out}} = B_{\text{in}} + \lceil K\log_2 M\rceil\]

This is the DC gain in bits: the CIC’s DC gain is $M^K$, so it adds $\log_2(M^K) = K\log_2 M$ bits, rounded up. For 16-bit input, $K=3$, $M=16$: $B_{\text{out}} = 16 + \lceil 3\log_2 16\rceil = 16 + 12 = 28$ bits. A 32-bit int32_t suffices. For larger $M$ or $K$, use 64-bit accumulators. (The ceiling is taken on the whole product $K\log_2 M$, not stage-by-stage; the two differ by up to one bit when $M$ is not a power of two.)

The overflows in the integrator are intentional: a CIC relies on modular (wrap-around) arithmetic, and the comb stage subtracts the wrap back out exactly. Unsigned integers are the safe choice, because C defines unsigned overflow as wrap-around; signed overflow (the int32_t above) is technically undefined behaviour, so either switch the integrators and combs to uint32_t (casting the comb output back to signed) or compile with -fwrapv. An overflow-trapping build (-ftrapv, UBSan) will fault on the wrap, so do not enable one for this code.

CIC on an FPGA: registers and adders

The CIC’s real home is gateware. On an FPGA it maps to almost nothing: each integrator is one register and one adder, each comb is one register and one subtractor, and there is not a single multiplier, so no DSP slices are consumed at all. Because every stage is a registered add, the whole chain pipelines trivially and runs at the sample clock even when that clock is in the hundreds of MHz. This is why the decimator at the front of a software-radio receiver is a CIC and not an FIR.

The overflow story from the callout above becomes easier here, not harder: a signed value in VHDL’s numeric_std wraps on overflow by default, which is exactly the modular arithmetic a CIC needs. You do nothing special; you just make sure the accumulator is wide enough and never ask the tools to saturate.

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;

-- Parameterised CIC decimator: K integrator + comb stages, decimate by M.
-- Two's-complement wrap-around in the integrators is INTENTIONAL and required:
-- numeric_std 'signed' addition wraps, and the combs undo the wrap exactly.
-- Do not saturate.
entity cic_decimator is
  generic (
    B_IN  : positive := 16;   -- input word width
    K     : positive := 3;    -- CIC order (number of stages)
    M     : positive := 8;    -- decimation factor
    B_OUT : positive := 25    -- = B_IN + ceil(K*log2 M); here 16 + ceil(3*log2 8) = 25
  );
  port (
    clk     : in  std_logic;                  -- high (input-rate) clock
    rst     : in  std_logic;                  -- synchronous reset, active high
    x_in    : in  signed(B_IN-1 downto 0);    -- one input sample per clk
    y_out   : out signed(B_OUT-1 downto 0);   -- decimated output
    y_valid : out std_logic                   -- '1' the clock a new y_out is ready
  );
end entity;

architecture rtl of cic_decimator is
  type word_array is array (natural range <>) of signed(B_OUT-1 downto 0);
  signal integ : word_array(0 to K-1) := (others => (others => '0'));
  signal comb  : word_array(0 to K-1) := (others => (others => '0'));
  signal count : integer range 0 to M-1 := 0;
begin
  process(clk)
    variable val, prev : signed(B_OUT-1 downto 0);
  begin
    if rising_edge(clk) then
      if rst = '1' then
        integ <= (others => (others => '0'));
        comb  <= (others => (others => '0'));
        count <= 0; y_valid <= '0'; y_out <= (others => '0');
      else
        -- Integrator cascade: every clock, at the high rate.
        val := resize(x_in, B_OUT);            -- sign-extend the input
        for k in 0 to K-1 loop
          -- integ(k) here reads its OLD (pre-update) value, so val becomes
          -- exactly the new integ(k) and cascades to the next stage this clock.
          integ(k) <= integ(k) + val;          -- wraps mod 2**B_OUT (intended)
          val := integ(k) + val;               -- feed the running sum to the next stage
        end loop;

        if count = M-1 then                    -- decimated tick: run the combs
          count <= 0; y_valid <= '1';
          for k in 0 to K-1 loop
            prev    := comb(k);
            comb(k) <= val;                    -- remember this stage's input
            val     := val - prev;             -- difference; wrap cancels the integ wrap
          end loop;
          y_out <= val;
        else
          count <= count + 1; y_valid <= '0';
        end if;
      end if;
    end if;
  end process;
end architecture;

The one design-time number is B_OUT. Size it from the register-growth bound above ($B_\text{in} + \lceil K\log_2 M\rceil$) and the integrators can overflow as violently as they like; the combs still reconstruct the exact result. That claim is not folklore here: the main page’s cic_decimate and a bit-accurate model of this exact structure are cross-checked in test_multirate.py, feeding a full-scale signal that drives the 3-stage integrator state to 42 bits and confirming the 25-bit wrapped output is identical to an ideal unbounded-integer CIC. (The example uses a within-clock combinational cascade for clarity; for the very highest clock rates you register between stages, which only adds a few clocks of latency and does not change the arithmetic.)

Not simulated in this repo’s CI

The VHDL above is a faithful transcription of the algorithm verified in Python, but the workshop has no HDL toolchain in CI, so it is checked by inspection and by the equivalent Python model, not by a VHDL simulator. Run it through your own ghdl/vendor flow before trusting it on hardware.

Why CIC first, FIR second

A CIC filter has a sinc-like frequency response with passband droop and poor stopband rejection. It is not suitable as a standalone filter for most applications. The standard approach is a two-stage pipeline:

CIC for coarse decimation (e.g., $\div 16$ or $\div 32$), cheap, runs at the high rate
FIR for fine decimation and passband compensation, more expensive per tap, but runs at the reduced rate

This combination gives the sharp cutoff of an FIR filter with the efficiency of a CIC for the high-rate heavy lifting.

STM32F4 (NUCLEO-F446RE): CMSIS-DSP decimation

CMSIS-DSP provides arm_fir_decimate_f32 which combines an FIR anti-aliasing filter with downsampling in a single optimised call. It computes only the output samples that will be kept, the efficiency gain of polyphase decomposition, without the user having to implement it.

FIR decimation

#include "arm_math.h"

#define BLOCK_SIZE_IN   256    // input block (high rate)
#define DECIM_FACTOR    4
#define BLOCK_SIZE_OUT  (BLOCK_SIZE_IN / DECIM_FACTOR)
#define NUM_TAPS        32     // anti-aliasing FIR length

static float32_t fir_state[NUM_TAPS + BLOCK_SIZE_IN - 1];
static float32_t fir_coeffs[NUM_TAPS];  // lowpass at 1/(2*DECIM_FACTOR)
static arm_fir_decimate_instance_f32 decim;

static float32_t input[BLOCK_SIZE_IN];
static float32_t output[BLOCK_SIZE_OUT];

void init_decimator(void) {
    // Design coefficients in Python, copy here
    // scipy.signal.firwin(32, 1.0/8, window='hamming')
    arm_fir_decimate_init_f32(&decim, NUM_TAPS, DECIM_FACTOR,
                               fir_coeffs, fir_state, BLOCK_SIZE_IN);
}

void process_block(void) {
    arm_fir_decimate_f32(&decim, input, output, BLOCK_SIZE_IN);
    // output[] now contains BLOCK_SIZE_OUT = 64 samples at fs/4
}

FIR interpolation

The inverse operation uses arm_fir_interpolate_f32:

#define INTERP_FACTOR   4
#define BLOCK_SIZE_LO   64
#define BLOCK_SIZE_HI   (BLOCK_SIZE_LO * INTERP_FACTOR)
#define INTERP_TAPS     32   // per phase: total filter length = INTERP_TAPS * INTERP_FACTOR

static float32_t interp_state[INTERP_TAPS + BLOCK_SIZE_LO - 1];
static float32_t interp_coeffs[INTERP_TAPS * INTERP_FACTOR];
static arm_fir_interpolate_instance_f32 interp;

static float32_t lo_rate[BLOCK_SIZE_LO];
static float32_t hi_rate[BLOCK_SIZE_HI];

void init_interpolator(void) {
    // CMSIS wants the TOTAL coefficient count (must be a multiple of INTERP_FACTOR),
    // not the per-phase count; it derives phaseLength = numTaps / L internally.
    arm_fir_interpolate_init_f32(&interp, INTERP_FACTOR, INTERP_TAPS * INTERP_FACTOR,
                                  interp_coeffs, interp_state, BLOCK_SIZE_LO);
}

void upsample_block(void) {
    arm_fir_interpolate_f32(&interp, lo_rate, hi_rate, BLOCK_SIZE_LO);
    // hi_rate[] now contains 256 samples at 4x the input rate
}

Tip

CMSIS-DSP’s numTaps argument is the total number of coefficients and must be a multiple of L; it computes the per-phase length numTaps / L internally (and arm_fir_interpolate_init_f32 returns ARM_MATH_LENGTH_ERROR otherwise). Here INTERP_TAPS (32) is our chosen per-phase length, so the call passes INTERP_TAPS * INTERP_FACTOR (128) as numTaps and the state buffer is sized phaseLength + blockSize - 1. Design the full 128-tap prototype in Python and pass all of it; CMSIS-DSP handles the polyphase decomposition.

CIC + FIR two-stage pipeline on STM32F4

For high decimation factors, combine a hand-written CIC first stage with CMSIS-DSP FIR second stage:

// Example: decimate from 1 MHz to 8 kHz (factor 125 = 25 x 5)
// Stage 1: CIC decimate by 25 (1 MHz -> 40 kHz)
// Stage 2: FIR decimate by 5  (40 kHz -> 8 kHz)

#define FIR_BLOCK_SIZE 60   // must be a multiple of the stage-2 decimation factor (5)
#define CIC_GAIN 15625   // M^K = 25^3 for the order-3 CIC

CicDecimatorN cic_stage;    // order-3 CIC, decimating by M = 25
arm_fir_decimate_instance_f32 fir_stage;
static float32_t fir_input_buf[FIR_BLOCK_SIZE];
static float32_t fir_output_buf[FIR_BLOCK_SIZE / 5];   // second stage decimates by 5
static int fir_idx = 0;

void process_sample_high_rate(int32_t adc_sample) {
    int32_t cic_out;
    if (cic_n_process(&cic_stage, adc_sample, &cic_out)) {
        // CIC output at 40 kHz, feed to FIR buffer
        fir_input_buf[fir_idx++] = (float32_t)cic_out / CIC_GAIN;
        if (fir_idx >= FIR_BLOCK_SIZE) {
            arm_fir_decimate_f32(&fir_stage, fir_input_buf,
                                  fir_output_buf, FIR_BLOCK_SIZE);
            fir_idx = 0;
            // fir_output_buf[] is now at 8 kHz
        }
    }
}

The CIC gain for order $K$ and decimation $M$ is $M^K$ (here $25^3 = 15625$). Divide by this after the CIC stage to normalise.

Performance budget (NUCLEO-F446RE, 180 MHz)

Stage	Rate	Operations per output sample	Est. cycles	Time
CIC order-3, $\div 25$	1 MHz in, 40 kHz out	3 adds per input = 75 adds/out	~100	0.6 us
FIR 32-tap, $\div 5$	40 kHz in, 8 kHz out	32 MACs per output (6-7 taps per phase $\times$ 5 phases)	~40	0.2 us
Total per 8 kHz output			~140	0.8 us
Available per output (8 kHz)			22,500	125 us
Utilisation				0.6%

Without the CIC, a single-stage 32-tap FIR at 1 MHz would cost 32 MACs per input sample = 32M MACs/s. With CIC first, the FIR runs at 40 kHz and costs only 256K MACs/s, a 125x reduction.

ESP32-S3: CIC + ESP-DSP FIR

The ESP32-S3 lacks CMSIS-DSP but provides ESP-DSP with optimised FIR functions. The CIC code is pure C and runs identically on both platforms.

ESP-DSP FIR decimation

#include "dsps_fir.h"

#define FIR_TAPS    32
#define DECIM       4
#define BLOCK_IN    256
#define BLOCK_OUT   (BLOCK_IN / DECIM)

static fir_f32_t fir;
static float coeffs[FIR_TAPS];
static float delay[FIR_TAPS];
static float input[BLOCK_IN];
static float output[BLOCK_IN];  // filtered at high rate

void init_fir(void) {
    dsps_fir_init_f32(&fir, coeffs, delay, FIR_TAPS);
}

void decimate_block(float *in, float *out, int n_in) {
    // Filter at high rate
    dsps_fir_f32(&fir, in, output, n_in);

    // Downsample: keep every DECIM-th sample
    for (int i = 0; i < n_in / DECIM; i++)
        out[i] = output[i * DECIM];
}

Note

ESP-DSP’s dsps_fir_f32 does not have a built-in decimation mode like CMSIS-DSP’s arm_fir_decimate_f32. The filter runs at the full input rate, and downsampling is done separately. This is less efficient (computes samples that are discarded), but for moderate decimation factors the overhead is acceptable. For large factors, use the CIC + FIR two-stage approach.

Two-stage pipeline on ESP32-S3

// Decimate from 48 kHz to 8 kHz (factor 6 = 3 x 2)
// Stage 1: CIC order-3 (CIC_ORDER), decimate by 3 (48 kHz -> 16 kHz)
// Stage 2: FIR 32-tap, decimate by 2 (16 kHz -> 8 kHz)

CicDecimatorN cic;
fir_f32_t fir;

void i2s_process_block(int32_t *samples, int n) {
    static float fir_buf[256];
    static float fir_out[256];
    static int fir_idx = 0;  // static: persists across calls

    for (int i = 0; i < n; i++) {
        // INMP441: 24-bit audio in bits 31:8 of the 32-bit I2S frame
        float x = (float)(samples[i] >> 8) / 8388608.0f;
        int32_t cic_out;
        if (cic_n_process(&cic, (int32_t)(x * 32768), &cic_out)) {
            if (fir_idx < 256) {  // bounds guard
                fir_buf[fir_idx++] = (float)cic_out / (3 * 3 * 3);  // CIC gain = M^K = 3^3 = 27
            }

            if (fir_idx >= 64) {
                // FIR filter
                dsps_fir_f32(&fir, fir_buf, fir_out, 64);
                // Downsample by 2
                for (int j = 0; j < 32; j++) {
                    float out_sample = fir_out[j * 2];
                    // out_sample is now at 8 kHz, process further
                }
                fir_idx = 0;
            }
        }
    }
}

Performance budget (ESP32-S3, 240 MHz)

Each stage runs at its own rate, so the figures below are normalised to one 8 kHz output sample: the CIC (48 kHz) runs 6 times per output, the FIR (16 kHz) twice, the downsampler once.

Stage	Rate	Operations per 8 kHz output	Est. cycles	Time
CIC order-3, $\div 3$	48 kHz in	3 adds/sample $\times$ 6 = 18 adds	~36	0.15 us
FIR 32-tap (ESP-DSP)	16 kHz	32 MACs/sample $\times$ 2 = 64 MACs	~120	0.50 us
Downsample $\div 2$	16 kHz → 8 kHz	1 copy	~2	0.01 us
Total per 8 kHz output			~158	~0.66 us
Available per output (8 kHz)			30,000	125 us
Utilisation				~0.5%

Real-time rate conversion: two clocks

Everything above processes blocks already sitting in memory. On real hardware a sample-rate converter sits between two clocks, and unlike the single-rate biquad loop, its input and output tick at different rates. That asymmetry is the whole design problem, and it is what the original practicum’s decimator and interpolator assignments were really about.

Decimator: high-rate in, low-rate out

A decimator samples fast and emits slow. A timer triggers the ADC at the high input rate, each conversion feeds the CIC or FIR, and an output appears only every $M$-th sample. The neat consequence is that you do not need a second clock for the output: writing the DAC exactly when the filter produces a sample paces it at $f_\text{in}/M$ automatically, one DAC write per $M$ ADC triggers.

// Real-time decimator on a NUCLEO-F446RE (the practicum's SRD assignment).
// TIM6 TRGO triggers ADC1 at the high input rate; the end-of-conversion
// interrupt feeds each sample to the CIC, which emits one output every M-th
// sample. Writing the DAC only when an output is ready paces it at f_in / M,
// so no second timer is needed.
//
// Handles from CubeMX: extern ADC_HandleTypeDef hadc1;
//                      extern DAC_HandleTypeDef hdac;
//                      extern TIM_HandleTypeDef htim6;

#define M        8     // decimation factor for this example
#define CIC_GAIN 512   // M^K = 8^3 for the order-3 CIC above

static CicDecimatorN cic;   // order-3 CIC decimating by M

void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc) {
    int32_t x = (int32_t)HAL_ADC_GetValue(hadc) - 2048;   // 12-bit, mid-rail removed
    int32_t y;
    if (cic_n_process(&cic, x, &y)) {
        y /= CIC_GAIN;                                     // normalise by M^K
        int32_t code = y + 2048;
        if (code < 0) code = 0; else if (code > 4095) code = 4095;
        HAL_DAC_SetValue(&hdac, DAC_CHANNEL_1, DAC_ALIGN_12B_R, (uint16_t)code);
    }
}

void start_decimator(void) {
    cic_n_init(&cic, M);
    HAL_DAC_Start(&hdac, DAC_CHANNEL_1);
    HAL_ADC_Start_IT(&hadc1);     // armed; TIM6 TRGO triggers each conversion at f_in
    HAL_TIM_Base_Start(&htim6);
}

Interpolator: low-rate in, high-rate out

The interpolator is the mirror image and the harder of the two, because the output is faster than the input. Each input sample drives the polyphase filter to produce $L$ output samples, and those have to leave at the high rate $L \cdot f_\text{in}$. You cannot derive that fast output clock from the slow input interrupt, so the DAC needs its own timer (or a timer-paced DAC DMA) running at $L \cdot f_\text{in}$, draining an output FIFO that the low-rate input keeps filled. The catch is starvation: if the input falls behind and the FIFO empties, the DAC repeats stale samples and the output clicks. The fix is to stay a block ahead, computing the next $L$ outputs as soon as an input arrives rather than on demand.

The shared lesson

With two clocks you always need a small FIFO between the rate-converting filter and whichever side runs faster, and both the ADC trigger and the DAC update must be hardware-timed for the reasons in DSP on a microcontroller. To check the result, generate a tone on a laptop, run it through the converter, and compare the output spectrum against the same conversion done in Python: the passband content should line up, shifted to the new rate, with the anti-aliasing or anti-imaging filter’s roll-off where you designed it.

Platform comparison

Feature	STM32F4 (NUCLEO-F446RE)	ESP32-S3
FIR decimation library	`arm_fir_decimate_f32` (polyphase, optimised)	`dsps_fir_f32` + manual downsample
FIR interpolation library	`arm_fir_interpolate_f32`	No built-in (manual implementation)
Fixed-point support	Q15/Q31 variants with SIMD	s16 FIR with SIMD (~2.4x speedup)
CIC implementation	Hand-written (same C on both)	Hand-written (same C on both)
ADC for high-rate input	Up to 2.4 Msps (12-bit, DMA)	~100 ksps practical (noisy ADC)
Best for	High-rate sensor decimation, SDR front-ends	Audio-rate multirate, I2S codec interfacing

Recommendation

STM32F4 for high-rate decimation: the Cortex-M4’s ADC can sample at MHz rates, and CMSIS-DSP’s polyphase FIR decimation avoids computing discarded samples. Choose this for sigma-delta post-processing, vibration analysis, or SDR.
ESP32-S3 for audio-rate conversion: 48 → 16 → 8 kHz pipelines for voice processing, or upsampling for DAC output. The s16 FIR SIMD gives a 2.4x speedup over float.

References

Hogenauer, Eugene B. 1981. “An Economical Class of Digital Filters for Decimation and Interpolation.” IEEE Transactions on Acoustics, Speech, and Signal Processing 29 (2): 155–62.

--- title: "Multirate on Hardware" subtitle: "Decimation, CIC filters, and polyphase FIR on ESP32-S3, STM32F4, and an FPGA (VHDL)" aliases: - /topics/multirate/embedded.html --- Multirate processing is where embedded DSP shines. A sigma-delta ADC oversampling at 3 MHz must be decimated to 48 kHz before any meaningful processing can happen, and the decimation filter runs at the *input* rate, so efficiency matters enormously. The CIC filter exists precisely because of this constraint: it achieves high decimation ratios using only additions, no multiplications. This page covers practical implementations of decimation, CIC, and polyphase FIR on microcontrollers. For the theory (spectral folding, polyphase decomposition, multi-stage design), see the [main multirate page](index.qmd). <hr> ## CIC decimator: the multiply-free filter The CIC (Cascaded Integrator-Comb) filter is the quintessential "exists because of hardware constraints" algorithm [@hogenauer1981cascaded]. A $K$-th order CIC decimating by $M$: $$H(z) = \left(\frac{1 - z^{-M}}{1 - z^{-1}}\right)^K$$ decomposes into integrators running at the high rate and comb sections at the low rate. The entire structure requires **zero multiplications**, only additions, subtractions, and delays. ### Bare-metal CIC implementation A first-order CIC ($K=1$) decimating by $M$: ```c #include <stdint.h> typedef struct { int32_t integrator; // runs at high rate int32_t comb_delay; // runs at low rate int count; int M; // decimation factor } CicDecimator; void cic_init(CicDecimator *cic, int M) { cic->integrator = 0; cic->comb_delay = 0; cic->count = 0; cic->M = M; } // Call at the HIGH sample rate. Returns 1 when a decimated output is ready. int cic_process(CicDecimator *cic, int32_t in, int32_t *out) { // Integrator: accumulate every input sample cic->integrator += in; if (++cic->count >= cic->M) { cic->count = 0; // Comb: difference at the low rate *out = cic->integrator - cic->comb_delay; cic->comb_delay = cic->integrator; return 1; // output ready } return 0; // no output this sample } ``` For a higher-order CIC ($K > 1$), cascade $K$ integrators at the high rate and $K$ combs at the low rate: ```c #define CIC_ORDER 3 typedef struct { int32_t integrators[CIC_ORDER]; int32_t comb_delays[CIC_ORDER]; int count; int M; } CicDecimatorN; void cic_n_init(CicDecimatorN *cic, int M) { for (int k = 0; k < CIC_ORDER; k++) { cic->integrators[k] = 0; cic->comb_delays[k] = 0; } cic->count = 0; cic->M = M; } int cic_n_process(CicDecimatorN *cic, int32_t in, int32_t *out) { // K integrators at the high rate int32_t val = in; for (int k = 0; k < CIC_ORDER; k++) { cic->integrators[k] += val; val = cic->integrators[k]; } if (++cic->count >= cic->M) { cic->count = 0; // K combs at the low rate for (int k = 0; k < CIC_ORDER; k++) { int32_t tmp = val; val = val - cic->comb_delays[k]; cic->comb_delays[k] = tmp; } *out = val; return 1; } return 0; } ``` ::: {.callout-warning title="Integer overflow in CIC"} The integrator word width must accommodate the worst-case accumulation. For a $K$-th order CIC decimating by $M$ with $B$-bit input, the required register width is: $$B_{\text{out}} = B_{\text{in}} + \lceil K\log_2 M\rceil$$ This is the DC gain in bits: the CIC's DC gain is $M^K$, so it adds $\log_2(M^K) = K\log_2 M$ bits, rounded up. For 16-bit input, $K=3$, $M=16$: $B_{\text{out}} = 16 + \lceil 3\log_2 16\rceil = 16 + 12 = 28$ bits. A 32-bit `int32_t` suffices. For larger $M$ or $K$, use 64-bit accumulators. (The ceiling is taken on the whole product $K\log_2 M$, not stage-by-stage; the two differ by up to one bit when $M$ is not a power of two.) The overflows in the integrator are **intentional**: a CIC relies on modular (wrap-around) arithmetic, and the comb stage subtracts the wrap back out exactly. Unsigned integers are the safe choice, because C defines unsigned overflow as wrap-around; signed overflow (the `int32_t` above) is technically undefined behaviour, so either switch the integrators and combs to `uint32_t` (casting the comb output back to signed) or compile with `-fwrapv`. An overflow-trapping build (`-ftrapv`, UBSan) will fault on the wrap, so do not enable one for this code. ::: ### CIC on an FPGA: registers and adders The CIC's real home is gateware. On an FPGA it maps to almost nothing: each integrator is **one register and one adder**, each comb is **one register and one subtractor**, and there is not a single multiplier, so no DSP slices are consumed at all. Because every stage is a registered add, the whole chain pipelines trivially and runs at the sample clock even when that clock is in the hundreds of MHz. This is why the decimator at the front of a software-radio receiver is a CIC and not an FIR. The overflow story from the callout above becomes *easier* here, not harder: a `signed` value in VHDL's `numeric_std` wraps on overflow by default, which is exactly the modular arithmetic a CIC needs. You do nothing special; you just make sure the accumulator is wide enough and never ask the tools to saturate. ```vhdl library ieee; use ieee.std_logic_1164.all; use ieee.numeric_std.all; -- Parameterised CIC decimator: K integrator + comb stages, decimate by M. -- Two's-complement wrap-around in the integrators is INTENTIONAL and required: -- numeric_std 'signed' addition wraps, and the combs undo the wrap exactly. -- Do not saturate. entity cic_decimator is generic ( B_IN : positive := 16; -- input word width K : positive := 3; -- CIC order (number of stages) M : positive := 8; -- decimation factor B_OUT : positive := 25 -- = B_IN + ceil(K*log2 M); here 16 + ceil(3*log2 8) = 25 ); port ( clk : in std_logic; -- high (input-rate) clock rst : in std_logic; -- synchronous reset, active high x_in : in signed(B_IN-1 downto 0); -- one input sample per clk y_out : out signed(B_OUT-1 downto 0); -- decimated output y_valid : out std_logic -- '1' the clock a new y_out is ready ); end entity; architecture rtl of cic_decimator is type word_array is array (natural range <>) of signed(B_OUT-1 downto 0); signal integ : word_array(0 to K-1) := (others => (others => '0')); signal comb : word_array(0 to K-1) := (others => (others => '0')); signal count : integer range 0 to M-1 := 0; begin process(clk) variable val, prev : signed(B_OUT-1 downto 0); begin if rising_edge(clk) then if rst = '1' then integ <= (others => (others => '0')); comb <= (others => (others => '0')); count <= 0; y_valid <= '0'; y_out <= (others => '0'); else -- Integrator cascade: every clock, at the high rate. val := resize(x_in, B_OUT); -- sign-extend the input for k in 0 to K-1 loop -- integ(k) here reads its OLD (pre-update) value, so val becomes -- exactly the new integ(k) and cascades to the next stage this clock. integ(k) <= integ(k) + val; -- wraps mod 2**B_OUT (intended) val := integ(k) + val; -- feed the running sum to the next stage end loop; if count = M-1 then -- decimated tick: run the combs count <= 0; y_valid <= '1'; for k in 0 to K-1 loop prev := comb(k); comb(k) <= val; -- remember this stage's input val := val - prev; -- difference; wrap cancels the integ wrap end loop; y_out <= val; else count <= count + 1; y_valid <= '0'; end if; end if; end if; end process; end architecture; ``` The one design-time number is `B_OUT`. Size it from the register-growth bound above ($B_\text{in} + \lceil K\log_2 M\rceil$) and the integrators can overflow as violently as they like; the combs still reconstruct the exact result. That claim is not folklore here: the [main page's](index.qmd) `cic_decimate` and a bit-accurate model of this exact structure are cross-checked in `test_multirate.py`, feeding a full-scale signal that drives the 3-stage integrator state to 42 bits and confirming the 25-bit wrapped output is *identical* to an ideal unbounded-integer CIC. (The example uses a within-clock combinational cascade for clarity; for the very highest clock rates you register between stages, which only adds a few clocks of latency and does not change the arithmetic.) ::: {.callout-note title="Not simulated in this repo's CI"} The VHDL above is a faithful transcription of the algorithm verified in Python, but the workshop has no HDL toolchain in CI, so it is checked by inspection and by the equivalent Python model, not by a VHDL simulator. Run it through your own `ghdl`/vendor flow before trusting it on hardware. ::: ### Why CIC first, FIR second A CIC filter has a sinc-like frequency response with passband droop and poor stopband rejection. It is not suitable as a standalone filter for most applications. The standard approach is a **two-stage pipeline**: 1. **CIC** for coarse decimation (e.g., $\div 16$ or $\div 32$), cheap, runs at the high rate 2. **FIR** for fine decimation and passband compensation, more expensive per tap, but runs at the reduced rate This combination gives the sharp cutoff of an FIR filter with the efficiency of a CIC for the high-rate heavy lifting. <hr> ## STM32F4 (NUCLEO-F446RE): CMSIS-DSP decimation CMSIS-DSP provides `arm_fir_decimate_f32` which combines an FIR anti-aliasing filter with downsampling in a single optimised call. It computes only the output samples that will be kept, the efficiency gain of polyphase decomposition, without the user having to implement it. ### FIR decimation ```c #include "arm_math.h" #define BLOCK_SIZE_IN 256 // input block (high rate) #define DECIM_FACTOR 4 #define BLOCK_SIZE_OUT (BLOCK_SIZE_IN / DECIM_FACTOR) #define NUM_TAPS 32 // anti-aliasing FIR length static float32_t fir_state[NUM_TAPS + BLOCK_SIZE_IN - 1]; static float32_t fir_coeffs[NUM_TAPS]; // lowpass at 1/(2*DECIM_FACTOR) static arm_fir_decimate_instance_f32 decim; static float32_t input[BLOCK_SIZE_IN]; static float32_t output[BLOCK_SIZE_OUT]; void init_decimator(void) { // Design coefficients in Python, copy here // scipy.signal.firwin(32, 1.0/8, window='hamming') arm_fir_decimate_init_f32(&decim, NUM_TAPS, DECIM_FACTOR, fir_coeffs, fir_state, BLOCK_SIZE_IN); } void process_block(void) { arm_fir_decimate_f32(&decim, input, output, BLOCK_SIZE_IN); // output[] now contains BLOCK_SIZE_OUT = 64 samples at fs/4 } ``` ### FIR interpolation The inverse operation uses `arm_fir_interpolate_f32`: ```c #define INTERP_FACTOR 4 #define BLOCK_SIZE_LO 64 #define BLOCK_SIZE_HI (BLOCK_SIZE_LO * INTERP_FACTOR) #define INTERP_TAPS 32 // per phase: total filter length = INTERP_TAPS * INTERP_FACTOR static float32_t interp_state[INTERP_TAPS + BLOCK_SIZE_LO - 1]; static float32_t interp_coeffs[INTERP_TAPS * INTERP_FACTOR]; static arm_fir_interpolate_instance_f32 interp; static float32_t lo_rate[BLOCK_SIZE_LO]; static float32_t hi_rate[BLOCK_SIZE_HI]; void init_interpolator(void) { // CMSIS wants the TOTAL coefficient count (must be a multiple of INTERP_FACTOR), // not the per-phase count; it derives phaseLength = numTaps / L internally. arm_fir_interpolate_init_f32(&interp, INTERP_FACTOR, INTERP_TAPS * INTERP_FACTOR, interp_coeffs, interp_state, BLOCK_SIZE_LO); } void upsample_block(void) { arm_fir_interpolate_f32(&interp, lo_rate, hi_rate, BLOCK_SIZE_LO); // hi_rate[] now contains 256 samples at 4x the input rate } ``` ::: {.callout-tip} CMSIS-DSP's `numTaps` argument is the **total** number of coefficients and must be a multiple of `L`; it computes the per-phase length `numTaps / L` internally (and `arm_fir_interpolate_init_f32` returns `ARM_MATH_LENGTH_ERROR` otherwise). Here `INTERP_TAPS` (32) is our chosen *per-phase* length, so the call passes `INTERP_TAPS * INTERP_FACTOR` (128) as `numTaps` and the state buffer is sized `phaseLength + blockSize - 1`. Design the full 128-tap prototype in Python and pass all of it; CMSIS-DSP handles the polyphase decomposition. ::: ### CIC + FIR two-stage pipeline on STM32F4 For high decimation factors, combine a hand-written CIC first stage with CMSIS-DSP FIR second stage: ```c // Example: decimate from 1 MHz to 8 kHz (factor 125 = 25 x 5) // Stage 1: CIC decimate by 25 (1 MHz -> 40 kHz) // Stage 2: FIR decimate by 5 (40 kHz -> 8 kHz) #define FIR_BLOCK_SIZE 60 // must be a multiple of the stage-2 decimation factor (5) #define CIC_GAIN 15625 // M^K = 25^3 for the order-3 CIC CicDecimatorN cic_stage; // order-3 CIC, decimating by M = 25 arm_fir_decimate_instance_f32 fir_stage; static float32_t fir_input_buf[FIR_BLOCK_SIZE]; static float32_t fir_output_buf[FIR_BLOCK_SIZE / 5]; // second stage decimates by 5 static int fir_idx = 0; void process_sample_high_rate(int32_t adc_sample) { int32_t cic_out; if (cic_n_process(&cic_stage, adc_sample, &cic_out)) { // CIC output at 40 kHz, feed to FIR buffer fir_input_buf[fir_idx++] = (float32_t)cic_out / CIC_GAIN; if (fir_idx >= FIR_BLOCK_SIZE) { arm_fir_decimate_f32(&fir_stage, fir_input_buf, fir_output_buf, FIR_BLOCK_SIZE); fir_idx = 0; // fir_output_buf[] is now at 8 kHz } } } ``` The CIC gain for order $K$ and decimation $M$ is $M^K$ (here $25^3 = 15625$). Divide by this after the CIC stage to normalise. ### Performance budget (NUCLEO-F446RE, 180 MHz) | Stage | Rate | Operations per output sample | Est. cycles | Time | |---|---|---|---|---| | CIC order-3, $\div 25$ | 1 MHz in, 40 kHz out | 3 adds per input = 75 adds/out | ~100 | 0.6 us | | FIR 32-tap, $\div 5$ | 40 kHz in, 8 kHz out | 32 MACs per output (6-7 taps per phase $\times$ 5 phases) | ~40 | 0.2 us | | **Total per 8 kHz output** | | | **~140** | **0.8 us** | | Available per output (8 kHz) | | | 22,500 | 125 us | | **Utilisation** | | | | **0.6%** | Without the CIC, a single-stage 32-tap FIR at 1 MHz would cost 32 MACs per *input* sample = 32M MACs/s. With CIC first, the FIR runs at 40 kHz and costs only 256K MACs/s, a **125x** reduction. <hr> ## ESP32-S3: CIC + ESP-DSP FIR The ESP32-S3 lacks CMSIS-DSP but provides ESP-DSP with optimised FIR functions. The CIC code is pure C and runs identically on both platforms. ### ESP-DSP FIR decimation ```c #include "dsps_fir.h" #define FIR_TAPS 32 #define DECIM 4 #define BLOCK_IN 256 #define BLOCK_OUT (BLOCK_IN / DECIM) static fir_f32_t fir; static float coeffs[FIR_TAPS]; static float delay[FIR_TAPS]; static float input[BLOCK_IN]; static float output[BLOCK_IN]; // filtered at high rate void init_fir(void) { dsps_fir_init_f32(&fir, coeffs, delay, FIR_TAPS); } void decimate_block(float *in, float *out, int n_in) { // Filter at high rate dsps_fir_f32(&fir, in, output, n_in); // Downsample: keep every DECIM-th sample for (int i = 0; i < n_in / DECIM; i++) out[i] = output[i * DECIM]; } ``` ::: {.callout-note} ESP-DSP's `dsps_fir_f32` does not have a built-in decimation mode like CMSIS-DSP's `arm_fir_decimate_f32`. The filter runs at the full input rate, and downsampling is done separately. This is less efficient (computes samples that are discarded), but for moderate decimation factors the overhead is acceptable. For large factors, use the CIC + FIR two-stage approach. ::: ### Two-stage pipeline on ESP32-S3 ```c // Decimate from 48 kHz to 8 kHz (factor 6 = 3 x 2) // Stage 1: CIC order-3 (CIC_ORDER), decimate by 3 (48 kHz -> 16 kHz) // Stage 2: FIR 32-tap, decimate by 2 (16 kHz -> 8 kHz) CicDecimatorN cic; fir_f32_t fir; void i2s_process_block(int32_t *samples, int n) { static float fir_buf[256]; static float fir_out[256]; static int fir_idx = 0; // static: persists across calls for (int i = 0; i < n; i++) { // INMP441: 24-bit audio in bits 31:8 of the 32-bit I2S frame float x = (float)(samples[i] >> 8) / 8388608.0f; int32_t cic_out; if (cic_n_process(&cic, (int32_t)(x * 32768), &cic_out)) { if (fir_idx < 256) { // bounds guard fir_buf[fir_idx++] = (float)cic_out / (3 * 3 * 3); // CIC gain = M^K = 3^3 = 27 } if (fir_idx >= 64) { // FIR filter dsps_fir_f32(&fir, fir_buf, fir_out, 64); // Downsample by 2 for (int j = 0; j < 32; j++) { float out_sample = fir_out[j * 2]; // out_sample is now at 8 kHz, process further } fir_idx = 0; } } } } ``` ### Performance budget (ESP32-S3, 240 MHz) Each stage runs at its own rate, so the figures below are normalised to one 8 kHz output sample: the CIC (48 kHz) runs 6 times per output, the FIR (16 kHz) twice, the downsampler once. | Stage | Rate | Operations per 8 kHz output | Est. cycles | Time | |---|---|---|---|---| | CIC order-3, $\div 3$ | 48 kHz in | 3 adds/sample $\times$ 6 = 18 adds | ~36 | 0.15 us | | FIR 32-tap (ESP-DSP) | 16 kHz | 32 MACs/sample $\times$ 2 = 64 MACs | ~120 | 0.50 us | | Downsample $\div 2$ | 16 kHz → 8 kHz | 1 copy | ~2 | 0.01 us | | **Total per 8 kHz output** | | | **~158** | **~0.66 us** | | Available per output (8 kHz) | | | 30,000 | 125 us | | **Utilisation** | | | | **~0.5%** | <hr> ## Real-time rate conversion: two clocks Everything above processes blocks already sitting in memory. On real hardware a sample-rate converter sits between two clocks, and unlike the [single-rate biquad loop](../09-biquad/embedded.qmd#real-time-io-clocking-samples-in-and-out), its input and output tick at *different* rates. That asymmetry is the whole design problem, and it is what the original practicum's decimator and interpolator assignments were really about. ### Decimator: high-rate in, low-rate out A decimator samples fast and emits slow. A timer triggers the ADC at the high input rate, each conversion feeds the CIC or FIR, and an output appears only every $M$-th sample. The neat consequence is that you do not need a second clock for the output: writing the DAC exactly when the filter produces a sample paces it at $f_\text{in}/M$ automatically, one DAC write per $M$ ADC triggers. ```c // Real-time decimator on a NUCLEO-F446RE (the practicum's SRD assignment). // TIM6 TRGO triggers ADC1 at the high input rate; the end-of-conversion // interrupt feeds each sample to the CIC, which emits one output every M-th // sample. Writing the DAC only when an output is ready paces it at f_in / M, // so no second timer is needed. // // Handles from CubeMX: extern ADC_HandleTypeDef hadc1; // extern DAC_HandleTypeDef hdac; // extern TIM_HandleTypeDef htim6; #define M 8 // decimation factor for this example #define CIC_GAIN 512 // M^K = 8^3 for the order-3 CIC above static CicDecimatorN cic; // order-3 CIC decimating by M void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc) { int32_t x = (int32_t)HAL_ADC_GetValue(hadc) - 2048; // 12-bit, mid-rail removed int32_t y; if (cic_n_process(&cic, x, &y)) { y /= CIC_GAIN; // normalise by M^K int32_t code = y + 2048; if (code < 0) code = 0; else if (code > 4095) code = 4095; HAL_DAC_SetValue(&hdac, DAC_CHANNEL_1, DAC_ALIGN_12B_R, (uint16_t)code); } } void start_decimator(void) { cic_n_init(&cic, M); HAL_DAC_Start(&hdac, DAC_CHANNEL_1); HAL_ADC_Start_IT(&hadc1); // armed; TIM6 TRGO triggers each conversion at f_in HAL_TIM_Base_Start(&htim6); } ``` ### Interpolator: low-rate in, high-rate out The interpolator is the mirror image and the harder of the two, because the output is *faster* than the input. Each input sample drives the polyphase filter to produce $L$ output samples, and those have to leave at the high rate $L \cdot f_\text{in}$. You cannot derive that fast output clock from the slow input interrupt, so the DAC needs its own timer (or a timer-paced DAC DMA) running at $L \cdot f_\text{in}$, draining an output FIFO that the low-rate input keeps filled. The catch is starvation: if the input falls behind and the FIFO empties, the DAC repeats stale samples and the output clicks. The fix is to stay a block ahead, computing the next $L$ outputs as soon as an input arrives rather than on demand. ### The shared lesson With two clocks you always need a small **FIFO** between the rate-converting filter and whichever side runs faster, and both the ADC trigger and the DAC update must be hardware-timed for the reasons in [DSP on a microcontroller](../12-dsp-on-a-microcontroller.qmd#getting-samples-in-and-out). To check the result, generate a tone on a laptop, run it through the converter, and compare the output spectrum against the same conversion done in Python: the passband content should line up, shifted to the new rate, with the anti-aliasing or anti-imaging filter's roll-off where you designed it. <hr> ## Platform comparison | Feature | STM32F4 (NUCLEO-F446RE) | ESP32-S3 | |---|---|---| | FIR decimation library | `arm_fir_decimate_f32` (polyphase, optimised) | `dsps_fir_f32` + manual downsample | | FIR interpolation library | `arm_fir_interpolate_f32` | No built-in (manual implementation) | | Fixed-point support | Q15/Q31 variants with SIMD | s16 FIR with SIMD (~2.4x speedup) | | CIC implementation | Hand-written (same C on both) | Hand-written (same C on both) | | ADC for high-rate input | Up to 2.4 Msps (12-bit, DMA) | ~100 ksps practical (noisy ADC) | | Best for | High-rate sensor decimation, SDR front-ends | Audio-rate multirate, I2S codec interfacing | ### Recommendation - **STM32F4** for high-rate decimation: the Cortex-M4's ADC can sample at MHz rates, and CMSIS-DSP's polyphase FIR decimation avoids computing discarded samples. Choose this for sigma-delta post-processing, vibration analysis, or SDR. - **ESP32-S3** for audio-rate conversion: 48 → 16 → 8 kHz pipelines for voice processing, or upsampling for DAC output. The s16 FIR SIMD gives a 2.4x speedup over float.