PSD Estimation on Hardware

Averaging periodograms on a microcontroller

A Welch PSD estimate on a microcontroller is three things: an FFT, a squared-magnitude accumulation buffer, and enough memory to hold the average. This page shows a fixed-point implementation for the two default platforms (ESP32-S3 and NUCLEO-F446RE, per ADR-005), covering the memory budget, the accumulation-in-place trick, and how to verify your result against the known quantization noise floor.

The computation in one sentence

Collect $K$ overlapping windowed segments, compute the real FFT of each, square and sum the magnitudes into an accumulation buffer, divide by $K$ at the end, and normalise for the window and sample rate. The result is a Welch PSD estimate in fixed-point units.

Memory budget

The dominant cost is the accumulation buffer. For a segment length $L$:

Buffer	Size	Example ($L = 256$)
Input segment (time domain)	$L$ × `float` or `int16`	1 KB (float32) or 512 B (int16)
FFT output (frequency)	$L$ × `float`	1 KB
Accumulation (PSD sum)	$(L/2+1)$ × `float` or `int32`	516 B (float) or 516 B (int32)
Window coefficients	$L$ × `float`	1 KB (can be stored in flash)

Total: roughly $2.5L$ words in RAM ($L + L + L/2$, with the window table in flash). For $L = 256$ with float32, about 2.5 KB, comfortable on both platforms. For $L = 1024$, about 10 KB: fine on the ESP32-S3 (512 KB SRAM), tighter on the F446RE (128 KB SRAM) but workable if the rest of the application is lean.

Fixed-point accumulation

If your FFT is fixed-point (CMSIS-DSP arm_rfft_q15 on the F446RE), square the real and imaginary parts, sum in int32_t, and convert to float only at the end when you divide by $K$ and normalise. This avoids accumulating floating-point round-off over many segments. (The ESP32-S3 path below uses ESP-DSP’s float FFT instead; with a float accumulator the round-off concern is smaller but the same structure applies.)

// CMSIS-DSP fixed-point Welch accumulation (Q15 FFT -> int32 accumulator).
// Assumes: L = segment length, K = number of segments.
// Input: q15_input[L] = windowed time-domain segment in Q1.15.
// Output: psd_accum[L/2+1] = summed squared magnitudes (int32).

#include "arm_math.h"

#define L        256
#define L_BINS   (L/2 + 1)

static q15_t         windowed[L];       // windowed segment, Q1.15
static q15_t         fft_out[L * 2];    // real/imag interleaved, Q?.15
static int32_t       psd_accum[L_BINS]; // squared magnitude sum
static arm_rfft_instance_q15 rfft_inst;

void welch_init(void) {
    arm_rfft_init_q15(&rfft_inst, L, 0, 1);  // forward RFFT, bit-reversed
    memset(psd_accum, 0, sizeof(psd_accum));
}

void welch_accumulate_segment(const q15_t *segment, const q15_t *window) {
    // Apply window
    arm_mult_q15(segment, window, windowed, L);

    // Real FFT
    arm_rfft_q15(&rfft_inst, windowed, fft_out);

    // Square and accumulate.  CMSIS-DSP arm_rfft_q15 packed format:
    //   fft_out[0] = Re(DC),  fft_out[1] = Re(Nyquist)
    //   for k=1..L/2-1:  fft_out[2k] = Re(k), fft_out[2k+1] = Im(k)
    int32_t acc;

    // DC bin (i=0): real part only.
    acc = (int32_t)fft_out[0] * fft_out[0];
    psd_accum[0] += acc;

    // Interior bins (i=1 .. L/2-1): real and imaginary parts.
    for (int i = 1; i < L/2; i++) {
        q15_t re = fft_out[2 * i];
        q15_t im = fft_out[2 * i + 1];
        acc = ((int32_t)re * re + (int32_t)im * im);
        psd_accum[i] += acc;
    }

    // Nyquist bin (i = L/2): real part at fft_out[1].
    acc = (int32_t)fft_out[1] * fft_out[1];
    psd_accum[L/2] += acc;
}

