121б. Audio on a microcontroller: from a PDM microphone to MFCC and command recognition

Sections 116–120 recognised vibration: a window of samples, a spectrum, a few features, a small model. With audio the chain is the same, but longer at both ends. At the start there is a digital microphone that delivers not samples but a stream of bits at a rate of megahertz. At the end there is not one answer per window but a stream of answers 50 times a second, from which you must pick the moment "the command was spoken" and not trigger on an unrelated sound.

The section goes the whole way in integers: PDM decimation with CIC and FIR filters, MFCC features, an int8 network and a decision on the stream. The code is our own, in C, without libraries; every stage is checked against an independent computation in numpy.

About the data, honestly. The section contains no recordings of real speech. The "commands" here are synthetic sounds: four classes (a rising tone, a falling tone, two short beeps, a vowel with two formants) with random pitch, duration, loudness and noise, and a fifth class "other": noise, unrelated tones, cut-off commands. The numbers in this section show that the chain works and that the integers reproduce the floating-point results. They say nothing about accuracy on real speech.

Prerequisites: PDM microphone, DFSDM and the Sinc filter — section 18б; FFT, window and leakage — section 117; FIR filters and Q15 — section 81а; a single bin without an FFT — section 81б; training — section 118; int8 quantization, the Q31 multiplier and shift — section 119. The DFSDM, MDF and ADF blocks by family — appendix H.

What it is called in the documentation

Term	Meaning
PDM (pulse density modulation)	pulse density modulation: one bit per clock
decimation	decimation: filter, then lower the sample rate
CIC (cascaded integrator-comb)	a cascade of integrators and combs; the same as Sinc in DFSDM
droop	roll-off of the frequency response inside the passband
frame, hop	the window and the step between windows
pre-emphasis	pre-emphasis: boosting the high frequencies
mel scale, mel filterbank	the mel scale, a set of triangular mel filters
MFCC (mel-frequency cepstral coefficients)	mel-frequency cepstral coefficients
DCT-II	discrete cosine transform, type II
block floating point	block normalisation: one common shift for the whole window
keyword spotting (KWS)	recognising short commands in a stream
posterior smoothing	smoothing the network outputs over time
false trigger	a false trigger
refractory period	a pause after a trigger

121б.1 The whole chain

Numbers of the chain:

• Microphone: 1 bit per clock at 3.072 MHz. Down to 16 kHz is a factor of 192.
• Decimation in two stages: CIC by 48 (to 64 kHz) and FIR by 4 (to 16 kHz).
• Window: 400 samples (25 ms), zero-padded to a 512-point FFT. The classic hop for speech is 10 ms: 98 windows per second of recording. Here the hop is 20 ms (320 samples): 49 windows per second, half the computation and memory.
• From each window, 10 MFCC coefficients. One second of sound is 49 × 10 = 490 int8 numbers: this is the network input.
• The 490 → 24 → 5 network answers after each new window, 50 times a second.

Where to compute it. The STM32F407 has no DFSDM block (the header stm32f407xx.h has not a single line about it): the microphone bits are received over I2S or SPI with DMA as ordinary bytes and decimated in software; that is the function pdm_decimate below. The STM32U585 (Cortex-M33 core) has the MDF1 block with six filters and ADF1 with one: the hardware computes the CIC, as DFSDM does in section 18б, and the software gets samples. Everything after 16 kHz does not depend on the family.

121б.2 From PDM to 16 kHz: CIC and FIR

Section 18б used a single Sinc filter straight down to the required rate. It has two drawbacks: roll-off in the upper part of the band and weak suppression of what folds back into the audio band during decimation. So here there are two stages.

CIC: four integrators at the PDM rate and four combs after decimation by 48. It has no multiplications, only additions. The output range is 48⁴ = 5308416, less than 2²³ = 8388608: the computation fits in 32 bits with headroom. The integrators overflow and wrap around, and that is normal: the combs subtract in the same modulo 2³² arithmetic, and the result is correct as long as the output range fits in the word.

FIR with 64 taps already runs at 64 kHz and keeps every 4th sample. It is a windowed sinc with an 8 kHz cutoff (Hamming window, section 81а). The CIC gain is removed in the coefficients: 48⁴ = 81 · 2¹⁶, so the taps are multiplied by 2²³ / 81 and the result is shifted by 23. The sum of the taps is 103560 (the exact 2²³ / 81 is 103563.1), the largest is 25264: they fit in int16.

/* Байты PDM старшим битом вперёд; бит 1 = +1, бит 0 = -1.
   Возвращает число готовых отсчётов 16 кГц, записанных в pcm. */
uint32_t pdm_decimate(pdm_t *st, const uint8_t *pdm, uint32_t nbytes, int16_t *pcm)
{
    uint32_t nout = 0;
    for (uint32_t i = 0; i < nbytes; i++) {
        uint32_t byte = pdm[i];
        for (int b = 7; b >= 0; b--) {
            st->i1 += ((byte >> b) & 1u) ? 1u : 0xFFFFFFFFu;   /* четыре интегратора */
            st->i2 += st->i1;
            st->i3 += st->i2;
            st->i4 += st->i3;
            if (++st->nbit < CIC_R)
                continue;
            st->nbit = 0;
            uint32_t d1 = st->i4 - st->c1;                       /* четыре гребёнки */
            st->c1 = st->i4;
            uint32_t d2 = d1 - st->c2;
            st->c2 = d1;
            uint32_t d3 = d2 - st->c3;
            st->c3 = d2;
            uint32_t d4 = d3 - st->c4;
            st->c4 = d3;
            st->line[st->pos] = (int32_t)d4;                     /* |d4| <= 48^4 < 2^23 */
            st->pos = (uint8_t)((st->pos + 1) & (FIR_N - 1));
            if (++st->phase < FIR_DEC)
                continue;
            st->phase = 0;
            int64_t acc = 0;
            for (uint32_t k = 0; k < FIR_N; k++)                 /* line[pos] — самый старый */
                acc += (int64_t)fir_q[k] * st->line[(st->pos + k) & (FIR_N - 1)];
            OPS(FIR_N);
            int32_t v = (int32_t)((acc + (1 << 22)) >> 23);
            pcm[nout++] = (int16_t)(v > 32767 ? 32767 : v < -32768 ? -32768 : v);
        }
    }
    return nout;
}

