Signal processing

Signal processing#

The package provides tools for spectral analysis of seismic data. The main class is the ShortTimeFourierTransform class, which extends the ShortTimeFFT class to provide a more user-friendly interface for the Short-Time Fourier Transform with streams.

Other functions are provided to normalize seismic traces in the spectral domain, such as the modulus_division() and smooth_modulus_division() functions. Also, we provide the smooth_envelope_division() function to normalize seismic traces by a smooth version of its envelope.

class covseisnet.signal.ShortTimeFourierTransform(window_duration: float = 2.0, window_step: None | float = None, window_function: str = 'hann', sampling_rate: float = 1.0, **kwargs)[source]#

Bases: ShortTimeFFT

Short-Time Fourier Transform instance.

This class extends the scipy.signal.ShortTimeFFT class to provide a more user-friendly interface for the Short-Time Fourier Transform with ObsPy and Covseisnet Streams.

Parameters:

  • window_duration (float, optional) -- The duration of the window in seconds. Default is 2.0 seconds.

  • window_step (float, optional) -- The step between windows in seconds. Default is half the window duration.

  • window_function (str, optional) -- The window function to use. Default is the Hann window.

  • sampling_rate (float, optional) -- The sampling rate of the data. Default is 1.0 Hz. In most methods the sampling rate is set automatically from the data.

  • **kwargs (dict, optional) -- Additional keyword arguments are passed to the scipy.signal.ShortTimeFFT.stft() method.

Notes

By default, the scipy.signal.ShortTimeFFT class pads the data with zeros on both sides to avoid edge effects, and enable the analysis of edge samples (see the tutorial on the Short-Time Fourier Transform in the SciPy documentation). Because the main purpose of this class is to analyse the spatial coherence of seismic data, this default behaviour is disabled. Indeed, the edges appear wrongly coherent because of the zero-padding. This is why the representation of the Short-Time Fourier Transform indicates the number of samples and windows that are out-of-bounds. These samples are discarded from the analysis. Note that this is not a real issue since the covariance matrix analysis is performed on sliding windows.

Examples

>>> from covseisnet.signal import ShortTimeFourierTransform
>>> stft = ShortTimeFourierTransform(
...     window_duration=2.0,
...     window_step=1.0,
...     window_function="hann",
...     sampling_rate=100.0,
... )
>>> print(stft)
Short-Time Fourier Transform instance
    Sampling rate: 100.0 Hz
    Frequency range: 0.50 to 50.00 Hz
    Frequency resolution: 0.5 Hz
    Frequency bins: 101
    Window length: 200 samples
    Window step: 100 samples
    Out-of-bounds: 101 sample(s), 1 window(s)

See also

scipy.signal.ShortTimeFFT

property sampling_rate: float#

The sampling rate of the Short-Time Fourier Transform.

This property returns the sampling rate of the Short-Time Fourier Transform in Hertz.

Returns:

float -- The sampling rate of the Short-Time Fourier Transform in Hertz.

property frequencies: ndarray#

The frequencies of the Short-Time Fourier Transform.

This property returns the frequencies of the Short-Time Fourier Transform in Hertz. The frequencies are calculated from the sampling rate and the number of frequency bins.

Returns:

numpy.ndarray -- The frequencies of the Short-Time Fourier Transform in Hertz.

times(stats: Stats, **kwargs) → ndarray[source]#

The times of the window centers in seconds.

This method returns the times of the window centers in seconds. The times are calculated from the start time of the trace and the window step.

Parameters:

stats (Stats) -- The stats object of the trace.

Returns:

numpy.ndarray -- The times of the window centers in matplotlib datenum format.

valid_window_bounds(stats: Stats, **kwargs) → dict[source]#

The bounds of the Short-Time Fourier Transform.

This method calculate the valid window bounds of the Short-Time Fourier Transform in seconds. The bounds are calculated from the start time of the trace, the window duration, and the window step.

Parameters:

stats (Stats) -- The stats object of the trace.

Returns:

tuple -- The lower and upper bounds of the Short-Time Fourier Transform in seconds.

transform(trace: Trace, detrend: Literal['linear', 'constant'] = 'linear', **kwargs) → tuple[source]#

Short-time Fourier Transform of a trace.

This method calculates the Short-Time Fourier Transform of a trace using the scipy.signal.ShortTimeFFT.stft_detrend() method. The method returns the times of the window centers in matplotlib datenum format, the frequencies of the spectrogram, and the short-time spectra of the trace.

