81а. Digital filters: FIR, IIR and CMSIS-DSP

In section 81 you suppressed sensor noise with a moving average and exponential smoothing — the filters that are the easiest to write. When you need precision — pass everything up to 100 Hz and suppress everything above 200 — you use real digital filters. There are two kinds: FIR (finite impulse response) and IIR (infinite impulse response). On STM32 you do not write them by hand: the CMSIS-DSP library has ready-made functions optimized for the Cortex-M core, for floating-point numbers and for fixed-point fractions.

All numbers in this section come from a run: the CMSIS-DSP 1.16.2 functions from the pack were built with gcc on a PC (_verify/dsp81a/dsp81a.c), and sine waves of four frequencies were passed through each filter, with the gain measured on the steady-state signal. Separately, without the library, the frequency response was computed from the coefficients in Python — both calculations agreed to the fifth digit. The check is _verify/dsp81a_check.py.

Terms used in the documentation

Term	Meaning
finite impulse response (FIR) filter	the output is a weighted sum of the inputs
infinite impulse response (IIR) filter	past outputs also take part in the sum
biquad, second-order section	a second-order stage
taps	the coefficients of a FIR filter
group delay	the delay of the signal envelope through the filter
feedforward, feedback coefficients	coefficients on the inputs and on the past outputs
postShift	the shift of the result after accumulation in Q15/Q31
cutoff frequency	the frequency where the passband ends

81а.1 Two kinds of filters

The moving average from section 81 is a FIR filter with all coefficients equal to 1/N, and exponential smoothing is a first-order IIR filter.

81а.2 FIR: design and check

Sampling rate 1000 Hz, cutoff 100 Hz, 31 coefficients. The window method: the ideal low-pass filter (the sinc function) is multiplied by a Hamming window, and the coefficients are normalized so that they sum to 1:

for (int n = 0; n < NTAPS; n++) {
    double m = n - (NTAPS - 1) / 2.0;
    double h = (m == 0) ? 2 * fc / FS : sin(2 * M_PI * fc / FS * m) / (M_PI * m);
    h *= 0.54 - 0.46 * cos(2 * M_PI * n / (NTAPS - 1));
    fir_c[n] = (float32_t)h;
    sum += h;
}

The coefficients are symmetric — the middle one is 0.19966, the edges are zero. A run through arm_fir_f32 and arm_fir_q15:

• At the cutoff frequency a filter designed with the window method attenuates by half (−6 dB), not by −3 dB as analog filters do — that is how the cutoff frequency is defined here.
• −70 dB at 200 Hz is a hit close to a zero of the frequency response. It is more honest to look at the side lobes: the first one is −51.4 dB at 166.7 Hz, and above 200 Hz none rises above −62.6 dB — suppression of more than 1000 times.
• Q15 is almost as good as float: the coefficients are 16-bit, but the function accumulates in a 64-bit accumulator.

Delay. An impulse at the input gives a peak at the output after 15 samples = 15 ms — exactly (N − 1) / 2. A symmetric FIR has the same delay for all frequencies: the waveform is preserved, it is just shifted. For a controller (section 47) 15 ms of lag in the loop is a lot; there you use a shorter filter or an IIR.

81а.3 IIR: one stage instead of 31 coefficients

A second-order stage (biquad) with a Butterworth response at 100 Hz, coefficients from the bilinear transform:

A run through arm_biquad_cascade_df1_f32:

Five multiplications instead of 31 — but the roll-off is gentler: −14 dB at 200 Hz versus −70 for the FIR. For a steeper one, chain several stages: that is why the function is called cascade, and the number of stages is set at initialization.

81а.4 The sign trap

This is the main pitfall of CMSIS-DSP. The library computes a biquad like this:

y[n] = b0x[n] + b1x[n−1] + b2x[n−2] + a1y[n−1] + a2y[n−2]

— with a plus in front of the feedback. Most filter design programs (and the formulas above) use a minus. So in the CMSIS coefficient array you put {b0, b1, b2, −a1, −a2}:

float32_t cm_ok[5] = {(float32_t)c[0], (float32_t)c[1], (float32_t)c[2], (float32_t)-c[3], (float32_t)-c[4]};