Frequency response of the chain from the formulas (CIC: a ratio of sines to the fourth power, FIR: a sum of the taps with twiddles):

The CIC roll-off at 8 kHz is only 0.90 dB: at 64 kHz its zeros are still far away. The FIR transition band runs from 7 to 10 kHz; whatever lies between 8 and 9 kHz is weakly suppressed and folds into the top kilohertz. For command recognition that is acceptable; for audio recording you take a longer filter.

Test run. There is no microphone; a model replaces it: a 2nd-order sigma-delta modulator of the same kind as in section 18б, but in integers. A tone with an amplitude of 0.45 of the modulator full scale, 768000 bits, 4000 samples at the output:

The output of the C function matched an integer model in numpy (convolutions instead of integrators and combs) in all 4000 samples of each tone. The largest magnitude at the CIC output is 2386140, within 48⁴. The signal/noise ratio here is a property of the modulator model: a real microphone has its own noise, and it is worse.

Cost. For 4 samples at 16 kHz the function performed 3072 additions in the integrators, 64 subtractions in the combs and 256 multiply-accumulates. Per second that is 12288000 additions, 256000 subtractions and 1024000 multiplications. The integrators are the main load: four additions for every bit from the microphone. This is an argument for a hardware filter (DFSDM, MDF, ADF) or for a lower microphone clock rate.

121б.3 MFCC: what each step does

MFCC is a way to compress the spectrum of a window into about ten numbers that describe its shape, not the loudness and not the pitch of the voice. The steps:

1. Pre-emphasis y[n] = x[n] − 0.97 · x[n−1]. Speech falls off towards high frequencies; the filter lifts the top and removes the DC component, so that the low bands do not drown out the rest.
2. Hamming window 0.54 − 0.46 · cos(2πn / 399), against leakage (section 117; a Hann window works just as well).
3. FFT of 512 points: a bin is 31.25 Hz, 257 bins up to 8 kHz.
4. Power of a bin: re² + im². The phase is discarded.
5. Mel filters. The ear resolves frequencies unevenly: finely at the bottom, coarsely at the top. The mel scale reproduces this: mel = 2595 · lg(1 + f / 700). 20 Hz is 31.75 mel, 1000 Hz is 999.99 mel, 8000 Hz is 2840.02 mel. On the mel scale, 42 points are placed evenly (step 68.49 mel): the edges and centres of 40 triangular filters. In hertz the first centres are 65.1, 113.1 and 164.0 Hz, the last ones 7004.2 and 7487.0 Hz. The energy of a band is the sum of the bin powers weighted by the triangle. There were 257 numbers; now there are 40.
6. Logarithm of the energy: loudness is perceived in decibels, and multiplying the signal by a number turns into a shift by a constant.
7. DCT-II over the 40 logarithms: c[i] = √(2/40) · Σ L[m] · cos(π · i · (m + 0.5) / 40). Neighbouring bands are strongly correlated; the cosine transform gathers the shape of the spectrum into the first coefficients. You take the first 10: c[0] is the overall loudness, the rest are ever finer details of the envelope. The fine comb of voice harmonics stays in the discarded ones.

The mel filter weights in the tables are computed by this fragment of the generator (gen_tables.py):