void welch_finalise(float *psd_float, float fs, int K, float window_power) {
    // Convert accumulated Q30 sums to float PSD [units²/Hz].
    // PSD = (psd_accum / K) / (fs * L * window_power)
    //    with interior bins ×2 for one-sided spectrum.
    float scale = 1.0f / (K * fs * L * window_power);
    for (int i = 0; i < L_BINS; i++) {
        psd_float[i] = (float)psd_accum[i] * scale;
    }
    // Double interior bins (one-sided)
    for (int i = 1; i < L/2; i++) {
        psd_float[i] *= 2.0f;
    }
}

Scaling warning: Q15 FFT output scaling depends on the CMSIS-DSP version. Recent versions divide by 2 at each radix-2 stage, so the output is scaled by $1/L$. Check your library’s documentation: the normalisation must match.

ESP32-S3: float Welch with ESP-DSP

The ESP32-S3 has a hardware FPU, so the natural path is ESP-DSP’s float radix-2 FFT. ESP-DSP has no dedicated real FFT: pack the real segment into the real slots of a complex buffer, run the complex FFT, bit-reverse, and split with dsps_cplx2reC_fc32. The twiddle table must be initialised once at startup with dsps_fft2r_init_fc32 before any FFT call.

// ESP-DSP float Welch accumulation for the ESP32-S3.
#include "esp_dsp.h"

#define L        256
#define L_BINS   (L/2 + 1)

static float y_cf[L * 2];        // complex work buffer: re/im interleaved
static float psd_accum_f[L_BINS];

void welch_init_esp32(void) {
    // Once, at startup: generate the twiddle table (check the return code).
    ESP_ERROR_CHECK(dsps_fft2r_init_fc32(NULL, L));
    memset(psd_accum_f, 0, sizeof(psd_accum_f));
}

void welch_accumulate_segment_esp32(const float *segment, const float *window) {
    for (int i = 0; i < L; i++) {
        y_cf[2 * i]     = segment[i] * window[i];
        y_cf[2 * i + 1] = 0.0f;
    }
    dsps_fft2r_fc32(y_cf, L);
    dsps_bit_rev_fc32(y_cf, L);
    dsps_cplx2reC_fc32(y_cf, L);   // unpack to L/2+1 one-sided complex bins

    for (int i = 0; i < L_BINS; i++) {
        float re = y_cf[2 * i];
        float im = y_cf[2 * i + 1];
        psd_accum_f[i] += re * re + im * im;
    }
}

The finalisation step is the same as welch_finalise() above (divide by $K$, normalise by $f_s L$ and the window power, double the interior bins), minus the Q15-to-float conversion.

Verifying against the quantization floor

The easiest validation on an MCU: short the ADC input to a quiet reference (internal bandgap or a well-decoupled mid-rail divider), collect enough samples for $K \approx 30$ segments, compute the Welch PSD, and check that:

The PSD is flat (white), no peaks above the confidence interval.
The integrated power matches the measured time-domain variance: $\int_0^{f_s/2} \hat{S}(f)\ df \approx \sigma_x^2$.
If the ADC is quantization-limited, $\sigma_x^2 \approx Q^2/12$ within a factor of 2. A larger variance means thermal or interference dominance.

// Verify: integrate the PSD and compare with time-domain variance.
float psd_integral = 0.0f;
float df = fs / L;  // bin width
for (int i = 1; i < L/2; i++) {
    psd_integral += psd_float[i] * df;
}
// psd_integral should approximately equal time_domain_variance.
// If not, check the window normalisation and the one-sided/doubling factor.

Practical limits

Segment length vs memory. On an ESP32-S3 with ESP-DSP (dsps_fft2r_fc32, float FFT), $L = 512$ is a good default: gives $\Delta f \approx 20$ Hz at 10 kS/s, uses about 5 KB RAM, and allows $K \approx 38$ segments (50% overlap) from a 1-second noise record at 10 kS/s.

