Wavelets on Hardware

The DWT as a filter-bank tree, and the lifting scheme for integer arithmetic

The continuous wavelet transform on the main page is an analysis tool: redundant, scale-by-scale, the wavelet equivalent of running a spectrogram. You would not put it on a microcontroller. The discrete wavelet transform is the embedded-friendly one, and it is friendly for a reason already in this workshop: it is a filter bank, the same two short filters and a downsample-by-two you already know, applied as a tree.

This page shows two ways to land the DWT on the metal. The first is the direct one, a recursion of FIR-plus-decimate stages, mapped onto CMSIS-DSP on the NUCLEO-F446RE and ESP-DSP on the ESP32-S3. The second is the lifting scheme, which factors the same transform into in-place predict-and-update steps that need no floating point at all, run on the smallest 8-bit parts, and invert exactly in integer arithmetic.

One level is two filters and a downsample

Recall the analysis step: convolve the signal with the low-pass scaling filter \(h\) and the high-pass wavelet filter \(g\), then keep every second output. In C that is a short, exact transcription of the Python module, periodic at the boundary so the lengths halve cleanly:

// One analysis level. in[] has n samples (n even); dec_lo / dec_hi are the
// length-L analysis filters. Writes the n/2 approximation coefficients to cA[]
// and the n/2 detail coefficients to cD[]. Periodic boundary handling.
void dwt_step(const float *in, int n,
              const float *dec_lo, const float *dec_hi, int L,
              float *cA, float *cD)
{
    int half = n / 2;
    for (int m = 0; m < half; m++) {
        float a = 0.0f, d = 0.0f;
        for (int k = 0; k < L; k++) {
            int idx = (2 * m + k) % n;     // periodic wrap
            a += in[idx] * dec_lo[k];
            d += in[idx] * dec_hi[k];
        }
        cA[m] = a;
        cD[m] = d;
    }
}

A full \(J\)-level transform calls this \(J\) times, each time on the previous approximation band, which is half as long as the one before. The details are kept; only the shrinking approximation is fed back:

// In-place multilevel DWT. work[] holds n samples; after the call the layout is
// [a_J | d_J | d_{J-1} | ... | d_1], coarsest approximation first, matching the
// Python wavedec() convention. scratch[] is an n-sample temporary. The detail
// MUST be written to scratch, not back into work[], while dwt_step is still
// reading work[]: the periodic taps reach into the upper half, so writing the
// detail there would corrupt later reads of the same level.
void wavedec(float *work, int n, const float *dec_lo, const float *dec_hi,
             int L, int levels, float *scratch)
{
    int len = n;
    for (int j = 0; j < levels; j++) {
        int half = len / 2;
        dwt_step(work, len, dec_lo, dec_hi, L, scratch, scratch + half);
        for (int i = 0; i < len; i++) work[i] = scratch[i];  // approx then detail
        len = half;                                          // recurse on approx
    }
}

The cost of the whole tree is the key fact. Each level is \(O(\text{len})\), and the lengths halve, so the total is a geometric series that sums to about \(2LN\) multiply-accumulates for an \(N\)-sample, \(L\)-tap transform: linear in \(N\), against the FFT’s \(N\log N\).

N = 1024            # samples per block
L = 4               # Daubechies-4 taps
J = 5               # decomposition levels
fs = 16000.0