What happens if you put the coefficients in as they are: the run at all four frequencies gave infinity. The filter's denominator changed sign, and one of the poles ended up at −1.43 — outside the unit circle; the filter is unstable and the output grows without limit. With the correct sign the poles are complex, with radius 0.642 — inside the circle, and the filter is stable.

The stability check for a single stage is simple: |a2| < 1 and |a1| < 1 + a2 (in the minus notation).

81а.5 IIR in Q15: a coefficient greater than one

Q15 is a fraction from −1 to almost +1. But −a1 = 1.143 is greater than one and does not fit in Q15. The way out is to divide all coefficients by 2 and shift the accumulated result left by 1 bit: that is the postShift = 1 parameter at initialization.

The order of the coefficients in Q15 is also different: {b0, 0, b1, b2, a1, a2} — between b0 and b1 there is a zero word: this way, on Cortex-M4, the coefficients are fetched in pairs by 16-bit SIMD instructions.

1105, 0, 2210, 1105, 18727, -6763      postShift = 1

The result is almost the same as with float:

On a core without floating point (Cortex-M0, M3) Q15 is the only fast way; on M4 and M7 with an FPU, float is usually no slower and does not need this fuss.

81а.6 What it looks like in firmware

1. The coefficients are computed in advance — on a PC or once at startup — and stored in a const array.
2. The filter state is an array: for a FIR of length NTAPS + BLOCK − 1, for a biquad 4 numbers per stage. It lives for as long as the filter runs.
3. Samples are processed in blocks: DMA fills half of the ADC buffer (section 14), and the "half complete" interrupt converts the codes of that half into the filter format and calls arm_fir_f32 while DMA fills the other half. ADC codes are unsigned 12-bit: for float they are converted as (code − 2048) / 2048.0f, for Q15 as (int16_t)((code − 2048) * 16). You cannot pass the ADC buffer to the filter directly.
4. The block size is set at initialization and must match what you pass on every call.

81а.7 Pitfalls

The IIR "blew up" — the output goes to infinity. The feedback coefficients were put in without changing the sign. In CMSIS it is {b0, b1, b2, −a1, −a2}.

The Q15 biquad outputs garbage. The zero element of {b0, 0, b1, b2, a1, a2} was ignored, or a coefficient greater than 1 was not divided by 2 with postShift set.

The cutoff frequency is "wrong". A window-method FIR is at −6 dB at the cutoff frequency, a Butterworth filter at −3 dB. When you compare filters, say which level you mean.

The filter lags. A FIR with N coefficients delays by (N − 1)/2 samples. In a control loop this can make the system oscillate.

The first samples look strange. The filter state starts at zero — the output settles over a time equal to the FIR length (or several time constants for an IIR).

The state buffer is too small. For a FIR it is NTAPS + BLOCK − 1, not NTAPS: otherwise the function corrupts the adjacent memory.

81а.8 Summary

• FIR: always stable, linear phase, delay (N − 1)/2 — for a filter with 31 coefficients at 1 kHz this is 15 ms; −6 dB at the cutoff frequency, above 200 Hz suppression of at least 62.6 dB.
• IIR biquad: 5 multiplications, −3.01 dB at 100 Hz, a gentler roll-off; for a steeper one — a cascade of stages.
• In CMSIS the feedback has a "+" sign: the coefficients are −a1, −a2. Otherwise there is a pole at −1.43 and infinity at the output.
• Q15: coefficients greater than 1 — divide by 2 and set postShift = 1; a zero word after b0. Accuracy is the same as float to the fourth digit.

Full sources

dsp81a.c — FIR and IIR on CMSIS-DSP: float and Q15.

/* Раздел 81а: КИХ- и БИХ-фильтры на CMSIS-DSP 1.16.2, собранной gcc на компьютере.
 *
 * Частота отсчётов 1000 Гц, полоса пропускания до 100 Гц:
 *  1. КИХ 31 коэффициент, окно Хэмминга — arm_fir_f32 и arm_fir_q15;
 *  2. БИХ второго порядка (Баттерворт, 100 Гц) — arm_biquad_cascade_df1_f32 и _q15;
 *  3. ловушка знака: коэффициенты обратной связи в CMSIS со знаком «+».
 * Для каждого фильтра — усиление на 50, 100, 200 и 300 Гц по установившемуся синусу. */