Fixed-point noise floor. With Q15 arithmetic, the squared magnitude accumulator (int32_t) overflows if the sum exceeds $2^{31}-1 \approx 2.1 \times 10^9$. The worst case is brutal: a full-scale bin ($\pm 32767$) contributes $\mathrm{re}^2 + \mathrm{im}^2 \approx 2.1 \times 10^9$ in a single segment, so int32_t offers no worst-case headroom at all. In practice arm_rfft_q15’s internal $1/L$ scaling keeps bins far below full scale (a full-scale input tone lands near $\pm 16384$, whose square is $\approx 2.7 \times 10^8$, overflowing around $K = 8$; broadband noise spreads power over $L/2$ bins and allows $K$ in the hundreds), but you should not budget on “in practice”. Either: (a) right-shift each squared magnitude by a few bits before accumulating (reducing precision), or (b) use int64_t for the accumulator (cheap on Cortex-M4: SMLAL), which removes the ceiling entirely.

Real-time Welch. A running Welch estimate updates the PSD as new segments arrive by subtracting the oldest segment’s contribution and adding the newest, a circular accumulation buffer. The memory cost doubles (you must store each segment’s squared magnitudes separately), but the PSD updates continuously instead of requiring a full re-collection. This is the approach for a real-time spectrum display or an adaptive noise floor tracker.

When not to bother

If all you need is “is the noise white?”, compute the lag-1 ACF on the MCU (see the random processes embedded page). A flat ACF at nonzero lags answers the question in a few multiply-accumulates per sample, with no FFT and no accumulation buffer. Reserve the Welch estimate for when you need the spectral shape: the slope of $1/f^\alpha$ noise, the width of a resonant peak, or the flatness of a noise floor you intend to whiten.

The companion topic on noise generation closes the loop: generate noise with a specified PSD, then verify with Welch that you got what you asked for.

References