# per level: low- and high-pass each produce len/2 outputs of L MACs each
dwt_macs = sum(2 * (N >> j) * L for j in range(1, J + 1))
geometric_bound = 2 * L * N                        # the closed-form ceiling 2LN
# the FFT it is usually weighed against: a radix-2 FFT runs log2(N) stages of
# N/2 butterflies each (one complex multiply plus two complex adds per butterfly)
fft_butterflies = (N // 2) * int(np.log2(N))

assert dwt_macs < geometric_bound                  # the series sits below 2LN
print(f"db4 {J}-level DWT: ~{dwt_macs:,} real MACs  (bound 2LN = {geometric_bound:,})")
print(f"radix-2 FFT for scale: {fft_butterflies:,} complex butterflies ({N}-point)")
print(f"one {N}-sample block at {fs/1000:.0f} kHz spans {1000*N/fs:.0f} ms")
db4 5-level DWT: ~7,936 real MACs  (bound 2LN = 8,192)
radix-2 FFT for scale: 5,120 complex butterflies (1024-point)
one 1024-sample block at 16 kHz spans 64 ms

The headline is the \(O(N)\) scaling against the FFT’s \(O(N\log N)\), not the raw count at one block size: a complex butterfly is several real MACs, so at \(N=1024\) the two transforms are in the same ballpark of arithmetic. Where the DWT actually pulls ahead on a microcontroller is elsewhere: no twiddle tables in flash, no bit-reversal pass, and it streams block by block. It needs the \(N\)-sample block plus an \(N\)-sample scratch (the lifting form below is truly in place), and fits comfortably inside the per-block budget on either reference part.

On the NUCLEO-F446RE: CMSIS-DSP decimating FIR

The filter-and-downsample pair is exactly what CMSIS-DSP’s decimating FIR was built for: arm_fir_decimate_f32 runs an FIR filter and keeps every \(M\)-th output in one call, with \(M = 2\) here. One instance carries the low-pass scaling filter, a second the high-pass wavelet filter, and one level is two calls.

#include "arm_math.h"

#define NUM_TAPS   4
#define BLOCK_SIZE 1024

// db4 analysis filters in natural (time) order (this page's convention; see
// _DB4_LO / _qmf_highpass in wavelets.py). Low-pass is the scaling filter;
// high-pass is its QMF mirror. ESP-DSP's dsps_fir_f32 below uses this order.
static const float32_t coeffs_lo[NUM_TAPS] =
    { 0.48296291f, 0.83651630f, 0.22414387f, -0.12940952f };
static const float32_t coeffs_hi[NUM_TAPS] =
    { -0.12940952f, -0.22414387f, 0.83651630f, -0.48296291f };
// CMSIS-DSP reads taps TIME-REVERSED, {h[N-1], ..., h[0]}. db4 is asymmetric, so
// passing the natural-order arrays would realise the wrong (time-reversed) filter
// and break the match to the Python module. These reversed copies are what
// arm_fir_decimate_init_f32 must receive.
static const float32_t coeffs_lo_rev[NUM_TAPS] =
    { -0.12940952f, 0.22414387f, 0.83651630f, 0.48296291f };
static const float32_t coeffs_hi_rev[NUM_TAPS] =
    { -0.48296291f, 0.83651630f, -0.22414387f, -0.12940952f };

// CMSIS state buffer per instance: (numTaps + blockSize - 1) floats, zeroed.
static float32_t state_lo[NUM_TAPS + BLOCK_SIZE - 1];
static float32_t state_hi[NUM_TAPS + BLOCK_SIZE - 1];
arm_fir_decimate_instance_f32 lo, hi;

// Call once at startup. numTaps must be a multiple of the decimation factor M;
// db4's 4 taps and M = 2 satisfy that. CMSIS reads coefficients time-reversed,
// so pass the *_rev arrays, not the natural-order ones.
void dwt_level_cmsis_init(void)
{
    arm_fir_decimate_init_f32(&lo, NUM_TAPS, 2, (float32_t *)coeffs_lo_rev, state_lo, BLOCK_SIZE);
    arm_fir_decimate_init_f32(&hi, NUM_TAPS, 2, (float32_t *)coeffs_hi_rev, state_hi, BLOCK_SIZE);
}

// One analysis level on the M4F: two decimating FIR calls, one per branch.
void dwt_level_cmsis(const float *in, int n, float *cA, float *cD)
{
    arm_fir_decimate_f32(&lo, (float *)in, cA, n);   // low-pass + /2 -> n/2 out
    arm_fir_decimate_f32(&hi, (float *)in, cD, n);   // high-pass + /2 -> n/2 out
}

CMSIS does the boundary differently from the periodic wrap above: it carries a state buffer holding the tail of the previous block, which is the right behaviour for a continuous stream and is why the embedded result matches the desktop only up to the edge transient. With the single-cycle MAC and the FPU on the F446RE, the per-block cost is the few thousand MACs the budget above predicts, well under the block period.

On the ESP32-S3: ESP-DSP, or just the portable loop

The ESP32-S3 has no CMSIS-DSP, but it has ESP-DSP and a single-precision FPU, so two routes are open. The simplest is that the portable dwt_step above compiles and runs unchanged: it is ordinary float arithmetic, and the inner product over L taps is small. For more speed, ESP-DSP provides dsps_fir_f32 (a SIMD-accelerated FIR on the Xtensa LX7); run the low- and high-pass filters with it and decimate by taking every second output.

#include "esp_dsp.h"

// ESP-DSP's FIR instance type is fir_f32_t (the function is dsps_fir_f32), which
// reads taps in natural (time) order. It carries a delay line that must persist
// for the filter's lifetime. Reuse the natural-order db4 coeffs_lo / coeffs_hi
// from the CMSIS block above (NOT the time-reversed coeffs_*_rev that CMSIS needs).
static float delay_lo[NUM_TAPS] = {0};
static float delay_hi[NUM_TAPS] = {0};
fir_f32_t fir_lo, fir_hi;

void dwt_level_esp_init(void)   // call once at startup
{
    dsps_fir_init_f32(&fir_lo, (float *)coeffs_lo, delay_lo, NUM_TAPS);
    dsps_fir_init_f32(&fir_hi, (float *)coeffs_hi, delay_hi, NUM_TAPS);
}

// One analysis level on the ESP32-S3. dsps_fir_f32 gives the full-rate FIR
// output; keep every second sample to complete the downsample-by-2.
void dwt_level_esp(const float *in, int n, float *cA, float *cD,
                   float *full)   // full[] is an n-sample scratch
{
    dsps_fir_f32(&fir_lo, (float *)in, full, n);
    for (int m = 0; m < n / 2; m++) cA[m] = full[2 * m];
    dsps_fir_f32(&fir_hi, (float *)in, full, n);
    for (int m = 0; m < n / 2; m++) cD[m] = full[2 * m];
}

Filtering at full rate and throwing away half the outputs wastes work the decimating CMSIS call avoids; the polyphase trick from multirate recovers it by computing only the outputs you keep, and is worth doing if the DWT is the hot loop. For most block sizes the plain loop is already fast enough on the LX7.

The lifting scheme: integer wavelets with no FPU

Both routes above use floating point and produce non-integer coefficients. On an 8-bit part with no FPU, or whenever you need the transform to be losslessly reversible (image and data compression), the answer is the lifting scheme (Sweldens 1998). Lifting factors the filter bank into alternating predict and update steps that overwrite the array in place, halving the arithmetic and, crucially, staying exactly invertible even when each step rounds.

For the Haar wavelet the factorisation is a one-liner per pair. The detail is the odd sample minus a prediction from the even one; the approximation is the even sample updated by half the detail, which is just their (floored) average:

#include <stdint.h>

// In-place integer Haar via lifting (the reversible "S-transform"). After the
// call, even indices hold approximations, odd indices hold details. n even.
void haar_lift_fwd(int32_t *x, int n)
{
    for (int i = 0; i < n; i += 2) {
        x[i + 1] -= x[i];               // predict: detail = odd - even
        x[i]     += x[i + 1] >> 1;      // update:  approx = even + floor(detail/2)
    }
}

// Exact inverse: undo update, then undo predict. Bit-for-bit lossless, because
// the same floored detail/2 that was added back here was subtracted there.
void haar_lift_inv(int32_t *x, int n)
{
    for (int i = 0; i < n; i += 2) {
        x[i]     -= x[i + 1] >> 1;      // undo update
        x[i + 1] += x[i];               // undo predict
    }
}

The reversibility is the point and is worth seeing. The forward update adds floor(detail/2); the inverse subtracts the same quantity before the detail is touched, so whatever the floor discarded is restored exactly. No rounding error accumulates across levels, which a floating-point transform cannot promise. The cost is one subtract, one shift, and one add per pair, with no multiplies and no float unit, so the integer Haar lifting runs on an ATmega328P as happily as on a Cortex-M4F.

Code 
# verify the integer lifting is exactly reversible, the property the C relies on
def haar_lift_fwd(x):
    x = x.astype(np.int64).copy()
    x[1::2] -= x[0::2]                  # predict
    x[0::2] += x[1::2] >> 1             # update (arithmetic shift = floor for >>)
    return x

def haar_lift_inv(x):
    x = x.copy()
    x[0::2] -= x[1::2] >> 1
    x[1::2] += x[0::2]
    return x

rng = np.random.default_rng(0)
sig = rng.integers(-1000, 1000, size=64)
assert np.array_equal(haar_lift_inv(haar_lift_fwd(sig)), sig)   # bit-exact round trip
print("integer Haar lifting: forward then inverse recovers the input exactly")
integer Haar lifting: forward then inverse recovers the input exactly

Longer wavelets lift too. Daubechies-4 factors into two predict-update pairs and a scaling, and the same integer, in-place, reversible recipe applies; the predict and update coefficients are just longer. For a no-FPU build the Haar lifting above is usually enough for a coarse multiscale split (activity detection, edge finding, a cheap denoiser); reach for the longer wavelet only when its smoothness earns its keep.

When to use which

Target Route Why
NUCLEO-F446RE (M4F) arm_fir_decimate_f32, two instances per level FPU + single-cycle MAC + decimation built in; no wasted full-rate samples
ESP32-S3 (LX7) portable dwt_step, or dsps_fir_f32 + decimate FPU present; SIMD FIR if the DWT is the hot loop
ATmega328P / M0+ (no FPU) integer Haar lifting no float, in place, exactly reversible; the only sane option without an FPU
Lossless compression, any part lifting (Haar or db4) integer-to-integer and bit-exact invertible, which the float filter bank is not

The thread through all four is the one from the main page: the DWT is a filter bank, and a filter bank is something this workshop has already put on hardware. Wavelets do not need new embedded machinery, only the multirate machinery arranged as a tree.

References

Sweldens, Wim. 1998. “The Lifting Scheme: A Construction of Second Generation Wavelets.” SIAM Journal on Mathematical Analysis 29 (2): 511–46.