Quick Start

This guide walks through the main DSPkit workflows on a simulated 2DOF structural system — the same test signal used throughout the example scripts.

The test signal

DSPkit ships with a built-in 2DOF spring-mass-damper simulator:

from dspkit._testing import generate_2dof, natural_frequencies_2dof

fs = 1000.0   # Hz
t, a1, a2 = generate_2dof(
    duration=60.0,
    fs=fs,
    noise_std=1.0,
    output="acceleration",   # "displacement" | "velocity" | "acceleration"
    seed=42,
)

fn1, fn2 = natural_frequencies_2dof()
# fn1 ≈ 8.6 Hz,  fn2 ≈ 20.8 Hz

Spectral analysis

import dspkit as dsp

# Single-sided FFT amplitude spectrum
freqs, amp = dsp.fft_spectrum(a1, fs, window="hann", scaling="amplitude")
dsp.plot_fft(freqs, amp, xlim=(0, 80), title="FFT spectrum")

# Welch PSD
freqs, Pxx = dsp.psd(a1, fs, nperseg=4096)
dsp.plot_psd(freqs, Pxx, xlim=(0, 80))

# Magnitude-squared coherence between channels
freqs, Cxy = dsp.coherence(a1, a2, fs, nperseg=4096)
dsp.plot_coherence(freqs, Cxy, xlim=(0, 80))

# Autocorrelation
lags, acf = dsp.autocorrelation(a1, fs=fs, normalize=True, max_lag=0.5)
dsp.plot_autocorrelation(lags, acf, n_samples=len(a1), xlabel="Lag [s]")

Choosing the segment length, and when to stop using Welch

nperseg can be wrong in two opposite directions, and only one of them is obvious. Too few averages leaves a noisy spectrum and an inflated coherence; too short a segment to resolve the sharpest peak smears it, which biases coherence down and fabricates residual power. segment_advice weighs both on your own record and says which — including "this record cannot do both", which is a statement about the data rather than the parameters.

a = dsp.segment_advice(x1, fs, nperseg=1024)
print(a["verdict"], "—", a["advice"])
# resolution bandwidth, averages, bias floor and the sharpest peak found

When the record is genuinely too short, an autoregressive spectrum fits the whole of it at once instead of averaging segments, and can separate peaks Welch cannot:

freqs, Pxx, info = dsp.ar_psd(x1, fs)          # order chosen by AIC
freqs, Pxx, info = dsp.ar_psd(x1, fs, order=30)
print(info["order"], info["reflection_stable"])

The order is part of the answer, not a tuning knob: too low merges peaks, too high invents them. Check any peak that only the AR estimate shows against psd before believing it.

What one signal cannot explain about another

r = dsp.error_spectrum(target, predictors, fs, nperseg=2048)
r["error_spectrum"]        # (1 - coherence) * target PSD, in the target's units
r["unexplained_fraction"]  # the same thing as one number

Coherence says where a model fails; only this says whether that matters. A band where coherence collapses costs nothing if there is no power in it, which is the usual case.

Filtering

# Zero-phase Butterworth bandpass around first mode
a1_bp = dsp.bandpass(a1, fs, low=fn1 - 3, high=fn1 + 3, order=4)

# Notch filter to remove 50 Hz mains hum
a1_notched = dsp.notch(a1, fs, freq=50.0, q=30.0)

# Decimate to 100 Hz with anti-aliasing
a1_dec, fs_dec = dsp.decimate(a1, fs, target_fs=100.0)

Time-frequency analysis

import numpy as np

# Short-Time Fourier Transform
f, t, Zxx = dsp.stft(a1[:4096], fs, nperseg=256)
dsp.plot_spectrogram(f, t, Zxx, ylim=(0, 80))

# CWT scalogram (analytic Morlet, FFT-based)
analysis_freqs = np.geomspace(2.0, fs / 4, num=100)
f, t, W = dsp.cwt_scalogram(a1[:4096], fs, freqs=analysis_freqs)
dsp.plot_scalogram(f, t, W, ylim=(2, 80))

# Wigner-Ville (O(N²) — keep N small)
f, t, WVD = dsp.wigner_ville(a1[:512], fs)
dsp.plot_wvd(f, t, WVD, ylim=(0, 80))

Hilbert transform