--- title: "PSD Estimation on Hardware" subtitle: "Averaging periodograms on a microcontroller" bibliography: ../../references.bib --- A Welch PSD estimate on a microcontroller is three things: an FFT, a squared-magnitude accumulation buffer, and enough memory to hold the average. This page shows a fixed-point implementation for the two default platforms (ESP32-S3 and NUCLEO-F446RE, per [ADR-005](../../docs/adr/005-embedded-platforms.md)), covering the memory budget, the accumulation-in-place trick, and how to verify your result against the known quantization noise floor. <hr> ## The computation in one sentence Collect $K$ overlapping windowed segments, compute the real FFT of each, square and sum the magnitudes into an accumulation buffer, divide by $K$ at the end, and normalise for the window and sample rate. The result is a Welch PSD estimate in fixed-point units. <hr> ## Memory budget The dominant cost is the accumulation buffer. For a segment length $L$: | Buffer | Size | Example ($L = 256$) | |---|---|---| | Input segment (time domain) | $L$ × `float` or `int16` | 1 KB (float32) or 512 B (int16) | | FFT output (frequency) | $L$ × `float` | 1 KB | | Accumulation (PSD sum) | $(L/2+1)$ × `float` or `int32` | 516 B (float) or 516 B (int32) | | Window coefficients | $L$ × `float` | 1 KB (can be stored in flash) | Total: roughly $2.5L$ words in RAM ($L + L + L/2$, with the window table in flash). For $L = 256$ with `float32`, about 2.5 KB, comfortable on both platforms. For $L = 1024$, about 10 KB: fine on the ESP32-S3 (512 KB SRAM), tighter on the F446RE (128 KB SRAM) but workable if the rest of the application is lean. <hr> ## Fixed-point accumulation If your FFT is fixed-point (CMSIS-DSP `arm_rfft_q15` on the F446RE), square the real and imaginary parts, sum in `int32_t`, and convert to float only at the end when you divide by $K$ and normalise. This avoids accumulating floating-point round-off over many segments. (The ESP32-S3 path below uses ESP-DSP's float FFT instead; with a float accumulator the round-off concern is smaller but the same structure applies.) ```c // CMSIS-DSP fixed-point Welch accumulation (Q15 FFT -> int32 accumulator). // Assumes: L = segment length, K = number of segments. // Input: q15_input[L] = windowed time-domain segment in Q1.15. // Output: psd_accum[L/2+1] = summed squared magnitudes (int32). #include "arm_math.h" #define L 256 #define L_BINS (L/2 + 1) static q15_t windowed[L]; // windowed segment, Q1.15 static q15_t fft_out[L * 2]; // real/imag interleaved, Q?.15 static int32_t psd_accum[L_BINS]; // squared magnitude sum static arm_rfft_instance_q15 rfft_inst; void welch_init(void) { arm_rfft_init_q15(&rfft_inst, L, 0, 1); // forward RFFT, bit-reversed memset(psd_accum, 0, sizeof(psd_accum)); } void welch_accumulate_segment(const q15_t *segment, const q15_t *window) { // Apply window arm_mult_q15(segment, window, windowed, L); // Real FFT arm_rfft_q15(&rfft_inst, windowed, fft_out); // Square and accumulate. CMSIS-DSP arm_rfft_q15 packed format: // fft_out[0] = Re(DC), fft_out[1] = Re(Nyquist) // for k=1..L/2-1: fft_out[2k] = Re(k), fft_out[2k+1] = Im(k) int32_t acc; // DC bin (i=0): real part only. acc = (int32_t)fft_out[0] * fft_out[0]; psd_accum[0] += acc; // Interior bins (i=1 .. L/2-1): real and imaginary parts. for (int i = 1; i < L/2; i++) { q15_t re = fft_out[2 * i]; q15_t im = fft_out[2 * i + 1]; acc = ((int32_t)re * re + (int32_t)im * im); psd_accum[i] += acc; } // Nyquist bin (i = L/2): real part at fft_out[1]. acc = (int32_t)fft_out[1] * fft_out[1]; psd_accum[L/2] += acc; } void welch_finalise(float *psd_float, float fs, int K, float window_power) { // Convert accumulated Q30 sums to float PSD [units²/Hz]. // PSD = (psd_accum / K) / (fs * L * window_power) // with interior bins ×2 for one-sided spectrum. float scale = 1.0f / (K * fs * L * window_power); for (int i = 0; i < L_BINS; i++) { psd_float[i] = (float)psd_accum[i] * scale; } // Double interior bins (one-sided) for (int i = 1; i < L/2; i++) { psd_float[i] *= 2.0f; } } ``` **Scaling warning:** Q15 FFT output scaling depends on the CMSIS-DSP version. Recent versions divide by 2 at each radix-2 stage, so the output is scaled by $1/L$. Check your library's documentation: the normalisation must match. <hr> ## ESP32-S3: float Welch with ESP-DSP The ESP32-S3 has a hardware FPU, so the natural path is ESP-DSP's float radix-2 FFT. ESP-DSP has no dedicated real FFT: pack the real segment into the real slots of a complex buffer, run the complex FFT, bit-reverse, and split with `dsps_cplx2reC_fc32`. The twiddle table must be initialised **once** at startup with `dsps_fft2r_init_fc32` before any FFT call. ```c // ESP-DSP float Welch accumulation for the ESP32-S3. #include "esp_dsp.h" #define L 256 #define L_BINS (L/2 + 1) static float y_cf[L * 2]; // complex work buffer: re/im interleaved static float psd_accum_f[L_BINS]; void welch_init_esp32(void) { // Once, at startup: generate the twiddle table (check the return code). ESP_ERROR_CHECK(dsps_fft2r_init_fc32(NULL, L)); memset(psd_accum_f, 0, sizeof(psd_accum_f)); } void welch_accumulate_segment_esp32(const float *segment, const float *window) { for (int i = 0; i < L; i++) { y_cf[2 * i] = segment[i] * window[i]; y_cf[2 * i + 1] = 0.0f; } dsps_fft2r_fc32(y_cf, L); dsps_bit_rev_fc32(y_cf, L); dsps_cplx2reC_fc32(y_cf, L); // unpack to L/2+1 one-sided complex bins for (int i = 0; i < L_BINS; i++) { float re = y_cf[2 * i]; float im = y_cf[2 * i + 1]; psd_accum_f[i] += re * re + im * im; } } ``` The finalisation step is the same as `welch_finalise()` above (divide by $K$, normalise by $f_s L$ and the window power, double the interior bins), minus the Q15-to-float conversion. <hr> ## Verifying against the quantization floor The easiest validation on an MCU: short the ADC input to a quiet reference (internal bandgap or a well-decoupled mid-rail divider), collect enough samples for $K \approx 30$ segments, compute the Welch PSD, and check that: 1. The PSD is flat (white), no peaks above the confidence interval. 2. The integrated power matches the measured time-domain variance: $\int_0^{f_s/2} \hat{S}(f)\ df \approx \sigma_x^2$. 3. If the ADC is quantization-limited, $\sigma_x^2 \approx Q^2/12$ within a factor of 2. A larger variance means thermal or interference dominance. ```c // Verify: integrate the PSD and compare with time-domain variance. float psd_integral = 0.0f; float df = fs / L; // bin width for (int i = 1; i < L/2; i++) { psd_integral += psd_float[i] * df; } // psd_integral should approximately equal time_domain_variance. // If not, check the window normalisation and the one-sided/doubling factor. ``` <hr> ## Practical limits **Segment length vs memory.** On an ESP32-S3 with ESP-DSP (`dsps_fft2r_fc32`, float FFT), $L = 512$ is a good default: gives $\Delta f \approx 20$ Hz at 10 kS/s, uses about 5 KB RAM, and allows $K \approx 38$ segments (50% overlap) from a 1-second noise record at 10 kS/s. **Fixed-point noise floor.** With Q15 arithmetic, the squared magnitude accumulator (`int32_t`) overflows if the sum exceeds $2^{31}-1 \approx 2.1 \times 10^9$. The worst case is brutal: a full-scale bin ($\pm 32767$) contributes $\mathrm{re}^2 + \mathrm{im}^2 \approx 2.1 \times 10^9$ in a *single* segment, so `int32_t` offers no worst-case headroom at all. In practice `arm_rfft_q15`'s internal $1/L$ scaling keeps bins far below full scale (a full-scale input tone lands near $\pm 16384$, whose square is $\approx 2.7 \times 10^8$, overflowing around $K = 8$; broadband noise spreads power over $L/2$ bins and allows $K$ in the hundreds), but you should not budget on "in practice". Either: (a) right-shift each squared magnitude by a few bits before accumulating (reducing precision), or (b) use `int64_t` for the accumulator (cheap on Cortex-M4: `SMLAL`), which removes the ceiling entirely. **Real-time Welch.** A running Welch estimate updates the PSD as new segments arrive by subtracting the oldest segment's contribution and adding the newest, a circular accumulation buffer. The memory cost doubles (you must store each segment's squared magnitudes separately), but the PSD updates continuously instead of requiring a full re-collection. This is the approach for a real-time spectrum display or an adaptive noise floor tracker. <hr> ## When not to bother If all you need is "is the noise white?", compute the lag-1 ACF on the MCU (see the [random processes embedded page](../random-processes/embedded.qmd)). A flat ACF at nonzero lags answers the question in a few multiply-accumulates per sample, with no FFT and no accumulation buffer. Reserve the Welch estimate for when you need the spectral *shape*: the slope of $1/f^\alpha$ noise, the width of a resonant peak, or the flatness of a noise floor you intend to whiten. The companion topic on **[noise generation](../noise-generation/)** closes the loop: generate noise with a specified PSD, then verify with Welch that you got what you asked for. ## References ::: {#refs} :::