Spectral Features on Hardware

Five descriptors from one FFT, a few passes over the magnitude bins

The expensive part of spectral feature extraction on a microcontroller is the FFT, and you have already paid for it: the spectrogram front end produces a magnitude spectrum every hop. The features themselves are cheap reductions over those bins, a few sums and one or two divides per frame. The whole family rides for almost free on top of a transform you were running anyway.

This page assumes the real-time block-FFT loop from the STFT embedded page: a ring buffer for overlap, a windowed frame, a vendor FFT, and a one-sided magnitude array mag[NBINS]. That FFT is the platform-specific part (CMSIS-DSP on the NUCLEO-F446RE, ESP-DSP on the ESP32-S3, both per ADR-005); it is covered on the STFT page. The feature code here is plain portable float C that runs unchanged on either chip, so the C below is shown once and a short ESP32-S3 note covers what each vendor library adds.

---

They all read the same array

Look back at the definitions on the main page and a pattern jumps out: every per-frame feature is a weighted sum over the same mag[] bins.

Feature	What to accumulate over the bins
Centroid	\(\sum_k f_k\,X[k]\) and \(\sum_k X[k]\)
Spread	also \(\sum_k f_k^2\,X[k]\) (then combine with the centroid)
Rolloff	the running cumulative of \(X[k]\) until it passes \(0.85\sum_k X[k]\)
Flatness	\(\sum_k \ln(X[k]^2)\) and \(\sum_k X[k]^2\)

So one pass over mag[] collects every accumulator the centroid, spread, and flatness need; a second short pass finds the rolloff once the total is known. There is no second FFT, no extra buffer, just arithmetic over bins already in memory.

A practical detail first: skip the DC bin (k = 0). A small DC or near-DC offset, which survives even after windowing, sits at \(f = 0\) and drags the centroid down without telling you anything about the sound. Start every loop at k = 1.

One pass for centroid, spread, and flatness

The natural form is a single float loop (no CMSIS-DSP call is needed: this is a multi-moment accumulation, not a transform). The bin frequency is \(f_k = k\,f_s/L\); accumulate in bin units and scale once at the end, so the loop carries no multiply by fs/L.

#include <math.h>                 // logf, sqrtf, expf

// mag[]   : NBINS one-sided magnitudes from the FFT (see the STFT page)
// returns : centroid and spread in Hz, flatness in [0,1]
void spectral_shape(const float *mag, int nbins, float fs, int nfft,
                    float *centroid_hz, float *spread_hz, float *flatness)
{
    float s0 = 0.0f;             // sum mag
    float s1 = 0.0f;             // sum k*mag        (centroid numerator, in bins)
    float s2 = 0.0f;             // sum k*k*mag      (for the spread)
    float sum_pow = 0.0f;        // sum mag^2        (flatness denominator)
    float sum_lnp = 0.0f;        // sum ln(mag^2)    (flatness numerator, in log)

    for (int k = 1; k < nbins; k++) {     // skip DC at k = 0
        float m = mag[k];
        s0 += m;
        s1 += k * m;
        s2 += (float)k * k * m;
        float p = m * m + 1e-20f;         // floor so ln() never sees zero
        sum_pow += p;
        sum_lnp += logf(p);
    }

    float bin_hz = fs / (float)nfft;      // width of one bin
    float mu_bins = (s0 > 0.0f) ? s1 / s0 : 0.0f;
    float var_bins = (s0 > 0.0f) ? s2 / s0 - mu_bins * mu_bins : 0.0f;
    if (var_bins < 0.0f) var_bins = 0.0f; // guard rounding at near-zero spread

    *centroid_hz = mu_bins * bin_hz;
    *spread_hz   = sqrtf(var_bins) * bin_hz;

    int n = nbins - 1;                     // bins actually summed (DC skipped)
    float gmean = expf(sum_lnp / n);       // geometric mean of power
    float amean = sum_pow / n;             // arithmetic mean of power
    *flatness = gmean / amean;             // the 1e-20f floor keeps amean > 0
}

Two numerical points carry over from the desktop code. The spread uses the identity \(\sigma^2 = \overline{f^2} - \mu^2\), which is one pass but can go slightly negative from rounding when the true spread is near zero (a pure tone); clamp it before the sqrtf. And the flatness logarithm needs a floor on the power, or a silent bin sends logf to -inf and the geometric mean collapses to zero. That same floor means a fully silent frame returns a flatness of 1 (every floored bin is equal, so geometric and arithmetic means coincide), matching the Python module rather than returning a special-cased zero.

