Digital signal processing with scipy.signal: filter design, spectral analysis, peak detection, PSD estimation, and adaptive filtering for engineering workflows.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
Digital signal processing with scipy.signal: filter design, spectral analysis, peak detection, PSD estimation, and adaptive filtering for engineering workflows.
A comprehensive skill for digital signal processing (DSP) using scipy.signal.
Covers FIR/IIR filter design, time-series filtering, spectral analysis, peak
detection, power spectral density estimation, cross-correlation, matched filters,
and adaptive filtering. Designed for real-world applications such as biomedical
signal processing and vibration/fault analysis.
Core Functions
1. Filter Design
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
defdesign_bandpass_filter(
lowcut: float,
highcut: float,
fs: float,
order: int = 4,
filter_type: str = "butter",
rp: float = 1.0,
rs: float = 60.0,
) -> tuple[np.ndarray, np.ndarray]:
"""
Design a bandpass IIR or FIR filter.
Parameters
----------
lowcut : float
Lower cutoff frequency in Hz.
highcut : float
Upper cutoff frequency in Hz.
fs : float
Sampling frequency in Hz.
order : int
Filter order (default 4).
filter_type : str
One of 'butter', 'cheby1', 'cheby2', 'ellip', 'firwin'.
rp : float
Maximum ripple in the passband (dB), used by cheby1/ellip.
rs : float
Minimum attenuation in the stopband (dB), used by cheby2/ellip.
Returns
-------
b, a : array_like
Numerator and denominator polynomials of the IIR filter.
For FIR (firwin), a = [1.0].
"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
if filter_type == "butter":
b, a = signal.butter(order, [low, high], btype="bandpass")
elif filter_type == "cheby1":
b, a = signal.cheby1(order, rp, [low, high], btype="bandpass")
filter_type == :
b, a = signal.cheby2(order, rs, [low, high], btype=)
filter_type == :
b, a = signal.ellip(order, rp, rs, [low, high], btype=)
filter_type == :
order % != :
order +=
b = signal.firwin(order + , [low, high], pass_zero=)
a = np.array([])
:
ValueError()
b, a
() -> [np.ndarray, np.ndarray]:
w0 = freq / ( * fs)
b, a = signal.iirnotch(w0, quality)
b, a
() -> [np.ndarray, np.ndarray]:
nyq = * fs
normal_cutoff = cutoff / nyq
filter_type == :
b, a = signal.butter(order, normal_cutoff, btype=)
filter_type == :
order % == :
order +=
b = signal.firwin(order, normal_cutoff)
a = np.array([])
:
ValueError()
b, a
elif
"cheby2"
"bandpass"
elif
"ellip"
"bandpass"
elif
"firwin"
# FIR bandpass via window method; order must be even for bandpass
if
2
0
1
1
False
1.0
else
raise
f"Unsupported filter_type: {filter_type!r}"
return
def
design_notch_filter
freq: float, fs: float, quality: float = 30.0
tuple
"""
Design a notch (band-stop) filter to remove a single frequency.
Parameters
----------
freq : float
Frequency to attenuate in Hz.
fs : float
Sampling frequency in Hz.
quality : float
Quality factor Q = freq / bandwidth.
Returns
-------
b, a : array_like
"""
defapply_filter(
data: np.ndarray,
b: np.ndarray,
a: np.ndarray,
method: str = "sosfilt",
) -> np.ndarray:
"""
Apply a digital filter to a 1-D signal.
Parameters
----------
data : np.ndarray
Input signal array (1-D).
b : np.ndarray
Numerator coefficients.
a : np.ndarray
Denominator coefficients.
method : str
'sosfilt' (recommended, numerically stable) or 'filtfilt' (zero-phase).
Returns
-------
np.ndarray
Filtered signal with the same length as `data`.
"""if method == "sosfilt":
sos = signal.tf2sos(b, a)
return signal.sosfiltfilt(sos, data)
elif method == "filtfilt":
return signal.filtfilt(b, a, data)
elif method == "lfilter":
return signal.lfilter(b, a, data)
else:
raise ValueError(f"Unknown method: {method!r}")
defapply_envelope_detection(
data: np.ndarray, fs: float, lowpass_cutoff: float = 10.0) -> np.ndarray:
"""
Compute the amplitude envelope using the Hilbert transform followed by
a lowpass filter.
Parameters
----------
data : np.ndarray
Bandpass-filtered signal.
fs : float
Sampling frequency in Hz.
lowpass_cutoff : float
Cutoff for the smoothing lowpass filter (Hz).
Returns
-------
np.ndarray
Smoothed amplitude envelope.
"""
analytic = signal.hilbert(data)
envelope = np.abs(analytic)
b, a = design_lowpass_filter(lowpass_cutoff, fs, order=4)
return apply_filter(envelope, b, a)
3. Spectral Analysis
defcompute_spectrogram(
data: np.ndarray,
fs: float,
window: str = "hann",
nperseg: int = 256,
noverlap: int | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Compute the Short-Time Fourier Transform (STFT) spectrogram.
Parameters
----------
data : np.ndarray
Input time-series signal.
fs : float
Sampling frequency in Hz.
window : str
Window function name (e.g. 'hann', 'hamming', 'blackman').
nperseg : int
Length of each STFT segment.
noverlap : int or None
Number of overlapping samples. Defaults to nperseg // 2.
Returns
-------
f : np.ndarray
Frequency bins (Hz).
t : np.ndarray
Time bins (seconds).
Sxx : np.ndarray
Power spectrogram in dB (shape: [freq_bins, time_bins]).
"""if noverlap isNone:
noverlap = nperseg // 2
f, t, Zxx = signal.stft(data, fs=fs, window=window, nperseg=nperseg, noverlap=noverlap)
Sxx = 20 * np.log10(np.abs(Zxx) + 1e-12)
return f, t, Sxx
defcompute_psd_welch(
data: np.ndarray,
fs: float,
nperseg: int = 512,
window: str = "hann",
) -> tuple[np.ndarray, np.ndarray]:
"""
Estimate power spectral density using Welch's method.
Parameters
----------
data : np.ndarray
Input signal.
fs : float
Sampling frequency in Hz.
nperseg : int
Length of each Welch segment.
window : str
Window function name.
Returns
-------
freqs : np.ndarray
Frequency bins (Hz).
psd : np.ndarray
Power spectral density (V²/Hz).
"""
freqs, psd = signal.welch(data, fs=fs, window=window, nperseg=nperseg)
return freqs, psd
defplot_frequency_response(
b: np.ndarray,
a: np.ndarray,
fs: float,
title: str = "Filter Frequency Response",
) -> None:
"""Plot magnitude and phase response using freqz."""
w, h = signal.freqz(b, a, worN=8000, fs=fs)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))
ax1.plot(w, 20 * np.log10(np.abs(h) + 1e-12))
ax1.set_xlabel("Frequency (Hz)")
ax1.set_ylabel("Magnitude (dB)")
ax1.set_title(title)
ax1.grid(True)
ax2.plot(w, np.angle(h, deg=True))
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("Phase (degrees)")
ax2.grid(True)
plt.tight_layout()
plt.show()
4. Peak Detection and Correlation
defdetect_peaks(
data: np.ndarray,
min_height: float | None = None,
min_distance: int = 1,
prominence: float | None = None,
) -> tuple[np.ndarray, dict]:
"""
Detect peaks in a 1-D signal using scipy.signal.find_peaks.
Parameters
----------
data : np.ndarray
Input signal.
min_height : float or None
Minimum peak height.
min_distance : int
Minimum number of samples between peaks.
prominence : float or None
Minimum peak prominence.
Returns
-------
peaks : np.ndarray
Indices of detected peaks.
properties : dict
Peak properties (heights, prominences, etc.).
"""
kwargs: dict = {"distance": min_distance}
if min_height isnotNone:
kwargs["height"] = min_height
if prominence isnotNone:
kwargs["prominence"] = prominence
peaks, properties = signal.find_peaks(data, **kwargs)
return peaks, properties
defcompute_cross_correlation(
x: np.ndarray, y: np.ndarray, normalize: bool = True) -> tuple[np.ndarray, np.ndarray]:
"""
Compute full cross-correlation between two signals.
Returns
-------
lags : np.ndarray
Lag values in samples.
corr : np.ndarray
Cross-correlation coefficients.
"""
corr = signal.correlate(x, y, mode="full")
lags = signal.correlation_lags(len(x), len(y), mode="full")
if normalize:
corr = corr / (np.std(x) * np.std(y) * len(x))
return lags, corr
defmatched_filter(template: np.ndarray, noisy_signal: np.ndarray) -> np.ndarray:
"""
Apply a matched filter (cross-correlate signal with time-reversed template).
Returns
-------
np.ndarray
Matched filter output (same length as noisy_signal).
"""
h = template[::-1]
return signal.convolve(noisy_signal, h, mode="same")
5. Adaptive Filtering (LMS)
deflms_adaptive_filter(
desired: np.ndarray,
input_signal: np.ndarray,
filter_order: int = 32,
mu: float = 0.01,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Least Mean Squares (LMS) adaptive filter.
Parameters
----------
desired : np.ndarray
Desired (reference) signal.
input_signal : np.ndarray
Input (noisy) signal to filter.
filter_order : int
Number of filter taps.
mu : float
Step size (learning rate). Must be small for stability (< 1 / (filter_order * P_x)).
Returns
-------
output : np.ndarray
Filtered output signal.
error : np.ndarray
Error signal (desired - output).
weights : np.ndarray
Final filter weight vector.
"""
n = len(input_signal)
weights = np.zeros(filter_order)
output = np.zeros(n)
error = np.zeros(n)
for i inrange(filter_order, n):
x = input_signal[i - filter_order : i][::-1]
y = np.dot(weights, x)
e = desired[i] - y
weights += 2 * mu * e * x
output[i] = y
error[i] = e
return output, error, weights
Example 1: EMG Signal Processing Pipeline
This example simulates an electromyography (EMG) recording and runs a full
processing pipeline: bandpass filtering, power-line notch removal, and
amplitude envelope extraction.
Example 2: Vibration Analysis for Rotating Machinery
This example performs fault-frequency detection on a vibration signal from a
rotating machine. It uses FFT, spectrogram, and automated peak detection to
identify bearing defect frequencies.
Filter stability: Always use sosfilt / sosfiltfilt for high-order
filters. Direct form II (lfilter) can suffer from numerical instability for
orders above 6–8.
Zero-phase filtering: filtfilt applies the filter twice (forward and
backward), eliminating phase distortion but doubling the effective order.
Avoid for online/real-time applications.
Nyquist constraint: Cutoff frequencies must be strictly between 0 and
fs/2. Always validate inputs before calling design functions.
FIR vs IIR: FIR filters (firwin) are inherently stable and have linear
phase, but require higher orders to achieve sharp transitions. IIR filters
(butter, ellip) achieve steeper roll-off at lower order but introduce phase
distortion.
LMS convergence: The step size mu must satisfy 0 < mu < 1 / (M * Px)
where M is filter order and Px is the input signal power. Start small
(e.g., mu = 0.001) and increase cautiously.
Window selection for STFT: Hann window is a good default. Blackman-Harris
offers lower sidelobes for detecting weak tones near strong ones.
Welch PSD segments: Larger nperseg gives finer frequency resolution but
fewer averages (higher variance). Typical values: nperseg = fs (1-second
segments) for stationary signals.
