ppsd-rs

Seismic PPSD (Probabilistic Power Spectral Density) computation in Rust, matching ObsPy output. Zero unsafe, zero C dependencies.

What is PPSD?

PPSD (McNamara & Buland, 2004) is the standard method for characterizing ambient seismic noise at a station. It computes power spectral density from continuous waveforms, removes the instrument response, and accumulates statistics over time. PPSD is used by IRIS, FDSN, BMKG, and most seismic networks for station quality monitoring.

This crate implements the exact same algorithm as ObsPy's PPSD, validated against ObsPy output within 0.5 dB tolerance. It uses stationxml-rs types directly for instrument response evaluation -- no intermediate types, no conversion boilerplate.

Features

• Full PPSD pipeline matching ObsPy's PPSD.__process step-by-step
• Welch's method (mlab.psd compatible) with cosine taper and linear detrend
• Instrument response evaluation from stationxml_rs::Response -- walks PAZ, FIR, and Coefficients stages automatically
• FIR filter support with DTFT evaluation, DC normalization (matching evalresp), and symmetry expansion (None/Even/Odd)
• Konno-Ohmachi smoothing with configurable bandwidth
• Period binning with octave-based averaging (ObsPy style)
• Direct stationxml-rs integration -- takes &Response from any FDSN StationXML or SC3ML inventory
• Zero unsafe -- pure Rust math, no FFI
• Zero C dependencies -- compiles anywhere rustc runs

Quick Start

Add to your Cargo.toml:

[dependencies]
ppsd-rs = "0.1"
stationxml-rs = "0.2"

Compute PSD for one segment

use ppsd_rs::{process_segment, setup_period_binning, welch_psd};

// Load inventory (auto-detects FDSN StationXML or SC3ML)
let inv = stationxml_rs::read_from_file("IA.JAGI.xml").unwrap();
let channel = &inv.networks[0].stations[0].channels[0];
let response = channel.response.as_ref().unwrap();
let sample_rate = channel.sample_rate;

// ObsPy-compatible parameters
let ppsd_length = 3600.0;  // 1-hour segments
let nfft = (ppsd_length * sample_rate) as usize;
let nlap = nfft / 2;  // 50% overlap

// Period binning (matching ObsPy defaults)
let freqs: Vec<f64> = (1..=nfft / 2)
    .map(|i| i as f64 * sample_rate / nfft as f64)
    .collect();
let psd_periods: Vec<f64> = freqs.iter().rev().map(|&f| 1.0 / f).collect();
let (bin_left, bin_right, _centers) =
    setup_period_binning(&psd_periods, 0.04, 0.3);

// Process one segment (3600s of data)
let samples: Vec<f64> = load_your_data(); // your data source
let result = process_segment(
    &samples, sample_rate, nfft, nlap,
    response, &psd_periods, &bin_left, &bin_right,
).unwrap();

// result is Vec<Option<f64>> -- dB values per period bin
for (i, db) in result.iter().enumerate() {
    if let Some(val) = db {
        println!("bin {i}: {val:.1} dB");
    }
}

Evaluate instrument response only

use ppsd_rs::eval_response;

let inv = stationxml_rs::read_from_file("station.xml").unwrap();
let response = inv.networks[0].stations[0].channels[0]
    .response.as_ref().unwrap();

let freqs: Vec<f64> = (1..=100).map(|i| i as f64 * 0.1).collect();
let h_squared = eval_response(response, &freqs).unwrap();
// h_squared[i] = |H(f_i)|² including all stages (PAZ + FIR + ADC)

API Overview

use ppsd_rs::{
    // Full pipeline
    process_segment,           // one PPSD segment → dB values
    // Response evaluation
    eval_response,             // stationxml_rs::Response → |H(f)|²
    // Individual steps
    welch_psd,                 // Welch's method (mlab.psd)
    cosine_taper,              // Tukey window
    fft_power,                 // one-sided power spectrum
    konno_ohmachi_smooth,      // spectral smoothing
    period_bin_average,        // octave binning
    setup_period_binning,      // bin edge computation
    // Types
    PsdConfig, PsdResult,      // configuration and output
    PpsdError, Result,         // error handling
};

Response Evaluation

eval_response() walks all stages in stationxml_rs::Response and computes the combined |H(f)|²:

FIR coefficients are DC-normalized to match evalresp behavior. Symmetry expansion handles None (all coefficients), Even (mirror), and Odd (mirror with negation).

Architecture

src/
  lib.rs           -- single-file library: types, response eval, PSD math, tests

Design Decisions

• Direct stationxml-rs coupling: takes &Response directly -- no intermediate InstrumentResponse type. Both crates share the same author; updates propagate naturally.
• Single file: the crate is focused enough that one lib.rs keeps everything discoverable without module navigation overhead.
• Result<T> over panic: response evaluation can fail (empty response, missing decimation, unsupported stage) -- errors are explicit, not hidden behind defaults.
• DC normalization: FIR coefficients are normalized to unit DC gain, matching evalresp's internal behavior. This is critical for correct dB output.
• &[f64] over stream: user provides pre-sliced sample arrays -- simple, composable, no hidden I/O or buffering.

TDD with ObsPy

Test vectors are generated by Python/ObsPy scripts, ensuring the Rust output matches ObsPy PPSD within 0.5 dB tolerance across all 7 pipeline stages:

cd pyscripts && uv run python -m pyscripts
cargo test

Development

cargo build                     # build
cargo test                      # all tests (requires test vectors)
cargo clippy -- -D warnings     # lint (strict)
cargo fmt -- --check            # format check
cargo doc --no-deps --open      # browse docs locally

References

• McNamara & Buland (2004) -- original PPSD method paper
• ObsPy PPSD -- Python reference implementation
• FDSN StationXML -- instrument response format specification
• evalresp -- IRIS response evaluation (FIR normalization behavior)

Sister Projects

• stationxml-rs -- FDSN StationXML and SC3ML reader/writer (provides Response types used by this crate)
• miniseed-rs -- miniSEED v2/v3 decoder and encoder
• seedlink-rs -- SeedLink protocol client/server

License

Apache-2.0