A second short pass for rolloff

Rolloff needs the grand total before it can find the 85% crossing, so it is a separate walk up the bins. It is the cheapest of the lot: an add and a compare per bin, no multiply.

float spectral_rolloff(const float *mag, int nbins, float fs, int nfft,
                       float roll_percent)
{
    float total = 0.0f;
    for (int k = 1; k < nbins; k++) total += mag[k];
    if (total <= 0.0f) return 0.0f;             // silent frame: rolloff undefined

    float threshold = roll_percent * total;     // e.g. 0.85f
    float cumulative = 0.0f;
    for (int k = 1; k < nbins; k++) {
        cumulative += mag[k];
        if (cumulative >= threshold)
            return (float)k * fs / (float)nfft; // first bin that reaches it
    }
    return fs * 0.5f;                            // unreachable for roll_percent <= 1
                                                 // and total > 0; defensive fallback
}

If the loop has already computed s0 (the magnitude total) in spectral_shape, pass it in and drop the first pass here; the two functions are split only for clarity.

Flux needs one frame of memory

Flux is the only feature that looks across frames, so it needs the previous magnitude spectrum kept in a static buffer. That is one extra NBINS array, and the cost of normalising each frame to unit energy if you want flux to ignore loudness.

#define NBINS (NFFT / 2 + 1)    // NFFT comes from the STFT page; same size as mag[]

// keeps the previous (normalised) magnitude spectrum between calls
float spectral_flux(const float *mag, int nbins)
{
    // prev[] is sized to the full one-sided spectrum; nbins must not exceed NBINS
    static float prev[NBINS] = {0};

    float norm = 0.0f;
    for (int k = 1; k < nbins; k++) norm += mag[k] * mag[k];
    norm = sqrtf(norm) + 1e-20f;

    float flux = 0.0f;
    for (int k = 1; k < nbins; k++) {
        float cur = mag[k] / norm;              // unit-energy normalised
        float d = cur - prev[k];
        flux += d * d;
        prev[k] = cur;
    }
    return sqrtf(flux);
}

Because prev[] is a file-scope static, this function is not reentrant; like the FFT loop it feeds, only one context may call it. The first call compares against an all-zero prev[], so discard the very first flux value (it is a cold-start artefact, exactly like the first spectrogram column).

The compute budget

Every feature is \(O(K)\) in the number of bins \(K = L/2 + 1\), and the per-bin work is tiny: a few single-cycle MACs for the moments, an add and compare for the rolloff, and (the only real cost) one logf per bin for the flatness. Set against the FFT that produced the bins and the hop period that bounds everything, the features are in the noise.

NFFT = 1024
K = NFFT // 2 + 1                       # one-sided bins
fs = 16000.0
HOP = 512
f_cpu = 180e6                           # F446RE core clock

# rough per-bin cost (cycles). MACs are ~1 cycle on the M4F; a software logf is
# the expensive operation: the M4F FPU has no transcendental unit, so logf is a
# newlib routine, typically 50-150 cycles. These are order-of-magnitude, not a
# benchmark: confirm on real hardware with DWT->CYCCNT.
mac_bins = K * 6                        # centroid+spread+flux: a handful of MACs/bin
rolloff_bins = K * 2                    # one add + compare/bin, two passes
logf_bins = K * 80                      # flatness: ~1 logf/bin, ~80 cycles (50-150)
feat_cycles = mac_bins + rolloff_bins + logf_bins

