File psd_core.h¶

FileList > doppler > psd > psd_core.h

Go to the source code of this file

PSD — averaging power-spectral-density estimator (Welch's method), its per-frame kernel, and spectral measurement suite. More...

• #include "doppler/clib_common.h"
• #include "doppler/dp_state.h"
• #include "doppler/jm_perf.h"
• #include "doppler/fft/fft_core.h"
• #include "doppler/acc_trace/acc_trace_core.h"
• #include "doppler/f32_buffer/f32_buffer_core.h"

Classes¶

Public Types¶

Public Functions¶

Macros¶

Detailed Description¶

Two objects live here. dp_psd_kernel_t turns ONE frame into nfft DC-centred power bins and reads them against full scale, with no averaging: the window, the transform and the dBFS reference. PSD holds one by value and folds each frame it makes into a running average; the Spectrogram holds one and emits each frame as a row. A caller that wants single-frame spectra holds a kernel, not a PSD.

A stateful, C-first periodogram averager (Welch's method). It composes the existing pieces of the library rather than re-implementing them: an dp_fft_state_t forward plan, a spectral window (Hann, Kaiser, Blackman-Harris or rectangular) with its coherent gain and ENBW, an dp_acc_trace_state_t per-bin power averager (mean / EMA / max-hold / min-hold), and the spectral free functions (dp_magnitude_db_cf32, dp_find_peaks_f32, dp_obw_from_power, dp_noise_floor_db) for the derived measurements.

Feed complex baseband frames with dp_psd_accumulate(); each length-n frame is windowed, FFT'd, converted to power, fftshifted to DC-centred order and folded into the running average. Then read: * dp_psd_psd_db() : averaged power spectrum, dB (peak reads tone power) * dp_psd_psd_dbhz() : averaged PSD, dB/Hz (ENBW / fs normalised) * dp_psd_band_power() / dp_psd_total_band_power() : integrated band power, dB * dp_psd_occupied_bw() : occupied bandwidth, Hz * dp_psd_noise_floor() / dp_psd_snr() / dp_psd_sfdr() : level statistics, dB

All spectra are DC-centred (fftshift), matching find_peaks_f32's bin -> frequency convention (bin i maps to (i - nfft/2)/nfft in normalised frequency, so spectral peaks are obtained idiomatically with find_peaks_f32(w.psd_db(), n_peaks, min_db)).

Lifecycle: create -> (accumulate / reset)* -> (measurement getters)* -> destroy

The kernel's: dp_psd_kernel_init() -> (dp_psd_frame_power / _linear / _db)* -> dp_psd_kernel_free().

Public Types Documentation¶

enum psd__core_8h_1adc29c2ff13d900c2f185ee95427fb06c¶

The window of dp_psd_create() anddp_psd_kernel_init() , and of the objects built on them: all four for the Spectrogram (which holds a kernel), a subset for carrier acquisition (dp_carrier_acq_create() says which).

enum psd__core_8h_1adc29c2ff13d900c2f185ee95427fb06c {
    DP_PSD_WINDOW_HANN = 0,
    DP_PSD_WINDOW_KAISER = 1,
    DP_PSD_WINDOW_BLACKMAN_HARRIS = 2,
    DP_PSD_WINDOW_RECT = 3
};

HANN is 0, so a zeroed configuration is Hann. The Python faces name the same windows "hann", "kaiser", "blackman-harris" and "rect", bound to these constants by name (just-makeit.toml, [[enum]] psd_window); carrier acquisition's face names its subset (carrier_acq_window).

Public Functions Documentation¶

function dp_psd_accumulate¶

Window, FFT and fold complex baseband frames into the average. Folds every frame the stream completes: frame k covers stream samples k*hop through k*hop + n - 1 (hop = n unless dp_psd_create_hop() said otherwise). A trailing partial frame is carried to the next call, not dropped: the framer holds it. A single call of whole frames at hop = n averages what dp_psd_accumulate always has.

void dp_psd_accumulate (
    dp_psd_state_t * state,
    const float _Complex * x,
    size_t x_len
) 

Parameters:

• state Must be non-NULL.
• x Complex baseband samples (cf32).
• x_len Number of samples in x.

>>> import numpy as np
>>> from doppler.spectral import PSD
>>> n = 64
>>> w = PSD(n=n, fs=1.0, window="hann", mode="mean")
>>> k = 8
>>> x = np.exp(2j*np.pi*k*np.arange(n)/n).astype(np.complex64)
>>> for _ in range(4):
...     w.accumulate(x)
>>> psd = w.psd_db()
>>> psd.shape
(64,)
>>> int(np.argmax(psd)) == n // 2 + k
True
>>> w.count
4

function dp_psd_accumulate_real¶

Window, zero-pad, FFT and fold real frames into the average. The real-input counterpart to dp_psd_accumulate() : each length-n frame is windowed, zero-padded to nfft, transformed and folded as a DC-centred two-sided power spectrum (a real frame is Hermitian, so the +k and -k bins carry equal power). Read the one-sided fold withdp_psd_power_onesided() .

void dp_psd_accumulate_real (
    dp_psd_state_t * state,
    const float * x,
    size_t x_len
) 

This call stays on whole frames at hop = n: it folds floor(n_in / n) frames and drops a trailing partial frame. It neither reads nor writes the carry of the complex accumulate call, so mixing the two on one object mixes two streams, and a hop below n is not used by this call.

Parameters:

• state Must be non-NULL.
• x Real samples (f32).
• x_len Number of samples in x.

function dp_psd_band_power¶

Integrated power per band in dB. bands is a flat array of[lo0, hi0, lo1, hi1, ...] band edges in Hz; the output holds one dB value per band (n_bands = bands_len / 2). Edges are clamped to the analysed span; a band fully outside the span integrates to the dB floor. Returns 0 (Python None) when there is nothing to report: before any frame is accumulated, or whenbands holds no complete lo/hi pair.

size_t dp_psd_band_power (
    dp_psd_state_t * state,
    const double * bands,
    size_t bands_len,
    float * out,
    size_t max_out
) 

Parameters:

• state Must be non-NULL.
• bands Flat [lo,hi,...] band edges, Hz.
• bands_len Number of edge values (2 * n_bands).
• out Destination, at least n_bands float32 elements.
• max_out Capacity of out in elements. Emission stops there, so the return value is the number actually written.

Returns:

The number of bands written, min(n_bands, max_out); 0 (Python None) when there is nothing to report.

>>> import numpy as np
>>> from doppler.spectral import PSD
>>> w = PSD(n=64, fs=1.0, window="hann", mode="mean")
>>> w.accumulate(np.ones(64, dtype=np.complex64))
>>> pb = w.band_power(np.array([-0.5, 0.0, 0.0, 0.5]))
>>> pb.shape
(2,)

function dp_psd_band_power_max_out¶

Output capacity hint for band_power(); 0 (binding sizes from bands).

size_t dp_psd_band_power_max_out (
    dp_psd_state_t * state
) 

function dp_psd_create¶

Create an averaging PSD estimator.

dp_psd_state_t * dp_psd_create (
    size_t n,
    double fs,
    int window,
    float beta,
    size_t pad,
    double full_scale,
    size_t bits,
    int mode,
    double alpha
) 

Parameters:

• n Window / frame length in samples. Must be >= 2.
• fs Sample rate in Hz (used for dB/Hz and band frequencies); finite and > 0.
• window The window: "hann" (the default), "kaiser", "blackman-harris" or "rect" (no taper: best resolution, worst leakage).
• beta Kaiser beta (ignored by every other window). One that makes the window non-finite NaN, or large enough to overflow it is refused.
• pad Zero-pad factor (>= 1; 0 is refused); nfft = next_pow_two(n * pad).
• full_scale Amplitude that reads 0 dBFS in the dB getters; finite and > 0. Ignored when bits > 0.
• bits ADC depth: when > 0, sets full_scale = 2^(bits-1) (the single definition of the dBFS reference); 0 = use full_scale directly. At most 64.
• mode Averaging mode: "mean", "exp", "maxhold" or "minhold".
• alpha EMA smoothing factor, exp mode only, where it must lie in (0, 1]: the rule is the AccTrace averager's, whose create refuses anything else, so this does too. Ignored, and so not checked, in the other modes.

Returns:

Heap-allocated state, or NULL on invalid argument or OOM. A window whose coherent gain sum(w) is 0 Hann at n = 2, whose symmetric form is all zeros is refused too: every reading divides by it. So is an n or n * pad too large for its buffers to be sized, rather than wrapped to a small one.

Note:

Caller must call dp_psd_destroy() when done.

>>> from doppler.spectral import PSD
>>> w = PSD(n=1024, fs=1.0e6, window="kaiser", beta=8.0, mode="mean")
>>> w.n, w.fs
(1024, 1000000.0)
>>> round(w.rbw / (w.fs / w.n), 3) == round(w.enbw, 3)
True

function dp_psd_create_hop¶

Create an averaging PSD estimator whose frames are spaced by a hop.

dp_psd_state_t * dp_psd_create_hop (
    size_t n,
    double fs,
    int window,
    float beta,
    size_t pad,
    double full_scale,
    size_t bits,
    int mode,
    double alpha,
    size_t hop
) 

Frame k covers stream samples k*hop through k*hop + n - 1. A hop equal to n gives non-overlapping frames. A trailing partial frame is carried to the next accumulate call, not dropped.

Parameters:

• n Window / frame length in samples. Must be >= 2.
• fs Sample rate in Hz (used for dB/Hz and band frequencies); finite and > 0.
• window The window: "hann" (the default), "kaiser", "blackman-harris" or "rect" (no taper: best resolution, worst leakage).
• beta Kaiser beta (ignored by every other window). One that makes the window non-finite NaN, or large enough to overflow it is refused.
• pad Zero-pad factor (>= 1; 0 is refused); nfft = next_pow_two(n * pad).
• full_scale Amplitude that reads 0 dBFS in the dB getters; finite and > 0. Ignored when bits > 0.
• bits ADC depth: when > 0, sets full_scale = 2^(bits-1) (the single definition of the dBFS reference); 0 = use full_scale directly. At most 64.
• mode Averaging mode: "mean", "exp", "maxhold" or "minhold".
• alpha EMA smoothing factor, exp mode only, where it must lie in (0, 1]: the rule is the AccTrace averager's, whose create refuses anything else, so this does too. Ignored, and so not checked, in the other modes.
• hop Samples between frame starts, 1..n. 0 means n, which gives non-overlapping frames; the default. A hop above n would leave samples between frames that are never kept, so it is refused.

Returns:

Heap-allocated state, or NULL on invalid argument or OOM. A window whose coherent gain sum(w) is 0 Hann at n = 2, whose symmetric form is all zeros is refused too: every reading divides by it. So is an n or n * pad too large for its buffers to be sized, rather than wrapped to a small one.

Note:

Caller must call dp_psd_destroy() when done.

>>> from doppler.spectral import PSD
>>> w = PSD(n=1024, fs=1.0e6, window="kaiser", beta=8.0, mode="mean", hop=512)
>>> w.n, w.fs
(1024, 1000000.0)
>>> round(w.rbw / (w.fs / w.n), 3) == round(w.enbw, 3)
True

function dp_psd_destroy¶

Destroy a PSD instance and release all memory.

void dp_psd_destroy (
    dp_psd_state_t * state
) 

Parameters:

• state May be NULL (no-op).

function dp_psd_frame_db¶

One frame in dBFS, against the estimator's own reference.

void dp_psd_frame_db (
    dp_psd_kernel_t * k,
    const float _Complex * x,
    float * out
) 

dp_psd_frame_power() divided by the window coherent gain squared and the full-scale reference, as 10*log10 with the -200 dB floor, so a full-scale tone on a bin reads 0 dB whatever the window. The result equals dp_psd_psd_db() after accumulating that one frame, bit for bit (pinned by test_psd_core.c). Does not touch a PSD's running average.

Parameters:

• k A kernel dp_psd_kernel_init() filled.
• x Exactly k->n complex samples (cf32).
• out Output, k->nfft floats, DC-centred.

dp_psd_kernel_t k;
if (dp_psd_kernel_init (&k, 8, DP_PSD_WINDOW_RECT, 0.0f, 1, 1.0) != DP_OK)
  return 1;
float _Complex x[8];
float db[8];
for (int i = 0; i < 8; i++)
  x[i] = cexpf (I * 2.0f * 3.14159265f * 2.0f * (float)i / 8.0f);
dp_psd_frame_db (&k, x, db);       // unit tone on bin 2, rectangular
const int bad = fabsf (db[4 + 2]) > 1e-4f; // bin 2 is at nfft/2 + 2
dp_psd_kernel_free (&k);           // freed on either path
if (bad)
  return 1;

function dp_psd_frame_linear¶

One frame in linear full-scale^2 units, against the estimator's own reference.

void dp_psd_frame_linear (
    dp_psd_kernel_t * k,
    const float _Complex * x,
    float * out
) 

dp_psd_frame_power() divided by the window coherent gain squared and the full-scale reference squared, cg^2 * full_scale^2: a full-scale tone on a bin reads 1.0 whatever the window. dp_psd_frame_db() is 10*log10 of the same quotient, with the -200 dB floor the two are read from one double, so they cannot disagree about the reference. This is the frame a consumer that averages for itself should take: the raw dp_psd_frame_power() sits 20*log10(sum(w)) above it, 54.18 dB for Hann at 1024 points (docs/design/spectrogram-measurements.md section 5.3). Does not touch a PSD's running average; not thread-safe on one kernel.

Parameters:

• k A kernel dp_psd_kernel_init() filled.
• x Exactly k->n complex samples (cf32).
• out Output, k->nfft floats, DC-centred.

dp_psd_kernel_t k;
if (dp_psd_kernel_init (&k, 8, DP_PSD_WINDOW_HANN, 0.0f, 1, 1.0) != DP_OK)
  return 1;
float _Complex x[8];
float lin[8];
for (int i = 0; i < 8; i++)
  x[i] = cexpf (I * 2.0f * 3.14159265f * 2.0f * (float)i / 8.0f);
dp_psd_frame_linear (&k, x, lin);  // unit tone on bin 2, Hann
const int bad = fabsf (lin[4 + 2] - 1.0f) > 1e-4f; // bin 2 at nfft/2 + 2
dp_psd_kernel_free (&k);           // freed on either path
if (bad)
  return 1;

function dp_psd_frame_power¶

One frame's DC-centred two-sided power: window, zero-pad, FFT, |X|^2.

void dp_psd_frame_power (
    dp_psd_kernel_t * k,
    const float _Complex * x,
    float * pwr
) 

The per-frame kernel dp_psd_accumulate() folds into its average, exposed without the average. The result is the UN-NORMALISED power the averager folds: no division by the window coherent gain squared, so a tone does not read at its true level (dp_psd_frame_db() applies that division). Do not compare this raw output with dp_psd_power_twosided(), which divides by cg^2 on the double path. It stays the raw |X|^2 because that is what the averager folds; a consumer that wants a frame against the dBFS reference in linear units to average itself, or in dB takes dp_psd_frame_linear() or dp_psd_frame_db(). Uses the kernel's scratch, so it is not thread-safe on one kernel; a PSD's kernel is &psd->k, and this does not touch that PSD's running average.

Parameters:

• k A kernel dp_psd_kernel_init() filled.
• x Exactly k->n complex samples (cf32).
• pwr Output, k->nfft floats, bin k (centred) at index nfft/2 + k.

function dp_psd_get_state¶

void dp_psd_get_state (
    const dp_psd_state_t * state,
    void * blob
) 

function dp_psd_kernel_free¶

Release what dp_psd_kernel_init() allocated, and zero the kernel.

void dp_psd_kernel_free (
    dp_psd_kernel_t * k
) 

Parameters:

• k A kernel dp_psd_kernel_init() was called on (any outcome), a zeroed one, or NULL, which does nothing (as every destroy).

function dp_psd_kernel_init¶

Fill a kernel: the window, its gains, the scratch and the plan.

int dp_psd_kernel_init (
    dp_psd_kernel_t * k,
    size_t n,
    int window,
    float beta,
    size_t pad,
    double full_scale
) 

k is zeroed first, so it is safe to call dp_psd_kernel_free() on it whatever this returns. On a refusal or an allocation failure nothing is left allocated.

Parameters:

• k The kernel to fill (caller-owned storage).
• n Frame length in samples; at least 2.
• window DP_PSD_WINDOW_HANN, _KAISER, _BLACKMAN_HARRIS or _RECT.
• beta Kaiser beta (ignored by every other window). One that makes the window non-finite is refused.
• pad Zero-pad factor (>= 1); nfft = next_pow_two(n * pad).
• full_scale Amplitude that reads 0 dBFS; finite and > 0.

Returns:

DP_OK; DP_ERR_INVALID for an argument out of range, a window with no gain at n (Hann at 2) or an n * pad too large to size; DP_ERR_MEMORY when an allocation fails.

function dp_psd_noise_floor¶

Noise-floor estimate: median of the averaged dB spectrum.

double dp_psd_noise_floor (
    dp_psd_state_t * state
) 

Returns:

Median dB level (0 if empty).

function dp_psd_occupied_bw¶

Occupied bandwidth in Hz holding fraction of the total power.

double dp_psd_occupied_bw (
    dp_psd_state_t * state,
    double fraction
) 

Calls obw_from_power on the averaged spectrum, so the search and its domain are that function's and the two cannot disagree.

Parameters:

• state Must be non-NULL.
• fraction Power fraction in the OPEN interval (0, 1), e.g. 0.99. 1 is excluded: all of the power in float is every bin that rounding residue reached, which has no stable answer.

Returns:

Occupied bandwidth in Hz (0 if empty or no power); NaN when fraction is outside (0, 1), before any frame too, rather than a silently clamped answer.

>>> import numpy as np
>>> from doppler.spectral import PSD
>>> w = PSD(n=64, fs=64.0, window="rect")
>>> w.accumulate(np.ones(64, dtype=np.complex64))  # all power at DC
>>> w.occupied_bw(0.99)                            # one bin, fs / 64
1.0
>>> w.occupied_bw(1.0)                             # outside (0, 1)
nan

function dp_psd_power_onesided¶

Averaged linear power, one-sided (length nfft/2 + 1). Folds the DC-centred two-sided estimate onto [0, fs/2] : the DC and Nyquist bins are kept as-is, every interior bin is the sum of its +k and -k halves (so a real-input tone reads 2*avg|X[k] |^2 / cg^2 there). Coherent-gain normalised; full_scale is NOT applied. Returns 0 before any accumulate.

size_t dp_psd_power_onesided (
    dp_psd_state_t * state,
    float * out,
    size_t max_out
) 

Parameters:

• state Must be non-NULL.
• out Destination, at least nfft/2 + 1 float32 elements.
• max_out Capacity of out in elements. Emission stops there, so the return value is the number actually written.

Returns:

min(nfft/2 + 1, max_out), or 0 if empty.

function dp_psd_power_onesided_max_out¶

Output capacity hint for dp_psd_power_onesided() ; equals nfft/2+1.

size_t dp_psd_power_onesided_max_out (
    dp_psd_state_t * state
) 

function dp_psd_power_twosided¶

Averaged linear power, DC-centred two-sided (length nfft). Coherent-gain normalised ( out[k] = avg|X[k] |^2 / cg^2); full_scale is NOT applied (callers that want a dBFS reference divide by full_scale^2). This is the raw spectral estimate the measurement kernels integrate over. Returns 0 (and writes nothing) before any frame is accumulated.

size_t dp_psd_power_twosided (
    dp_psd_state_t * state,
    float * out,
    size_t max_out
) 

Parameters:

• state Must be non-NULL.
• out Destination, at least nfft float32 elements.
• max_out Capacity of out in elements. Emission stops there, so the return value is the number actually written.

Returns:

min(nfft, max_out), or 0 if empty.

function dp_psd_power_twosided_max_out¶

Output capacity hint for dp_psd_power_twosided() ; equals nfft.

size_t dp_psd_power_twosided_max_out (
    dp_psd_state_t * state
) 

function dp_psd_psd_db¶

Averaged power spectrum in dB, DC-centred. Normalised by window coherent gain so a full-scale tone reads its true power (a unit-amplitude tone peaks near 0 dB). Returns 0 (Python None) before any frame is accumulated.

size_t dp_psd_psd_db (
    dp_psd_state_t * state,
    float * out,
    size_t max_out
) 

Parameters:

• state Must be non-NULL.
• out Destination, at least nfft float32 elements.
• max_out Capacity of out in elements. Emission stops there, so the return value is the number actually written.

Returns:

min(nfft, max_out), or 0 if empty.

function dp_psd_psd_db_max_out¶

Output capacity hint for psd_db(); equals nfft.

size_t dp_psd_psd_db_max_out (
    dp_psd_state_t * state
) 

function dp_psd_psd_dbhz¶

Averaged power spectral density in dB/Hz, DC-centred. Normalised by fs * sum(w^2) (ENBW-aware), the standard one-sided-free PSD scaling. Differs from psd_db() by the constant10*log10(cg^2/(fs*s2)) . Returns 0 (Python None) before any frame is accumulated.

size_t dp_psd_psd_dbhz (
    dp_psd_state_t * state,
    float * out,
    size_t max_out
) 

>>> import numpy as np
>>> from doppler.spectral import PSD
>>> w = PSD(n=32, fs=2.0, window="hann", mode="mean")
>>> w.accumulate(np.ones(32, dtype=np.complex64))
>>> a = w.psd_db(); b = w.psd_dbhz()
>>> bool(np.allclose(a - b, (a - b)[0]))   # offset is a constant
True

function dp_psd_psd_dbhz_max_out¶

Output capacity hint for psd_dbhz(); equals nfft.

size_t dp_psd_psd_dbhz_max_out (
    dp_psd_state_t * state
) 

function dp_psd_reset¶

Discard the running average and any carried partial frame; the next accumulate re-seeds the average and starts the stream again.

void dp_psd_reset (
    dp_psd_state_t * state
) 

Parameters:

• state Must be non-NULL.

function dp_psd_set_state¶

int dp_psd_set_state (
    dp_psd_state_t * state,
    const void * blob
) 

function dp_psd_sfdr¶

Spurious-free dynamic range in dB from the two strongest peaks.

double dp_psd_sfdr (
    dp_psd_state_t * state,
    float min_db
) 

Parameters:

• state Must be non-NULL.
• min_db Minimum peak level considered, dB.

Returns:

Carrier-minus-highest-spur level in dB (0 if fewer than two peaks).

function dp_psd_snr¶

In-band SNR in dB: peak level in [lo_hz, hi_hz] minus the noise floor.

double dp_psd_snr (
    dp_psd_state_t * state,
    double lo_hz,
    double hi_hz
) 

Parameters:

• state Must be non-NULL.
• lo_hz Band lower edge, Hz.
• hi_hz Band upper edge, Hz.

Returns:

SNR in dB (0 if empty).

function dp_psd_state_bytes¶

size_t dp_psd_state_bytes (
    const dp_psd_state_t * state
) 

function dp_psd_total_band_power¶

Total integrated power across all bands in dB.

double dp_psd_total_band_power (
    dp_psd_state_t * state,
    const double * bands,
    size_t bands_len
) 

Parameters:

• state Must be non-NULL.
• bands Flat [lo,hi,...] band edges, Hz.
• bands_len Number of edge values (2 * n_bands).

Returns:

Total band power in dB (dB floor if empty).

Macro Definition Documentation¶

define PSD_STATE_MAGIC¶

#define PSD_STATE_MAGIC `DP_FOURCC ('P','S','D',' ')`

define PSD_STATE_VERSION¶

#define PSD_STATE_VERSION `2u`

The documentation for this class was generated from the following file native/inc/doppler/psd/psd_core.h