Prior to the Short-Time Fourier Transform, the method removes the linear trend from the trace to avoid edge effects. The method also discards the out-of-bounds samples and windows that are padded with zeros by the scipy.signal.ShortTimeFFT.

Parameters:
Returns:

  • times (numpy.ndarray) -- The times of the window centers in matplotlib datenum format.

  • frequencies (numpy.ndarray) -- The frequencies of the spectrogram.

  • short_time_spectra (numpy.ndarray) -- The short-time spectra of the trace with shape (n_frequencies, n_times).

Examples

>>> import covseisnet as csn
>>> stft = csn.ShortTimeFourierTransform(
...     window_duration=2.0,
...     window_step=1.0,
...     window_function="hann",
...     sampling_rate=100.0,
... )
>>> stream = csn.read()
>>> trace = stream[0]
>>> times, frequencies, short_time_spectra = stft.transform(trace)
>>> print(times.shape, frequencies.shape, short_time_spectra.shape)
(28,) (101,) (101, 28)
map_transform(stream: Stream, **kwargs) → tuple[source]#

Transform a stream into the spectral domain.

This method transforms a stream into the spectral domain by applying the Short-Time Fourier Transform to each trace of the stream. The method returns the times of the spectrogram in matplotlib datenum format, the frequencies of the spectrogram, and the spectrogram of the stream.

This method is basically a wrapper around the transform() method that applies the Short-Time Fourier Transform to each trace of the stream.

Note that the stream must be ready for processing, i.e., it must pass the quality control from the property is_ready_to_process.

Parameters:

stream (Stream) -- The stream to transform.

Returns:

  • times (numpy.ndarray) -- The times of the spectrogram in matplotlib datenum format.

  • frequencies (numpy.ndarray) -- The frequencies of the spectrogram.

  • short_time_spectra (numpy.ndarray) -- The spectrogram of the stream with shape (n_trace, n_frequencies, n_times).

  • **kwargs (dict, optional) -- Additional keyword arguments are passed to the transform() method.

Examples

>>> import covseisnet as csn
>>> stft = csn.ShortTimeFourierTransform(
...     window_duration=2.0,
...     window_step=1.0,
...     window_function="hann",
...     sampling_rate=100.0,
... )
>>> stream = csn.read()
>>> times, frequencies, short_time_spectra = stft.map_transform(stream)
>>> print(times.shape, frequencies.shape, short_time_spectra.shape)
(28,) (101,) (3, 101, 28)
covseisnet.signal.modulus_division(x: int | float | complex | ndarray, epsilon=1e-10) → ndarray[source]#

Division of a number by the absolute value or its modulus.

Given a complex number (or array) \(x = a e^{i\phi}\), where \(a\) is the modulus and \(\phi\) the phase, the function returns the unit-modulus complex number such as

\[\mathbb{C} \ni \tilde{x} = \frac{x}{|x| + \epsilon} \approx e^{i\phi}\]

This method normalizes the input complex number by dividing it by its modulus plus a small epsilon value to prevent division by zero, effectively scaling the modulus to 1 while preserving the phase.

Note that this function also work with real-valued arrays, in which case the modulus is the absolute value of the real number. It is useful to normalize seismic traces in the temporal domain. It writes

\[\mathbb{R} \ni \tilde{x} = \frac{x}{|x| + \epsilon} \approx \text{sign}(x)\]
Parameters:

  • x (numpy.ndarray) -- The complex-valued data to extract the unit complex number from.

  • epsilon (float, optional) -- A small value added to the modulus to avoid division by zero. Default is 1e-10.

Returns:

numpy.ndarray -- The unit complex number with the same phase as the input data, or the sign of the real number if the input is real-valued.

covseisnet.signal.smooth_modulus_division(x: ndarray, smooth: None | int = None, order: int = 1, epsilon: float = 1e-10) → ndarray[source]#

Modulus division of a complex number with smoothing.

Given a complex array \(x[n] = a[n] e^{i\phi[n]}\), where \(a[n]\) is the modulus and \(\phi[n]\) the phase, the function returns the normalized complex array \(\tilde{x}[n]\) such as

\[\tilde{x}[n] = \frac{x[n]}{\mathcal S a[n] + \epsilon}\]

where \(\mathcal S a[n]\) is a smoothed version of the modulus array \(a[n]\). The smoothing function is performed with the Savitzky-Golay filter, and the order and length of the filter are set by the smooth and order parameters.