# All three attributes in one call
env, phase, fi = dsp.hilbert_attributes(a1_bp, fs)

import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(10, 5), constrained_layout=True)
dsp.plot_signal(t, a1_bp, ax=axes[0], envelope=env, title="Bandpass + envelope")
axes[1].plot(t, fi); axes[1].set_ylabel("Instantaneous freq [Hz]")
plt.show()

EMD and HHT

# Decompose into IMFs
imfs, residue = dsp.emd(a1, max_imfs=8)

# Hilbert-Huang Transform
envs, inst_freqs = dsp.hht(imfs, fs)

# Marginal spectrum (adaptive PSD analogue)
freq_bins, marginal = dsp.hht_marginal_spectrum(envs, inst_freqs, fs, n_bins=512)

Peak detection

# Detect peaks in a PSD
peak_freqs, peak_vals, prominences = dsp.find_peaks(
    freqs, Pxx, distance_hz=5.0, max_peaks=5,
)
dsp.plot_peaks(freqs, Pxx, peak_freqs, peak_vals, db=True, xlim=(0, 80))

# Bandwidth and Q-factor
pf, bw, Q = dsp.peak_bandwidth(freqs, Pxx, peak_freqs=peak_freqs)

# Harmonic detection
hf, hv, orders = dsp.find_harmonics(freqs, amp, fundamental=25.0, n_harmonics=5)

SHM indicators

# Scalar indicators
se = dsp.spectral_entropy(freqs, Pxx)  # 0 = tonal, 1 = white noise
k = dsp.kurtosis(a1)                   # excess kurtosis (0 for Gaussian)
s = dsp.skewness(a1)                   # 0 for symmetric

# Time-varying indicators
times, rms_vals = dsp.rms_variation(a1, fs, segment_duration=10.0)
times, dom_freqs = dsp.frequency_shift(a1, fs, segment_duration=10.0)
times, energies = dsp.energy_variation(a1, fs, segment_duration=10.0)

dsp.plot_indicators(times, rms_vals, title="RMS Variation", ylabel="RMS")

Multi-sensor analysis

data = np.vstack([a1, a2])  # shape (n_channels, N)

# Correlation and coherence matrices
R = dsp.correlation_matrix(data)
freqs, C = dsp.coherence_matrix(data, fs, nperseg=4096)

# PSD matrix (input to FDD)
freqs, G = dsp.psd_matrix(data, fs, nperseg=4096)

dsp.plot_correlation_matrix(R, labels=["Ch1", "Ch2"])

FDD — Frequency Domain Decomposition

# Full OMA workflow
freqs, S, U = dsp.fdd_svd(data, fs, nperseg=4096)
dsp.plot_singular_values(freqs, S, db=True, xlim=(0, 80))

# Pick peaks → natural frequencies
peak_freqs, peak_idx = dsp.fdd_peak_picking(freqs, S, distance_hz=5.0, max_peaks=2)

# Extract mode shapes
modes = dsp.fdd_mode_shapes(U, peak_idx)
dsp.plot_mode_shape(modes[0], sensor_labels=["Mass 1", "Mass 2"])

# EFDD damping estimation
zeta, fn = dsp.efdd_damping(freqs, S, U, peak_idx, fs)

Probability and joint statistics

# PDF estimation
xi, density = dsp.pdf_estimate(a1)
dsp.plot_pdf(xi, density, hist_data=a1)

# Joint distribution
xc, yc, H = dsp.joint_histogram(a1, a2, bins=60)
dsp.plot_joint_histogram(xc, yc, H, xlabel="Ch1", ylabel="Ch2")

# Covariance and Mahalanobis distance
cov = dsp.covariance_matrix(data)
distances = dsp.mahalanobis(data)

Embedding plots in your own figures

Every plot_* function accepts an ax parameter:

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

fig = plt.figure(figsize=(12, 8), constrained_layout=True)
gs = gridspec.GridSpec(2, 2, figure=fig)

dsp.plot_psd(freqs, Pxx, ax=fig.add_subplot(gs[0, 0]))
dsp.plot_coherence(freqs, Cxy, ax=fig.add_subplot(gs[0, 1]))
dsp.plot_spectrogram(f, t, Zxx, ax=fig.add_subplot(gs[1, :]))
plt.show()