# the FFT that produced the bins: ~N/2 log2 N butterflies (order, see STFT page)
fft_butterflies = (NFFT // 2) * int(np.log2(NFFT))

hop_cycles = (HOP / fs) * f_cpu         # cycles available before the next frame
print(f"features ~{feat_cycles:,} cycles/frame  ({K} bins)")
print(f"FFT      ~{fft_butterflies:,} butterflies/frame")   # butterflies != cycles

print(f"budget    {hop_cycles:,.0f} cycles/hop ({1000*HOP/fs:.0f} ms at 180 MHz)")

# the whole feature set is dominated by the flatness log, and even so is a small
# slice of the per-hop budget
assert feat_cycles < 0.05 * hop_cycles  # under 5% of the hop budget
print(f"-> spectral features ~{100*feat_cycles/hop_cycles:.1f}% of the hop budget")

features ~45,144 cycles/frame  (513 bins)
FFT      ~5,120 butterflies/frame
budget    5,760,000 cycles/hop (32 ms at 180 MHz)
-> spectral features ~0.8% of the hop budget

The lesson the numbers tell: the whole feature set costs around one percent or less of the per-hop budget, so once the FFT fits, the features fit easily on top of it. (The butterfly count is not directly comparable to the feature cycle count, since each butterfly is several cycles; the honest comparison is each against the hop budget, which both clear with room to spare.) The flatness logarithm dominates the feature cost, so on a part without a fast logf it is the first thing to drop or approximate.

On the ESP32-S3

The C above is platform-agnostic: the accumulation loops are ordinary single-precision float arithmetic, and logf, sqrtf, and expf are the same newlib routines on the ESP32-S3’s Xtensa core as on the M4F. So the four functions compile and run unchanged. What differs between the two chips is the FFT that fills mag[] (CMSIS-DSP’s arm_rfft_fast_f32 versus ESP-DSP’s dsps_fft2r_fc32, both on the STFT page), not the feature math.

Where ESP-DSP earns its keep is the accumulations, which are dot products in disguise and map onto its SIMD-accelerated dsps_dotprod_f32:

• \(s_0 = \sum_k \texttt{mag}[k]\) is mag dotted with a vector of ones,
• \(s_1 = \sum_k k\,\texttt{mag}[k]\) is mag dotted with a precomputed bin-index ramp [0, 1, 2, ...],
• the flatness denominator \(\sum_k \texttt{mag}[k]^2\) is mag dotted with itself.

#include "esp_dsp.h"

static float ramp[NBINS];   // fill once at startup: ramp[k] = (float)k
static float ones[NBINS];   // fill once at startup: ones[k] = 1.0f

// Start at bin 1, length NBINS-1: the scalar loop skips DC (k=0), so the SIMD
// dot products must skip mag[0] too, or a surviving DC offset biases all three.
float s0, s1, sum_pow;
dsps_dotprod_f32(mag + 1, ones + 1, &s0, NBINS - 1);   // sum mag, DC excluded
dsps_dotprod_f32(mag + 1, ramp + 1, &s1, NBINS - 1);   // sum k*mag (centroid numerator)
dsps_dotprod_f32(mag + 1, mag + 1,  &sum_pow, NBINS - 1);

The flatness logarithm has no SIMD shortcut on either chip; it stays a per-bin logf, and remains the first feature to drop on a tight budget. The cumulative-sum rolloff is inherently sequential and gains nothing from a dot product, so it stays the scalar loop shown above. As always, confirm the real cost on hardware with a cycle counter rather than trusting the estimate.

When there is no FPU

The moments (centroid, spread, rolloff) are integer-friendly and survive the move to a no-FPU core untouched in spirit. Keep the magnitudes as int16_t or int32_t, accumulate the weighted sums in int64_t to avoid overflow (a 10-bit bin index times a 16-bit magnitude, summed over hundreds of bins, runs past 32 bits), and do the single final divide in fixed point. Working in bin index rather than hertz keeps the weights small integers and defers the one fs/L scaling to the end.

// integer centroid: returns the centroid in bin units, x256 fixed point.
// mag[] must be non-negative one-sided FFT magnitudes, not signed PCM samples.
uint32_t centroid_bins_q8(const int16_t *mag, int nbins)
{
    int64_t s0 = 0, s1 = 0;
    for (int k = 1; k < nbins; k++) {
        s0 += mag[k];
        s1 += (int64_t)k * mag[k];
    }
    return (s0 > 0) ? (uint32_t)((s1 << 8) / s0) : 0;  // bins, Q8
}

Flatness is the awkward one: the geometric mean is a logarithm, and an exact logf is costly without an FPU. The usual fixed-point substitute is an integer \(\log_2\) (a leading-zero count for the integer part, a small lookup table for the fraction), which carries the finite word-length effects of any fixed-point pipeline. If flatness is not pulling its weight as a feature, the cheapest decision is to drop it and keep the three moment features, which is often enough for a coarse voiced-versus-noise gate.

References