Parameters:

  • x (numpy.ndarray) -- The spectra to detrend. Must be of shape (n_frequencies, n_times).

  • smooth (int) -- Smoothing window size in points.

  • order (int) -- Smoothing order. Check the scipy function savgol_filter() function for more details.

Keyword Arguments:

epsilon (float, optional) -- A regularizer for avoiding zero division.

Returns:

The spectrum divided by the smooth modulus spectrum.

covseisnet.signal.smooth_envelope_division(x: ndarray, smooth: int, order: int = 1, epsilon: float = 1e-10) → ndarray[source]#

Normalize seismic traces by a smooth version of its envelope.

This function normalizes seismic traces by a smooth version of its envelope. The envelope is calculated with the Hilbert transform, and then smoothed with the Savitzky-Golay filter. The order and length of the filter are set by the smooth and order parameters.

Considering the seismic trace \(x(t)\), the normalized trace \(\hat x(t)\) is obtained with

\[\hat x(t) = \frac{x(t)}{\mathcal{A}x(t) + \epsilon}\]

where \(A\) is an smoothing operator applied to Hilbert envelope of the trace \(x(t)\).

Parameters:

  • x (numpy.ndarray) -- The seismic trace to normalize.

  • smooth (int) -- The length of the Savitzky-Golay filter for smoothing the envelope.

  • order (int) -- The order of the Savitzky-Golay filter for smoothing the envelope.

  • epsilon (float, optional) -- Regularization parameter in division, set to 1e-10 by default.

Returns:

numpy.ndarray -- The normalized seismic trace.

covseisnet.signal.entropy(x: ndarray, epsilon: float = 1e-10, axis: int = -1) → ndarray[source]#

Entropy calculation.

Entropy calculated from a given distribution of values. The entropy is defined as

\[h = - \sum_{n=0}^N n x_n \log(x_n + \epsilon)\]

where \(x_n\) is the distribution of values. This function assumes the distribution is normalized by its sum.

Parameters:

  • x (numpy.ndarray) -- The distribution of values.

  • epsilon (float, optional) -- The regularization parameter for the logarithm.

  • **kwargs (dict, optional) -- Additional keyword arguments passed to the numpy.sum() function. Typically, the axis along which the sum is performed.

covseisnet.signal.diversity(x: ndarray, epsilon: float = 1e-10, **kwargs) → ndarray[source]#

Shanon diversity index calculation.

Shanon diversity index calculated from a given distribution of values. The diversity index is defined as

\[D = \exp(h + \epsilon)\]

where \(h\) is the entropy of the distribution of values. This function assumes the distribution is normalized by its sum.

Parameters:

  • x (numpy.ndarray) -- The distribution of values.

  • epsilon (float, optional) -- The regularization parameter for the logarithm.

  • **kwargs (dict, optional) -- Additional keyword arguments passed to the numpy.sum() function. Typically, the axis along which the sum is performed.

covseisnet.signal.width(x: ndarray, **kwargs) → ndarray[source]#

Width calculation.

Width calculated from a given distribution of values. The width is defined as

\[\sigma = \sum_{n=0}^N n x_n\]

where \(x_n\) is the distribution of values. This function assumes the distribution is normalized by its sum.

Parameters:

  • x (numpy.ndarray) -- The distribution of values.

  • **kwargs (dict, optional) -- Additional keyword arguments passed to the numpy.sum() function. Typically, the axis along which the sum is performed.

covseisnet.signal.bandpass_filter(x, sampling_rate, frequency_band, filter_order=4)[source]#

Bandpass filter the signal.

Apply a Butterworth bandpass filter to the signal. Uses butter() and filtfilt() to avoid phase shift.

Parameters:

  • x (ndarray) -- The signal to filter.

  • sampling_rate (float) -- The sampling rate in Hz.

  • frequency_band (tuple) -- The frequency band to filter in Hz.

  • filter_order (int, optional) -- The order of the Butterworth filter.

Returns:

ndarray -- The filtered signal.

covseisnet.signal.hilbert_envelope(x: ndarray, **kwargs) → ndarray[source]#

Hilbert envelope calculation.

Calculate the envelope of a signal using the Hilbert transform. The envelope is defined as

\[\text{envelope}(x) = \sqrt{x^2 + \text{hilbert}(x)^2}\]
Parameters:

x (numpy.ndarray) -- The signal to calculate the envelope from.

Returns:

numpy.ndarray -- The envelope of the signal.