pts = [mel(FMIN) + (mel(FMAX) - mel(FMIN)) * i / (NMEL + 1) for i in range(NMEL + 2)]
first, band, wgt = None, [], []
for k in range(NFFT // 2 + 1):
    m = mel(k * FS / NFFT)
    if m < pts[0] or m >= pts[-1]:
        continue
    if first is None:
        first = k
    b = max(i for i in range(NMEL + 1) if pts[i] <= m)
    band.append(b)
    wgt.append(q15((m - pts[b]) / (pts[b + 1] - pts[b])))

The trick in it matters most for memory: the points are evenly spaced on the mel scale, so a bin between points b and b + 1 belongs to exactly two triangles: with weight w on the rising slope of one and 1 − w on the falling slope of its neighbour. One weight and one index per bin are enough: 255 bins take 765 bytes. The full matrix 40 × 257 in int16 would take 20560 bytes.

121б.4 Integers: spectrum, logarithm, tables

The whole window computation is one file. A radix-2 FFT, int32 data, Q15 twiddles from a sine table of 384 numbers:

/* Данные int32, повороты Q15. div2 = 0: без деления, спектр до 23 бит.
   div2 = 1: каждый каскад делит на 2 — спектр / 512 помещается в 16 бит. */
static void fft512(int32_t *re, int32_t *im, int div2)
{
    for (uint32_t i = 1, j = 0; i < NFFT; i++) {             /* перестановка по обратным битам */
        uint32_t bit = NFFT >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        if (i < j) {
            int32_t t = re[i]; re[i] = re[j]; re[j] = t;
            t = im[i]; im[i] = im[j]; im[j] = t;
        }
    }
    for (uint32_t len = 2, step = NFFT / 2; len <= NFFT; len <<= 1, step >>= 1) {
        uint32_t half = len >> 1;
        for (uint32_t k = 0; k < half; k++) {
            int32_t wr = sin_q15[k * step + NFFT / 4];       /* cos */
            int32_t wi = -sin_q15[k * step];
            for (uint32_t i = k; i < NFFT; i += len) {
                uint32_t j = i + half;
                int32_t tr = (int32_t)(((int64_t)wr * re[j] - (int64_t)wi * im[j]) >> 15);
                int32_t ti = (int32_t)(((int64_t)wr * im[j] + (int64_t)wi * re[j]) >> 15);
                re[j] = (re[i] - tr) >> div2;
                im[j] = (im[i] - ti) >> div2;
                re[i] = (re[i] + tr) >> div2;
                im[i] = (im[i] + ti) >> div2;
            }
        }
    }
    OPS(4 * (NFFT / 2) * 9);
}

The logarithm is computed without logf: the index of the top bit gives the integer part of log2, and the next six bits select a table node; between nodes it is a straight line segment.

/* log2(v) в Q12: целая часть — номер старшего бита, дробная — таблица из 65 узлов и отрезок между ними */
static int32_t log2_q12(uint64_t v)
{
    if (v == 0)
        return 0;
    int n = 63 - __builtin_clzll(v);
    uint32_t m = (uint32_t)((v << (63 - n)) >> 32);          /* 1.31: старший бит = 1 */
    uint32_t idx = (m >> 25) & 63u;
    uint32_t fr = (m >> 10) & 0x7FFFu;
    uint32_t a = log2_q15[idx], b = log2_q15[idx + 1];
    uint32_t l = a + (((b - a) * fr) >> 15);
    OPS(1);
    return (n << 12) + (int32_t)((l + 4) >> 3);
}

And the window itself, from samples to 10 coefficients:

/* y[n] = x[n] - 0,97·x[n-1] в Q15: без округления, до 31 бита */
static inline int32_t preemph(int32_t x, int32_t prev)
{
    return x * 32768 - PREEMPH_Q15 * prev;
}

/* x — 400 отсчётов окна, prev — отсчёт перед окном; cep — 10 коэффициентов */
void mfcc_q15(mfcc_q15_t *w, const int16_t *x, int16_t prev, int32_t *cep)
{
    /* 1. наибольший модуль после предыскажения и сдвиг, при котором он занимает 15 бит */
    int32_t peak = 0, p = prev;
    for (uint32_t n = 0; n < FRAME; n++) {
        int32_t y = preemph(x[n], p);
        p = x[n];
        if (y < 0)
            y = -y;
        if (y > peak)
            peak = y;
    }
    int sh = -1;                                             /* без нормировки: 17 бит -> 16 */
    if (!w->no_norm)
        sh = __builtin_clz((uint32_t)peak | 1u) - 2;         /* от -1 до 29 */
    w->shift = (int8_t)sh;
    /* 2. предыскажение, окно Хэмминга (в таблице половина окна) и нормировка одним сдвигом */
    p = prev;
    for (uint32_t n = 0; n < FRAME; n++) {
        int64_t y = (int64_t)preemph(x[n], p) * hamming_q15[n < FRAME / 2 ? n : FRAME - 1 - n];
        p = x[n];
        w->re[n] = (int32_t)((y + ((int64_t)1 << (29 - sh))) >> (30 - sh));
        w->im[n] = 0;
    }
    for (uint32_t n = FRAME; n < NFFT; n++)
        w->re[n] = w->im[n] = 0;
    OPS(3 * FRAME);
    /* 3. БПФ: сумма 400 отсчётов по 15 бит с окном — меньше 2^23 */
    fft512(w->re, w->im, w->fft16);
    /* 4. мощность бина и раскладка по двум соседним мел-полосам */
    for (uint32_t m = 0; m < NMEL; m++)
        w->e[m] = 0;
    for (uint32_t i = 0; i < MEL_BINS; i++) {
        uint32_t k = MEL_FIRST + i;
        uint64_t pw = (uint64_t)((int64_t)w->re[k] * w->re[k] + (int64_t)w->im[k] * w->im[k]);   /* < 2^46 */
        uint64_t up = (pw * (uint32_t)mel_w_q15[i]) >> 15;
        uint32_t b = mel_band[i];
        if (b < NMEL)
            w->e[b] += up;                                   /* подъём полосы b */
        if (b > 0)
            w->e[b - 1] += pw - up;                          /* спад полосы b - 1 */
    }
    OPS(3 * MEL_BINS);
    /* 5. логарифм; поправка возвращает сдвиг нормировки (и деление на 512 при fft16) */
    for (uint32_t m = 0; m < NMEL; m++) {
        int32_t l = log2_q12(w->e[m]) + ((w->fft16 ? 18 : 0) - 2 * sh) * 4096;
        w->logmel[m] = l < 0 ? 0 : l;
    }
    /* 6. ДКП-II: первые 10 коэффициентов */
    for (uint32_t i = 0; i < NCEP; i++) {
        int64_t acc = 0;
        for (uint32_t m = 0; m < NMEL; m++)
            acc += (int64_t)w->logmel[m] * dct_q15[i * NMEL + m];
        cep[i] = (int32_t)((acc + 16384) >> 15);
    }
    OPS(NCEP * NMEL);
}

Three decisions here are not obvious, and each one was checked by a test run.

Window normalisation. A quiet sound occupies the low bits, and what follows are multiplications with a shift by 15, so the low bits are lost. That is why, before the window, the largest magnitude is found and the whole window is shifted so that it occupies 15 bits. The shift sh is then returned in the logarithm by subtracting 2 · sh; this is the benefit of step 6: the scale has become an additive term.

Pre-emphasis without rounding. The product 0.97 · x is computed in Q15 and is not rounded to a whole sample: the difference takes up to 31 bits and is brought to 15 by a single shift together with the window. If you round earlier, on a quiet sound the rounding error is comparable to the signal itself.

Spectrum in int32. The usual trick for Q15 is to divide by 2 in every FFT stage so that the spectrum stays in 16 bits. That avoids overflow, but weak bins next to strong ones drown in rounding: after dividing by 512 only units are left of them. Here the spectrum grows freely: the sum of 400 samples of 15 bits with the window is less than 2²³, the power of a bin is less than 2⁴⁶, and the band energy accumulates in 64 bits.

Cross-check: 1200 recordings of 49 windows each, 58800 windows; the coefficients from C were converted to the natural logarithm and compared with numpy (float64, np.fft.rfft, the full filter matrix). For scale: the spread of the coefficients over the dataset is 18.46 for c[0], 4.65 for c[1], 1.47 for c[9], and the smallest is 1.07.

And the same 50 recordings at four loudness levels, three ways of computing. MFCC error, RMS and largest:

Without normalisation the error is larger than the coefficients themselves even at full loudness. Normalisation with a 16-bit spectrum makes it the same at all levels, but an RMS of about 0.5 is half the spread of the higher-order coefficients. The int32 spectrum gives an error tens of times smaller than the spread and does not depend on loudness. The price is 2048 bytes of RAM: the spectrum arrays take 4096 bytes instead of 2048.

Tables (all const, in flash): half of the Hamming window, 400 bytes; sine, 768; mel filters, 765; logarithm, 130; DCT, 800. In total 2863 bytes.

Computation per window: 11621 multiplications: 1200 for pre-emphasis, the search for the largest magnitude and the window, 9216 for the FFT (2304 butterflies of 4), 765 for power and mel bands, 40 for the logarithm, 400 for the DCT. The FFT is four fifths of all the computation; a real-valued input lets you compute it at half the length, which is not done here.

121б.5 An int8 network on top of the features

Feature normalisation. c[0] varies tens of times more than c[9], so each coefficient has its own mean and its own multiplier; they are computed from the training recordings and go into the firmware together with the weights (section 118: normalisation is part of the model). A feature is (c − mean) / deviation · 32, saturated at ±127: the network input is int8 right away, with scale 1/32 and zero point 0. Saturated features in training: 583 out of 392000.

The network is fully connected, 490 → 24 → 5, with ReLU. It was trained in numpy (Adam, 60 epochs, fixed seed) on features from the same mfcc_q15 function that runs on the device: the model sees exactly the numbers it will get. Dataset: 800 recordings for training and 400 held out, 160 and 80 per class.

Quantization follows section 119: int8 weights with one scale per layer, int32 biases, rescaling of the sum with a Q31 multiplier and a shift with a single rounding. The multipliers came out as 1882415126 with a shift of 41 and 1468209255 with a shift of 39. The largest sum in the first layer is 294014: int32 is more than enough.

/* нормировка: q = (c - среднее) / отклонение · 32, с насыщением; среднее и множитель — из обучения */
void kws_features(const int32_t *cep, int8_t *q)
{
    for (uint32_t i = 0; i < NCEP; i++) {
        int64_t v = ((int64_t)(cep[i] - kws_mu[i]) * kws_k[i] + (1 << 23)) >> 24;
        q[i] = (int8_t)(v > 127 ? 127 : v < -127 ? -127 : v);
    }
    OPS(NCEP);
}

/* пересчёт суммы в масштаб выхода: умножение на mult / 2^(31 + sh), одно округление */
static inline int32_t requant(int32_t acc, int32_t mult, int sh)
{
    return (int32_t)(((int64_t)acc * mult + ((int64_t)1 << (30 + sh))) >> (31 + sh));
}

void kws_infer(const int8_t *feat, int8_t *out)
{
    int8_t h[KWS_HID];
    for (uint32_t j = 0; j < KWS_HID; j++) {                 /* слой 1 с ReLU */
        const int8_t *w = &kws_w1[j * KWS_IN];
        int32_t acc = kws_b1[j];
        for (uint32_t i = 0; i < KWS_IN; i++)
            acc += w[i] * feat[i];                           /* нулевая точка входа 0 */
        int32_t v = requant(acc, KWS_M1, KWS_S1) + KWS_ZH;
        h[j] = (int8_t)(v < KWS_ZH ? KWS_ZH : v > 127 ? 127 : v);
    }
    for (uint32_t j = 0; j < KWS_OUT; j++) {                 /* слой 2 */
        const int8_t *w = &kws_w2[j * KWS_HID];
        int32_t acc = kws_b2[j];
        for (uint32_t i = 0; i < KWS_HID; i++)
            acc += w[i] * (h[i] - KWS_ZH);
        int32_t v = requant(acc, KWS_M2, KWS_S2);
        out[j] = (int8_t)(v < -128 ? -128 : v > 127 ? 127 : v);
    }
    OPS(KWS_HID * KWS_IN + KWS_OUT * KWS_HID + KWS_HID + KWS_OUT);
}

Result on 400 held-out recordings:

Confusion matrix of the int8 network in C (row: what it was, column: the answer):

Quantization did not change a single answer. The errors are mostly between "other" and the commands: "other" contains cut-off commands, and the boundary between "half a command is visible" and "a command" is blurred.

Size. 11909 parameters: weights 11880 bytes, biases 116, normalisation 80; 12076 bytes in total. In float, the weights and biases alone take 47636 bytes. One network answer is 11909 multiplications: 11880 multiply-accumulates and 29 in the rescaling. Almost all the weights are in the first layer, 490 × 24. A convolutional network with depthwise separable convolutions on the same 49 × 10 features needs fewer weights at the cost of more computation; it is the usual choice for real speech; here it was not trained.

Through the microphone. The first 20 held-out recordings of each class were passed through the whole chain: the microphone model, pdm_decimate, MFCC, the network.

The answer matched the direct one in 88 recordings. Almost all the losses are in the class "two": 9 correct out of 20. The chain changes the spectrum (the microphone model stretches the samples by linear interpolation, and the filter rolls off towards 8 kHz), and the network never saw this during training. The practical conclusion: train on recordings that went through the same microphone and the same filters as in the device.

121б.6 The stream: sliding window, smoothing, false triggers

On the device there are no recordings, only a stream. The function kws_push accepts samples as they arrive (even 16 at a time from a DMA handler), assembles a window, after each window shifts the feature buffer by one frame and runs the network:

/* Отсчёты 16 кГц по мере прихода. Возвращает число новых кадров; после каждого
   кадра признаки сдвинуты и k->out пересчитан. */
int kws_push(kws_t *k, const int16_t *pcm, uint32_t n)
{
    int frames = 0;
    for (uint32_t i = 0; i < n; i++) {
        k->buf[k->fill++] = pcm[i];
        if (k->fill < 1 + FRAME)
            continue;
        int32_t cep[NCEP];
        mfcc_q15(&k->w, &k->buf[1], k->buf[0], cep);
        memmove(k->feat, k->feat + NCEP, NFEAT - NCEP);      /* самый старый кадр уходит */
        kws_features(cep, &k->feat[NFEAT - NCEP]);
        kws_infer(k->feat, k->out);
        memmove(k->buf, k->buf + HOP, (1 + FRAME - HOP) * sizeof k->buf[0]);
        k->fill = 1 + FRAME - HOP;                           /* 80 отсчётов перекрытия и один перед окном */
        frames++;
    }
    return frames;
}

The network looks at the last second and answers every 20 ms. While a command is sliding into the window, the network sees a piece of it and answers anything; once the command is entirely inside the window, the answer is correct for many frames in a row. Hence three techniques against false triggers:

• smoothing: add up the network outputs over the last few frames;
• margin: trigger only if the best class beats the second one by at least a threshold;
• pause: after a trigger stay silent until the command leaves the window.

/* Сглаживание и решение. smooth — сколько последних ответов сети усреднять (до 16),
   margin — на сколько сумма лучшего класса должна обгонять второй (в единицах выхода, на один ответ),
   hold — сколько кадров молчать после срабатывания. Возвращает номер команды 1…4 или 0. */
int kws_decide(kws_post_t *p, const int8_t *out, int smooth, int margin, int hold)
{
    for (uint32_t c = 0; c < NCLASS; c++) {
        if (p->fill == smooth)
            p->sum[c] = (int16_t)(p->sum[c] - p->hist[p->pos][c]);
        p->hist[p->pos][c] = out[c];
        p->sum[c] = (int16_t)(p->sum[c] + out[c]);
    }
    p->pos = (uint8_t)((p->pos + 1) % smooth);
    if (p->fill < smooth)
        p->fill++;
    if (p->hold) {
        p->hold--;
        return 0;
    }
    if (p->fill < smooth)
        return 0;
    int best = 0, second = -32768;
    for (int c = 1; c < NCLASS; c++)
        if (p->sum[c] > p->sum[best])
            best = c;
    for (int c = 0; c < NCLASS; c++)
        if (c != best && p->sum[c] > second)
            second = p->sum[c];
    if (best == 0 || p->sum[best] - second < margin * smooth)
        return 0;
    p->hold = (uint16_t)hold;
    return best;
}

Test run: a 180 s stream, noise at −40 dB, 41 commands and 26 unrelated sounds (tones, noise bursts, a vowel with other formants) every 1.6…3 s, 8999 frames. The pause is 50 frames (1 s). A trigger counts as correct if it fell on a command of the same class or within 1.5 s after it.

What you can see:

• Without a margin the detector is useless: 24–26 commands out of 41 and as many false ones. It triggers on the very first piece of a command, names the wrong class (16–19 times), and with the pause blocks its own correct answer.
• The margin matters more than smoothing. Already with a margin of 40 and no smoothing: 41 out of 41 and 2 false.
• Smoothing lets you lower the margin: 8 frames and a margin of 20 give 41 out of 41 and 3 false, all on unrelated sounds.
• A large margin removes false triggers entirely but loses commands: with a margin of 60 and 8 frames, 21 out of 41 and not a single false one.

The threshold is chosen by the cost of an error. Here it was tuned on the same stream on which it was measured; on a different stream the numbers will be different.

Build for Cortex-M. The whole chain is one function per block of 384 bytes of PDM (1 ms):

/* pdm — 384 байта: 3072 бита, 16 отсчётов 16 кГц (1 мс) */
int audio_block(const uint8_t *pdm)
{
    int16_t pcm[16];
    uint32_t n = pdm_decimate(&pdm_state, pdm, 384, pcm);
    if (!kws_push(&kws_state, pcm, n))
        return 0;
    return kws_decide(&post_state, kws_state.out, 8, 40, 50);
}

arm-none-eabi-gcc 12.2.1, object file sizes in bytes:

Code is about 2 KB, data 15 KB, and almost all of it is the first-layer weights. RAM: the recognition state is 5888 bytes (of which 4584 are the MFCC working arrays, 802 the sample window, 490 the features), decimation 292, smoothing 94. The stack of each function is at most 80 bytes. An image without the startup file, with memmove and memset from the library: 17964 bytes of code and constants, 6276 bytes of RAM.

Time. One second of the stream takes 1176869 multiplications: 50 MFCC windows and 50 network answers. Plus PDM decimation from subsection 121б.2 if it is done in software. Cycles were not measured: the FFT multiplications are 64-bit products, the network ones are 8-bit multiply-accumulates, and they cost differently. If the computation is too much, the first step is to run the network not on every frame but on every second to fourth one.

121б.7 Pitfalls

A quiet command is not heard. There is no window normalisation: in Q15 the low bits are lost in the multiplications. Without it the MFCC error at −40 dB is 20.89, while the coefficients themselves spread from 1.07 to 18.46.

The features are noisy although there is no overflow. A 16-bit spectrum with division by 2 in every FFT stage: weak bins drown in rounding. The MFCC error is about 0.5 against 0.015 with an int32 spectrum.

Pre-emphasis rounded to a sample. On a quiet signal the rounding error of 0.97 · x is comparable to the difference between neighbouring samples. Compute the difference in Q15 and shift it once, together with the window.

Everything is correct on recordings, but not from the microphone. The network was trained on samples but runs after the modulator and the filters: 87 correct out of 100 instead of 97. Train on recordings passed through the same chain.

Triggers on the start of a command and with the wrong class. There is no margin: 19 wrong commands out of 27 false triggers. You need a threshold between the best and the second class and a pause after a trigger.

A threshold picked "by eye" loses commands. A margin of 60 with smoothing of 8 gives 21 commands out of 41. Tune the threshold with a run over a stream with known labels.

Feature normalisation differs from training. The mean and the multiplier are part of the model: compute them from the training recordings and put them in the same header as the weights.

The CIC integrators "overflow". That is how it should be: the subtractions in the combs use the same modulus. The only thing not allowed is for the output range (48⁴ for decimation by 48 and 4th order) to exceed the word. Do not saturate the integrators.

The mel filter matrix ate the flash. The full 40 × 257 matrix is 20560 bytes. A weight and a band index per bin take 765 bytes, with the same result.

The top of the spectrum is dirty after decimation. The FIR transition band: at 9 kHz the suppression is only 21.70 dB, and that folds into 7 kHz. Either use a longer filter or do not use the upper mel bands.

121б.8 Summary

• Chain: 3.072 MHz → 4th-order CIC by 48 → 64-tap FIR by 4 → 16 kHz; CIC range 48⁴ = 5308416, less than 2²³; gain at 1, 4 and 7 kHz is 0.02, −0.22 and −1.40 dB, as the formula predicts; the C output matched numpy bit for bit.
• Software decimation costs 12288000 additions and 1024000 multiplications per second; on the F407 it is done in software, on the U585 the CIC is computed by MDF and ADF.
• MFCC: window 400, hop 320, FFT 512, 40 mel bands 20…8000 Hz, log2 from a table, DCT-II, 10 coefficients; 11621 multiplications per window, tables 2863 bytes.
• Integers reproduce float if the window is normalised, pre-emphasis is not rounded and the spectrum is int32: MFCC error 0.0148 with a spread from 1.07; with a 16-bit spectrum about 0.5, without normalisation from 5 to 29.
• The 490 → 24 → 5 network in int8: 11909 parameters, 12076 bytes; on 400 synthetic recordings 389 correct both in float and in int8; C matched the Python computation bit for bit.
• Through the microphone model, 87 out of 100 instead of 97: train on recordings that went through the chain.
• 180 s stream: smoothing over 8 frames, margin 20, pause 1 s gives 41 commands out of 41 and 3 false; without a margin, 25 out of 41 and 24 false.
• Cortex-M4: code about 2 KB, tables and weights 15 KB, RAM 5888 + 292 + 94 bytes; 1176869 multiplications per second of the stream.
• The commands are synthetic: this is a test of the chain, not accuracy on speech.

Full sources

pdm.c: PDM decimation, a 4th-order CIC and a FIR filter.

/* Прореживание PDM: CIC 4-го порядка на 48 и КИХ на 64 отвода с прореживанием на 4. */
#include "audio.h"
#include "mfcc_tables.h"

_Static_assert(CIC_R * CIC_R * CIC_R * CIC_R == 81 * 65536, "48^4 = 81 * 2^16: усиление CIC");
_Static_assert((FIR_N & (FIR_N - 1)) == 0, "длина КИХ — степень двойки");

/* Байты PDM старшим битом вперёд; бит 1 = +1, бит 0 = -1.
   Возвращает число готовых отсчётов 16 кГц, записанных в pcm. */
uint32_t pdm_decimate(pdm_t *st, const uint8_t *pdm, uint32_t nbytes, int16_t *pcm)
{
    uint32_t nout = 0;
    for (uint32_t i = 0; i < nbytes; i++) {
        uint32_t byte = pdm[i];
        for (int b = 7; b >= 0; b--) {
            st->i1 += ((byte >> b) & 1u) ? 1u : 0xFFFFFFFFu;   /* четыре интегратора */
            st->i2 += st->i1;
            st->i3 += st->i2;
            st->i4 += st->i3;
            if (++st->nbit < CIC_R)
                continue;
            st->nbit = 0;
            uint32_t d1 = st->i4 - st->c1;                       /* четыре гребёнки */
            st->c1 = st->i4;
            uint32_t d2 = d1 - st->c2;
            st->c2 = d1;
            uint32_t d3 = d2 - st->c3;
            st->c3 = d2;
            uint32_t d4 = d3 - st->c4;
            st->c4 = d3;
            st->line[st->pos] = (int32_t)d4;                     /* |d4| <= 48^4 < 2^23 */
            st->pos = (uint8_t)((st->pos + 1) & (FIR_N - 1));
            if (++st->phase < FIR_DEC)
                continue;
            st->phase = 0;
            int64_t acc = 0;
            for (uint32_t k = 0; k < FIR_N; k++)                 /* line[pos] — самый старый */
                acc += (int64_t)fir_q[k] * st->line[(st->pos + k) & (FIR_N - 1)];
            OPS(FIR_N);
            int32_t v = (int32_t)((acc + (1 << 22)) >> 23);
            pcm[nout++] = (int16_t)(v > 32767 ? 32767 : v < -32768 ? -32768 : v);
        }
    }
    return nout;
}

mfcc_q15.c: MFCC in integers, a window of 400 samples → 10 coefficients.

/* MFCC в целых числах: окно 400 отсчётов int16 -> 10 коэффициентов (log2 в Q12). */
#include "audio.h"
#include "mfcc_tables.h"

#define PREEMPH_Q15 31785            /* 0,97 */

/* Данные int32, повороты Q15. div2 = 0: без деления, спектр до 23 бит.
   div2 = 1: каждый каскад делит на 2 — спектр / 512 помещается в 16 бит. */
static void fft512(int32_t *re, int32_t *im, int div2)
{
    for (uint32_t i = 1, j = 0; i < NFFT; i++) {             /* перестановка по обратным битам */
        uint32_t bit = NFFT >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        if (i < j) {
            int32_t t = re[i]; re[i] = re[j]; re[j] = t;
            t = im[i]; im[i] = im[j]; im[j] = t;
        }
    }
    for (uint32_t len = 2, step = NFFT / 2; len <= NFFT; len <<= 1, step >>= 1) {
        uint32_t half = len >> 1;
        for (uint32_t k = 0; k < half; k++) {
            int32_t wr = sin_q15[k * step + NFFT / 4];       /* cos */
            int32_t wi = -sin_q15[k * step];
            for (uint32_t i = k; i < NFFT; i += len) {
                uint32_t j = i + half;
                int32_t tr = (int32_t)(((int64_t)wr * re[j] - (int64_t)wi * im[j]) >> 15);
                int32_t ti = (int32_t)(((int64_t)wr * im[j] + (int64_t)wi * re[j]) >> 15);
                re[j] = (re[i] - tr) >> div2;
                im[j] = (im[i] - ti) >> div2;
                re[i] = (re[i] + tr) >> div2;
                im[i] = (im[i] + ti) >> div2;
            }
        }
    }
    OPS(4 * (NFFT / 2) * 9);
}

/* log2(v) в Q12: целая часть — номер старшего бита, дробная — таблица из 65 узлов и отрезок между ними */
static int32_t log2_q12(uint64_t v)
{
    if (v == 0)
        return 0;
    int n = 63 - __builtin_clzll(v);
    uint32_t m = (uint32_t)((v << (63 - n)) >> 32);          /* 1.31: старший бит = 1 */
    uint32_t idx = (m >> 25) & 63u;
    uint32_t fr = (m >> 10) & 0x7FFFu;
    uint32_t a = log2_q15[idx], b = log2_q15[idx + 1];
    uint32_t l = a + (((b - a) * fr) >> 15);
    OPS(1);
    return (n << 12) + (int32_t)((l + 4) >> 3);
}

/* y[n] = x[n] - 0,97·x[n-1] в Q15: без округления, до 31 бита */
static inline int32_t preemph(int32_t x, int32_t prev)
{
    return x * 32768 - PREEMPH_Q15 * prev;
}

/* x — 400 отсчётов окна, prev — отсчёт перед окном; cep — 10 коэффициентов */
void mfcc_q15(mfcc_q15_t *w, const int16_t *x, int16_t prev, int32_t *cep)
{
    /* 1. наибольший модуль после предыскажения и сдвиг, при котором он занимает 15 бит */
    int32_t peak = 0, p = prev;
    for (uint32_t n = 0; n < FRAME; n++) {
        int32_t y = preemph(x[n], p);
        p = x[n];
        if (y < 0)
            y = -y;
        if (y > peak)
            peak = y;
    }
    int sh = -1;                                             /* без нормировки: 17 бит -> 16 */
    if (!w->no_norm)
        sh = __builtin_clz((uint32_t)peak | 1u) - 2;         /* от -1 до 29 */
    w->shift = (int8_t)sh;
    /* 2. предыскажение, окно Хэмминга (в таблице половина окна) и нормировка одним сдвигом */
    p = prev;
    for (uint32_t n = 0; n < FRAME; n++) {
        int64_t y = (int64_t)preemph(x[n], p) * hamming_q15[n < FRAME / 2 ? n : FRAME - 1 - n];
        p = x[n];
        w->re[n] = (int32_t)((y + ((int64_t)1 << (29 - sh))) >> (30 - sh));
        w->im[n] = 0;
    }
    for (uint32_t n = FRAME; n < NFFT; n++)
        w->re[n] = w->im[n] = 0;
    OPS(3 * FRAME);
    /* 3. БПФ: сумма 400 отсчётов по 15 бит с окном — меньше 2^23 */
    fft512(w->re, w->im, w->fft16);
    /* 4. мощность бина и раскладка по двум соседним мел-полосам */
    for (uint32_t m = 0; m < NMEL; m++)
        w->e[m] = 0;
    for (uint32_t i = 0; i < MEL_BINS; i++) {
        uint32_t k = MEL_FIRST + i;
        uint64_t pw = (uint64_t)((int64_t)w->re[k] * w->re[k] + (int64_t)w->im[k] * w->im[k]);   /* < 2^46 */
        uint64_t up = (pw * (uint32_t)mel_w_q15[i]) >> 15;
        uint32_t b = mel_band[i];
        if (b < NMEL)
            w->e[b] += up;                                   /* подъём полосы b */
        if (b > 0)
            w->e[b - 1] += pw - up;                          /* спад полосы b - 1 */
    }
    OPS(3 * MEL_BINS);
    /* 5. логарифм; поправка возвращает сдвиг нормировки (и деление на 512 при fft16) */
    for (uint32_t m = 0; m < NMEL; m++) {
        int32_t l = log2_q12(w->e[m]) + ((w->fft16 ? 18 : 0) - 2 * sh) * 4096;
        w->logmel[m] = l < 0 ? 0 : l;
    }
    /* 6. ДКП-II: первые 10 коэффициентов */
    for (uint32_t i = 0; i < NCEP; i++) {
        int64_t acc = 0;
        for (uint32_t m = 0; m < NMEL; m++)
            acc += (int64_t)w->logmel[m] * dct_q15[i * NMEL + m];
        cep[i] = (int32_t)((acc + 16384) >> 15);
    }
    OPS(NCEP * NMEL);
}

kws.c: command recognition, a stream of samples → int8 features → network → decision.

/* Распознавание команды: поток отсчётов -> признаки int8 -> сеть 490 -> 24 -> 5 -> решение. */
#include <string.h>
#include "audio.h"
#include "kws_model.h"

_Static_assert(KWS_IN == NFEAT && KWS_OUT == NCLASS, "размеры сети и признаков");

void kws_init(kws_t *k)
{
    memset(k, 0, sizeof *k);
    k->fill = 1;                                   /* buf[0] = 0: перед первым окном тишина; */
                                                   /* признаки 0 — среднее по обучающим записям */
}

/* нормировка: q = (c - среднее) / отклонение · 32, с насыщением; среднее и множитель — из обучения */
void kws_features(const int32_t *cep, int8_t *q)
{
    for (uint32_t i = 0; i < NCEP; i++) {
        int64_t v = ((int64_t)(cep[i] - kws_mu[i]) * kws_k[i] + (1 << 23)) >> 24;
        q[i] = (int8_t)(v > 127 ? 127 : v < -127 ? -127 : v);
    }
    OPS(NCEP);
}

/* пересчёт суммы в масштаб выхода: умножение на mult / 2^(31 + sh), одно округление */
static inline int32_t requant(int32_t acc, int32_t mult, int sh)
{
    return (int32_t)(((int64_t)acc * mult + ((int64_t)1 << (30 + sh))) >> (31 + sh));
}

void kws_infer(const int8_t *feat, int8_t *out)
{
    int8_t h[KWS_HID];
    for (uint32_t j = 0; j < KWS_HID; j++) {                 /* слой 1 с ReLU */
        const int8_t *w = &kws_w1[j * KWS_IN];
        int32_t acc = kws_b1[j];
        for (uint32_t i = 0; i < KWS_IN; i++)
            acc += w[i] * feat[i];                           /* нулевая точка входа 0 */
        int32_t v = requant(acc, KWS_M1, KWS_S1) + KWS_ZH;
        h[j] = (int8_t)(v < KWS_ZH ? KWS_ZH : v > 127 ? 127 : v);
    }
    for (uint32_t j = 0; j < KWS_OUT; j++) {                 /* слой 2 */
        const int8_t *w = &kws_w2[j * KWS_HID];
        int32_t acc = kws_b2[j];
        for (uint32_t i = 0; i < KWS_HID; i++)
            acc += w[i] * (h[i] - KWS_ZH);
        int32_t v = requant(acc, KWS_M2, KWS_S2);
        out[j] = (int8_t)(v < -128 ? -128 : v > 127 ? 127 : v);
    }
    OPS(KWS_HID * KWS_IN + KWS_OUT * KWS_HID + KWS_HID + KWS_OUT);
}

/* Отсчёты 16 кГц по мере прихода. Возвращает число новых кадров; после каждого
   кадра признаки сдвинуты и k->out пересчитан. */
int kws_push(kws_t *k, const int16_t *pcm, uint32_t n)
{
    int frames = 0;
    for (uint32_t i = 0; i < n; i++) {
        k->buf[k->fill++] = pcm[i];
        if (k->fill < 1 + FRAME)
            continue;
        int32_t cep[NCEP];
        mfcc_q15(&k->w, &k->buf[1], k->buf[0], cep);
        memmove(k->feat, k->feat + NCEP, NFEAT - NCEP);      /* самый старый кадр уходит */
        kws_features(cep, &k->feat[NFEAT - NCEP]);
        kws_infer(k->feat, k->out);
        memmove(k->buf, k->buf + HOP, (1 + FRAME - HOP) * sizeof k->buf[0]);
        k->fill = 1 + FRAME - HOP;                           /* 80 отсчётов перекрытия и один перед окном */
        frames++;
    }
    return frames;
}

/* Сглаживание и решение. smooth — сколько последних ответов сети усреднять (до 16),
   margin — на сколько сумма лучшего класса должна обгонять второй (в единицах выхода, на один ответ),
   hold — сколько кадров молчать после срабатывания. Возвращает номер команды 1…4 или 0. */
int kws_decide(kws_post_t *p, const int8_t *out, int smooth, int margin, int hold)
{
    for (uint32_t c = 0; c < NCLASS; c++) {
        if (p->fill == smooth)
            p->sum[c] = (int16_t)(p->sum[c] - p->hist[p->pos][c]);
        p->hist[p->pos][c] = out[c];
        p->sum[c] = (int16_t)(p->sum[c] + out[c]);
    }
    p->pos = (uint8_t)((p->pos + 1) % smooth);
    if (p->fill < smooth)
        p->fill++;
    if (p->hold) {
        p->hold--;
        return 0;
    }
    if (p->fill < smooth)
        return 0;
    int best = 0, second = -32768;
    for (int c = 1; c < NCLASS; c++)
        if (p->sum[c] > p->sum[best])
            best = c;
    for (int c = 0; c < NCLASS; c++)
        if (c != best && p->sum[c] > second)
            second = p->sum[c];
    if (best == 0 || p->sum[best] - second < margin * smooth)
        return 0;
    p->hold = (uint16_t)hold;
    return best;
}

gen_tables.py: generator of the MFCC tables in Q15 and of the decimation filter.

# -*- coding: utf-8 -*-
"""Таблицы для MFCC в Q15 и фильтра прореживания: python3 gen_tables.py > mfcc_tables.h
Чистый Python, без numpy: результат одинаков на любой машине."""
import math
import sys

sys.stdout.reconfigure(encoding="utf-8", newline="\n")

FS, FRAME, NFFT, NMEL, NCEP = 16000, 400, 512, 40, 10
FMIN, FMAX = 20.0, 8000.0
FIR_N, FIR_FS, FIR_FC, FIR_S, CIC_GAIN = 64, 64000.0, 8000.0, 23, 81


def mel(f):
    return 2595.0 * math.log10(1.0 + f / 700.0)


def q15(v):
    return min(32767, int(math.floor(v * 32768.0 + 0.5)))


def emit(ctype, name, vals, per=10):
    print("static const %s %s[%d] = {" % (ctype, name, len(vals)))
    for i in range(0, len(vals), per):
        print("    " + ", ".join("%d" % v for v in vals[i:i + per]) + ",")
    print("};\n")


print("/* Создано gen_tables.py — не править руками. */")
print("#ifndef MFCC_TABLES_H\n#define MFCC_TABLES_H\n#include <stdint.h>\n")

# половина окна Хэмминга: w[n] = 0,54 − 0,46·cos(2πn/(N−1)), вторая половина зеркальна
emit("int16_t", "hamming_q15", [q15(0.54 - 0.46 * math.cos(2 * math.pi * n / (FRAME - 1))) for n in range(FRAME // 2)])

# синус на 3/4 периода: sin(2πk/512), k = 0…383; косинус — со сдвигом на 128
emit("int16_t", "sin_q15", [q15(math.sin(2 * math.pi * k / NFFT)) for k in range(3 * NFFT // 4)])

# мел-фильтры: 42 точки поровну по шкале мел; бин k лежит между точками band и band + 1
pts = [mel(FMIN) + (mel(FMAX) - mel(FMIN)) * i / (NMEL + 1) for i in range(NMEL + 2)]
first, band, wgt = None, [], []
for k in range(NFFT // 2 + 1):
    m = mel(k * FS / NFFT)
    if m < pts[0] or m >= pts[-1]:
        continue
    if first is None:
        first = k
    b = max(i for i in range(NMEL + 1) if pts[i] <= m)
    band.append(b)
    wgt.append(q15((m - pts[b]) / (pts[b + 1] - pts[b])))
print("#define MEL_FIRST %d   /* первый бин внутри 20…8000 Гц */" % first)
print("#define MEL_BINS  %d\n" % len(band))
emit("uint8_t", "mel_band", band, 16)
emit("int16_t", "mel_w_q15", wgt)

# log2(1 + i/64) в Q15, 65 узлов
emit("uint16_t", "log2_q15", [int(math.floor(math.log2(1.0 + i / 64.0) * 32768.0 + 0.5)) for i in range(65)])

# ДКП-II: sqrt(2/M)·cos(π·i·(m + 0,5)/M), i = 0…9, m = 0…39
emit("int16_t", "dct_q15", [int(math.floor(math.sqrt(2.0 / NMEL) * math.cos(math.pi * i * (m + 0.5) / NMEL) * 32768.0 + 0.5))
                             for i in range(NCEP) for m in range(NMEL)])

# КИХ после CIC: оконный sinc (Хэмминг), срез 8 кГц при 64 кГц; усиление 2^23/81 убирает 48^4/2^16 = 81
h = []
for n in range(FIR_N):
    t = n - (FIR_N - 1) / 2.0
    x = 2.0 * FIR_FC / FIR_FS * t
    s = 1.0 if x == 0 else math.sin(math.pi * x) / (math.pi * x)
    h.append(2.0 * FIR_FC / FIR_FS * s * (0.54 - 0.46 * math.cos(2 * math.pi * n / (FIR_N - 1))))
g = sum(h)
emit("int16_t", "fir_q", [int(math.floor(v / g * (1 << FIR_S) / CIC_GAIN + 0.5)) for v in h])
print("#endif")

Next: section 121в "STM32N6 and a battery-powered camera with a neural network: YOLO on the Neural-ART NPU, low power, an outdoor device": a camera on a battery, the NPU runs YOLO, the core parses the boxes, the device sleeps in Standby; the energy budget and the specification.