#include <math.h>
#include <stdio.h>
#include <string.h>

#include "arm_math.h"

#define FS      1000.0
#define NTAPS   31
#define BLOCK   32
#define NSAMP   4096          /* 4096 отсчётов: 4,1 с — хватает на установку и замер */

static float32_t fir_c[NTAPS];

/* КИХ методом окна: идеальный ФНЧ (sinc) × окно Хэмминга, сумма коэффициентов = 1. */
static void design_fir(double fc)
{
    double sum = 0;
    for (int n = 0; n < NTAPS; n++) {
        double m = n - (NTAPS - 1) / 2.0;
        double h = (m == 0) ? 2 * fc / FS : sin(2 * M_PI * fc / FS * m) / (M_PI * m);
        h *= 0.54 - 0.46 * cos(2 * M_PI * n / (NTAPS - 1));
        fir_c[n] = (float32_t)h;
        sum += h;
    }
    for (int n = 0; n < NTAPS; n++)
        fir_c[n] = (float32_t)(fir_c[n] / sum);
}

/* Биквад ФНЧ (Баттерворт, Q = 1/√2) по билинейному преобразованию. Возвращает b0 b1 b2 a1 a2
 * в привычной записи y = b·x − a1·y[n−1] − a2·y[n−2]. */
static void design_biquad(double fc, double c[5])
{
    double w0 = 2 * M_PI * fc / FS, al = sin(w0) / (2 * M_SQRT1_2), cw = cos(w0);
    double a0 = 1 + al;
    c[0] = (1 - cw) / 2 / a0;
    c[1] = (1 - cw) / a0;
    c[2] = (1 - cw) / 2 / a0;
    c[3] = -2 * cw / a0;
    c[4] = (1 - al) / a0;
}

static float32_t xin[NSAMP], yout[NSAMP];
static q15_t xq[NSAMP], yq[NSAMP];

static void make_tone(double f, double amp)
{
    for (int i = 0; i < NSAMP; i++) {
        xin[i] = (float32_t)(amp * sin(2 * M_PI * f * i / FS));
        xq[i] = (q15_t)lrint(xin[i] * 32767);
    }
}

/* Амплитуда на последних 1000 отсчётах (ровно 50…300 целых периодов для наших частот). */
static double amp_f(const float32_t *y)
{
    double s = 0;
    for (int i = NSAMP - 1000; i < NSAMP; i++)
        s += (double)y[i] * y[i];
    return sqrt(2 * s / 1000);
}

static double amp_q(const q15_t *y)
{
    double s = 0;
    for (int i = NSAMP - 1000; i < NSAMP; i++)
        s += ((double)y[i] / 32768) * ((double)y[i] / 32768);
    return sqrt(2 * s / 1000);
}

static const double FREQ[] = {50, 100, 200, 300};

