Signal Resampling and Downsampling in Seismology
Resampling changes the sampling rate of a signal, either by increasing it (upsampling) or decreasing it (downsampling). In seismology, downsampling is especially important when:
- Reducing file sizes for long time series
- Matching sampling rates between stations
- Preparing data for band-limited analyses
Why Be Careful with Downsampling?
Downsampling without preparation can lead to aliasing, where high-frequency energy folds back into lower frequencies and distorts the signal.
Key Steps to Proper Downsampling:
- Detrend the signal to remove DC offsets.
- Apply a lowpass filter to eliminate frequencies above the new Nyquist frequency.
- Resample in small steps, especially when reducing by large factors (e.g., 100 Hz → 1 Hz).
Recommended Downsampling Strategy
def safe_downsample(trace, target_rate, max_factor=10, pre_filter=True):
"""
Downsample a trace to target_rate in steps, applying filtering.
Parameters:
- trace: ObsPy Trace object
- target_rate: Final desired sampling rate (Hz)
- max_factor: Maximum downsample ratio per step (e.g., 10)
- pre_filter: Whether to apply lowpass filter before each step
"""
from obspy import Trace
import copy
tr = trace.copy()
while tr.stats.sampling_rate > target_rate:
current_rate = tr.stats.sampling_rate
factor = min(current_rate / target_rate, max_factor)
new_rate = current_rate / factor
tr.detrend("demean")
tr.taper(max_percentage=0.05)
if pre_filter:
f_nyquist = new_rate * 0.5
tr.filter("lowpass", freq=0.9 * f_nyquist, corners=4, zerophase=True)
tr.resample(new_rate)
return tr
Mathematical Concept of Fourier-Based Resampling
Fourier-based resampling operates in the frequency domain using the following key steps:
-
Compute the FFT of the time-domain signal: $$ X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-2\pi i \cdot kn/N} $$
-
Interpolate or zero-pad/truncate the spectrum to match the new length \( N' \), corresponding to the desired sampling rate.
-
Optionally apply a window \( W[k] \) (e.g., Hann) to reduce spectral leakage.
-
Perform the inverse FFT to return to the time domain: $$ x'[n] = \text{IFFT}(X'[k]) $$
Note: This method assumes periodicity of the signal and is well-suited for smoothly varying signals.
Adjusting Spectrum Length for Resampling
To change the sampling rate using Fourier methods, we adjust the number of frequency components before applying the inverse FFT. This process depends on whether we are downsampling or upsampling.
🔽 Downsampling: Truncate the Spectrum
When reducing the sampling rate, the new Nyquist frequency is lower. We must truncate high-frequency components that the new sampling rate cannot represent.
How?
# Assume X is the FFT of the original signal
X_truncated = X[:N_new // 2 + 1]
🔼 Upsampling: Zero-Pad the Spectrum
When increasing the sampling rate, the signal must contain more samples, so we insert zeros in the frequency domain.
How?
# Assume X is the FFT of the original signal
pad_length = N_new // 2 + 1 - len(X)
X_padded = np.pad(X, (0, pad_length), 'constant')
Important Note on Symmetry
When working with real signals and rfft
, ensure the spectrum retains Hermitian symmetry. Both truncation and padding should only be done on the one-sided spectrum and handled carefully when converting back to time domain with irfft
.
Summary Table
Step | Purpose |
---|---|
Detrend | Removes constant offset |
Lowpass Filter | Prevents aliasing by cutting high frequencies |
Resample | Adjusts time step and signal length |
Multi-Step | Ensures safe and smooth reduction |
FFT Adjustments | Truncate or pad spectrum appropriately |
This strategy ensures high-quality downsampling for long-duration or broadband seismic signals, protecting both time and frequency content.