TopolearnSensor, audio and physiological data need noise removal and frequency analysis before modelling.
scipy.signal covers the classic digital signal processing toolkit; scipy.fft provides fast Fourier transforms.
Filter design: signal.butter(N, Wn, btype="low", fs=fs, output="sos") designs a Butterworth filter of order N with cutoff Wn in Hz when you pass the sampling rate fs. Ask for second-order sections (output="sos"): the older (b, a) form becomes numerically unstable at higher orders.
Applying filters:
signal.sosfilt(sos, x): causal, as a real-time system would, but delays the signal,signal.sosfiltfilt(sos, x): runs forward and backward for zero phase shift, ideal for offline analysis.
Peaks: signal.find_peaks(x, height=..., distance=..., prominence=...) returns indices and properties. Prominence is usually the most useful criterion for noisy data.
Spectra: scipy.fft.rfft(x) with rfftfreq(n, d=1/fs) gives the frequency content of real signals. For noisy signals, signal.welch(x, fs=fs, nperseg=...) averages over segments to estimate the power spectral density with far less variance; signal.spectrogram or ShortTimeFFT show how it changes over time.
Remember the Nyquist limit: you can only see frequencies below fs / 2. Resample with signal.resample_poly, which filters to avoid aliasing, rather than taking every n-th sample.
Example
import numpy as np
from scipy import signal
from scipy.fft import rfft, rfftfreq
fs = 500 # samples per second
t = np.arange(0, 2, 1 / fs)
rng = np.random.default_rng(0)
x = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 80 * t) + 0.2 * rng.normal(size=t.size)
freqs = rfftfreq(t.size, d=1 / fs)
top = np.argsort(np.abs(rfft(x)))[-2:]
print(sorted(freqs[top].tolist())) # [5.0, 80.0] (dominant frequencies)
sos = signal.butter(4, 20, btype="low", fs=fs, output="sos")
clean = signal.sosfiltfilt(sos, x) # zero-phase low-pass at 20 Hz
peaks, _ = signal.find_peaks(clean, prominence=1)
print(len(peaks), t[peaks[:3]].round(3)) # 9 [0.244 0.444 0.65 ] (every ~0.2 s = 5 Hz)
f, pxx = signal.welch(x, fs=fs, nperseg=256)
print(f[np.argmax(pxx)]) # 5.859375 (near 5 Hz; bins are fs/nperseg ~ 2 Hz)Remember
- Design filters with
output="sos"and passfsto use Hz sosfiltfiltfilters without phase shift;sosfiltis causal- Use
welchfor noisy spectra and respect the Nyquist limit
Watch out for
- Downsampling by slicing
x[::10]without an anti-aliasing filter.