int main(void)
{
    design_fir(100);
    printf("fir_taps=%d", NTAPS);
    for (int n = 0; n < NTAPS; n++)
        printf("%s%.9f", n ? "," : "=", fir_c[n]);
    printf("\n");

    /* ── КИХ, float ── */
    static float32_t st_f[NTAPS + BLOCK - 1];
    static q15_t fir_q[NTAPS], st_q[NTAPS + BLOCK - 1];
    arm_float_to_q15(fir_c, fir_q, NTAPS);
    for (unsigned k = 0; k < 4; k++) {
        arm_fir_instance_f32 S;
        arm_fir_init_f32(&S, NTAPS, fir_c, st_f, BLOCK);
        make_tone(FREQ[k], 0.5);
        for (int i = 0; i < NSAMP; i += BLOCK)
            arm_fir_f32(&S, xin + i, yout + i, BLOCK);
        double gf = amp_f(yout) / 0.5;
        arm_fir_instance_q15 Q;
        arm_fir_init_q15(&Q, NTAPS, fir_q, st_q, BLOCK);
        for (int i = 0; i < NSAMP; i += BLOCK)
            arm_fir_q15(&Q, xq + i, yq + i, BLOCK);
        double gq = amp_q(yq) / 0.5;
        printf("fir f=%.0f gain_f32=%.5f db=%.2f gain_q15=%.5f\n", FREQ[k], gf, 20 * log10(gf), gq);
    }

    /* задержка КИХ: импульс на входе — пик на выходе */
    {
        arm_fir_instance_f32 S;
        arm_fir_init_f32(&S, NTAPS, fir_c, st_f, BLOCK);
        memset(xin, 0, sizeof xin);
        xin[0] = 1;
        arm_fir_f32(&S, xin, yout, BLOCK);
        arm_fir_f32(&S, xin + BLOCK, yout + BLOCK, BLOCK);
        int pk = 0;
        for (int i = 1; i < 2 * BLOCK; i++)
            if (yout[i] > yout[pk]) pk = i;
        printf("fir delay_samples=%d delay_ms=%.1f\n", pk, pk * 1000 / FS);
    }

    /* ── БИХ ── */
    double c[5];
    design_biquad(100, c);
    printf("biquad b0=%.9f b1=%.9f b2=%.9f a1=%.9f a2=%.9f\n", c[0], c[1], c[2], c[3], c[4]);
    /* CMSIS: y = b0x + b1x1 + b2x2 + a1y1 + a2y2 — обратная связь со знаком «+», поэтому −a1, −a2 */
    float32_t cm_ok[5] = {(float32_t)c[0], (float32_t)c[1], (float32_t)c[2], (float32_t)-c[3], (float32_t)-c[4]};
    float32_t cm_bad[5] = {(float32_t)c[0], (float32_t)c[1], (float32_t)c[2], (float32_t)c[3], (float32_t)c[4]};
    for (int variant = 0; variant < 2; variant++) {
        for (unsigned k = 0; k < 4; k++) {
            float32_t st[4] = {0};
            arm_biquad_casd_df1_inst_f32 S;
            arm_biquad_cascade_df1_init_f32(&S, 1, variant ? cm_bad : cm_ok, st);
            make_tone(FREQ[k], 0.5);
            for (int i = 0; i < NSAMP; i += BLOCK)
                arm_biquad_cascade_df1_f32(&S, xin + i, yout + i, BLOCK);
            double g = amp_f(yout) / 0.5;
            printf("biquad_%s f=%.0f gain=%.5f db=%.2f\n", variant ? "matlab_sign" : "cmsis_sign", FREQ[k], g, 20 * log10(g));
        }
    }

    /* ── БИХ в Q15: |a1| > 1 не помещается в Q15 — коэффициенты делят на 2, postShift = 1 ── */
    q15_t cq[6];
    {
        double half[5] = {c[0] / 2, c[1] / 2, c[2] / 2, -c[3] / 2, -c[4] / 2};
        cq[0] = (q15_t)lrint(half[0] * 32768);
        cq[1] = 0;                                  /* в Q15 между b0 и b1 — нулевое слово */
        cq[2] = (q15_t)lrint(half[1] * 32768);
        cq[3] = (q15_t)lrint(half[2] * 32768);
        cq[4] = (q15_t)lrint(half[3] * 32768);
        cq[5] = (q15_t)lrint(half[4] * 32768);
        printf("biquad_q15 coeffs=%d,%d,%d,%d,%d,%d postShift=1 (-a1=%.6f > 1)\n", cq[0], cq[1], cq[2], cq[3], cq[4], cq[5], -c[3]);
    }
    for (unsigned k = 0; k < 4; k++) {
        q15_t st[4] = {0};
        arm_biquad_casd_df1_inst_q15 S;
        arm_biquad_cascade_df1_init_q15(&S, 1, cq, st, 1);
        make_tone(FREQ[k], 0.5);
        for (int i = 0; i < NSAMP; i += BLOCK)
            arm_biquad_cascade_df1_q15(&S, xq + i, yq + i, BLOCK);
        double g = amp_q(yq) / 0.5;
        printf("biquad_q15 f=%.0f gain=%.5f db=%.2f\n", FREQ[k], g, 20 * log10(g));
    }
    printf("done=1\n");
    return 0;
}

Next: section 81б "The Goertzel algorithm: finding a tone without an FFT" — a single frequency without the whole spectrum: DTMF dialing, a window and q31.