The analysis of multi-feature extraction methodologies applied to seismic signal processing, based on the
implementation of research-validated techniques and analysis of experimental data comprising 56 non-overlapping
windows with 66 extracted features per window. The investigation encompasses temporal, spectral, and time-
frequency domain feature extraction following established methodologies for seismic data analysis, with particular
emphasis on the statistical significance and discriminative power of extracted features for machine learning applications
[3][4]
.
1. Introduction and Methodological Framework
1.1 Research Foundation and Scope
Seismic signal processing represents a critical subfield of digital signal processing that focuses on enhancing signals,
suppressing noise, and extracting meaningful features from ground motion data. The multi-feature extraction approach
implemented in this study follows established research methodologies that combine temporal statistical analysis,
spectral decomposition techniques, and time-frequency domain transformations to create comprehensive feature
representations suitable for automated classification systems.
The experimental dataset analyzed consists of 56 temporal windows of 1024 samples each, processed using
non-overlapping sliding window techniques as specified in seismic signal processing literature. This approach
ensures temporal independence between analysis segments while maintaining sufficient signal content for robust
feature extraction across multiple domains.
1.2 Feature Extraction Domains
The implemented methodology extracts features from three distinct signal processing domains, each providing
complementary information about seismic signal characteristics:
Temporal Domain (4 features): Statistical parameters including mean, standard deviation, skewness, and kurtosis
computed directly from detrended signal amplitudes.
Spectral Domain (22 features): Frequency-based characteristics derived through Welch power spectral density
estimation and Fast Fourier Transform analysis with frequency binning.
Time-Frequency Domain (40 features): Multi-resolution analysis using discrete wavelet transform with
Daubechies-4 wavelets and 5-level decomposition.
2. Temporal Domain Analysis
2.1 Statistical Parameter Characteristics
The temporal domain analysis reveals significant insights into the seismic signal amplitude distributions across the
56 analyzed windows. The statistical parameters demonstrate varying degrees of discriminative power, with standard
deviation showing the highest stability (CV = 0.401) among temporal features.
Mean Values: All temporal means equal zero due to proper signal detrending, confirming adherence to research
methodology specifications. This detrending removes DC bias and ensures subsequent statistical measures focus on
signal variations rather than absolute amplitude levels .
Standard Deviation Distribution: Values range from 1,432 to 8,449 with a mean of 4,161, representing the spread
of seismic amplitude around the mean. The coefficient of variation (0.401) indicates moderate variability across
windows, suggesting this parameter effectively discriminates between different seismic activity levels.
Skewness Characteristics: Distribution ranges from -2.86 to 1.27 with a mean of -0.26, indicating slight negative
asymmetry in seismic signal distributions. The high coefficient of variation (2.769) demonstrates substantial
variability, making skewness a potentially powerful discriminative feature for classification algorithms.
Kurtosis Analysis: Values span from -0.75 to 11.58 with a mean of 1.85, measuring signal distribution peakedness.
Higher kurtosis values indicate more impulsive, spiky signals characteristic of sudden seismic events, while lower
values suggest more uniform background noise patterns.
2.2 Feature Correlations and Relationships
The correlation analysis reveals a strong negative correlation (-0.587) between skewness and kurtosis,
indicating that more symmetric signals tend to exhibit higher peakedness. This relationship provides valuable
complementary information for machine learning algorithms, as the two parameters capture different aspects of
signal distribution characteristics while maintaining statistical independence.
3. Spectral Domain Feature Analysis
3.1 Welch Power Spectral Density Characteristics
The Welch PSD analysis, implemented with a 512-sample Hamming window, 50% overlap, and 1024 DFT
points, provides robust spectral energy estimates resistant to noise contamination . The analysis reveals
significant spectral characteristics consistent with seismic signal properties.
Spectral Energy Distribution: Welch mean power ranges from 54,341 to 1,977,830 with substantial variability (CV =
0.822), indicating diverse energy levels across different seismic events. The wide dynamic range suggests the
presence of both low-energy background activity and high-energy seismic events within the analyzed dataset.
Spectral Shape Parameters: The spectral centroid averages 5.51 Hz with a spread of 4.48 Hz, indicating that
seismic energy concentrates in the low-frequency range consistent with ground motion characteristics. This
frequency concentration aligns with established knowledge of seismic wave propagation through soil and rock layers.
3.2 FFT Frequency Bin Analysis
The FFT analysis using 1024 DFT points with 10 Hz frequency bins up to 100 Hz reveals the characteristic
frequency distribution pattern of seismic signals. The energy decreases exponentially with frequency, with the
dominant energy concentrated in the 0-10 Hz band (430,075 average magnitude).
Frequency Energy Distribution:
0-10 Hz: 430,075 (dominant band containing 60.8% of total energy)
10-20 Hz: 171,202 (24.2% of total energy)
20-30 Hz: 47,200 (6.7% of total energy)
Higher frequencies (30-100 Hz): Progressively decreasing energy content
This distribution pattern is characteristic of seismic signals where most energy concentrates below 20 Hz,
consistent with seismic wave propagation theory and ground motion characteristics. The dominant frequency analysis
shows a mean of 2.26 Hz, further confirming the low-frequency nature of the analyzed seismic events.
3.3 Spectral Feature Discriminative Power
The spectral features demonstrate high variability and discriminative potential for classification applications.
The coefficient of variation values (Welch mean: 0.822, Welch std: 0.821) indicate substantial inter-window
differences, making these features valuable for automated seismic event classification systems.
4. Time-Frequency Domain Analysis
4.1 Discrete Wavelet Transform Implementation
The DWT implementation using Daubechies-4 wavelets with 5-level decomposition provides multi-resolution
analysis capabilities essential for seismic signal characterization. This approach captures both transient and
sustained signal components across different frequency bands, enabling comprehensive time-frequency analysis.
4.2 Wavelet Energy Distribution Analysis
The energy distribution across wavelet decomposition levels reveals the multi-scale structure of seismic
signals:
Approximation Coefficients (L0): 64.1% of total energy, representing low-frequency trend components that
capture the fundamental ground motion characteristics.
Detail Coefficients Distribution:
L1 (highest frequency details): 22.0% - captures noise and high-frequency artifacts
L2 (mid-high frequency): 11.7% - contains seismic event characteristics
L3 (mid-frequency): 1.7% - represents structural ground motion patterns
L4-L5 (low-frequency details): <1% combined - background variations
This energy distribution is typical of seismic signals where low-frequency components dominate, but
significant mid-frequency content exists for event characterization. The concentration of energy in the
approximation and first detail levels (86.1% combined) demonstrates the effectiveness of wavelet
decomposition for seismic signal analysis.
4.3 Multi-Resolution Feature Significance
The wavelet features provide the most comprehensive signal characterization among the three domains. The relatively
low coefficient of variation for the approximation level (0.214) indicates stable low-frequency content, while
higher variation in detail levels (L1 CV: 0.411) suggests these components effectively capture transient seismic
events.
5. Data Quality and Statistical Validation
5.1 Data Integrity Assessment
The experimental dataset demonstrates high quality with no missing values and no infinite values, confirming proper
implementation of signal processing algorithms . The feature ranges span appropriate orders of magnitude, with
temporal features showing controlled variations and spectral features exhibiting the expected wide dynamic range
typical of seismic data.
5.2 Feature Stability and Reliability
The coefficient of variation analysis reveals varying degrees of feature stability:
High Stability Features (CV < 0.5):
Temporal standard deviation (0.401)
FFT 0-10 Hz bin (0.408)
Wavelet L0 energy ratio (0.214)
Moderate Stability Features (0.5 ≤ CV < 1.0):
Welch spectral parameters (0.821-0.822)
FFT 10-20 Hz bin (0.487)
High Variability Features (CV ≥ 1.0):
Temporal skewness (2.769)
Temporal kurtosis (1.070)
This variability distribution suggests that high-stability features provide consistent baseline characteristics,
while high-variability features offer strong discriminative power for classification tasks.
Code Implementaion
[Link]
welch(x=window, fs=sampling_rate, window='hamming', nperseg=512, noverLap=256, nfft=1024)
Estimate power spectral density using Welch's method.
Welch's method computes an estimate of the power spectral density by dividing the data into overlapping segments,
computing a modified periodogram for each segment and averaging the periodograms.
Parameters:
x : array_like
Time series of measurement values
fs : float, optional
Sampling frequency of the x time series. Defaults to 1.0.
window : str or tuple or array_like, optional
Desired window to use. If window is a string or tuple, it is passed to get_window to generate the window values, which are DFT-even
by default. If window is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann
window.
nperseg : int, optional
Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the
length of the window.
noverlap : int, optional
Number of points to overlap between segments. If None, noverlap = nperseg // 2 . Defaults to None.
nfft : int, optional
Length of the FFT used, if a zero padded FFT is desired. If None, the FFT length is nperseg. Defaults to None.
Returns:
f : ndarray
Array of sample frequencies.
Pxx : ndarray
Power spectral density or power spectrum of x .
Discrete Fourier Transform
Because the discrete Fourier transform separates its input into components that contribute at discrete frequencies, it has a
great number of applications in digital signal processing, e.g., for filtering, and in this context the discretized input to the transform is
customarily referred to as a signal, which exists in the time domain. The output is called a spectrum or transform and exists in the
frequency domain.
Implementation details
There are many ways to define the DFT, varying in the sign of the exponent, normalization, etc. In this implementation, the DFT is
defined as
{ }
n −1
mk
A k = ∑ ❑a m exp − 2 πi k=0 , … ,n − 1
m=0 n
The DFT is in general defined for complex inputs and outputs,
A single-frequency component at linear frequency f is represented by a complex exponential a m=exp {2 πifm Δ t },
where Δ t is the sampling interval.
If A=fft (a , n),
I. Then A [0] contains the zero-frequency term (the sum of the signal), which is always purely real for real inputs.
II. Then A [1 :n/2] contains the positive-frequency terms, in order of decreasingly negative frequency
III. A [n /2+1 :] contains the negative frequency terms, in order of decreasingly negative frequency.
IV. For an even number of input points, A [n /2 ] represents both positive and negative Nyquist frequency, and is also
purely real for real input.
V. For an odd number of input points, A [(n −1)/2] contains the largest positive frequency, while A [(n+1)/2]
contains the largest negative frequency.
VI. The routine [Link]( n ) returns an array giving the frequencies of corresponding elements in the output.
VII. The routine [Link](A) shifts transforms and their frequencies to put the zero-frequency components in the
middle
VIII. [Link](A) undoes that shift.
When the input a is a time-domain signal and A=fft (a) then:-
i. np ⋅|( A)| is its amplitude spectrum
ii. np ⋅|( A)|**2 is its power spectrum.
iii. The phase spectrum is obtained by np ⋅ a ngle(A ).
[Link]
[Link] (a, n=None)
This function computes the one-dimensional n -point discrete Fourier Transform (DFT) with the efficient Fast Fourier Transform
(FFT) algorithm.
Parameters:
a : array_like
Input array, can be complex.
n : int, optional
Length of the transformed axis of the output. If n is smaller than the length of the input, the input is cropped. If it is larger,
the input is padded with zeros. If n is not given, the length of the input along the axis specified by axis is used.
Returns:
out : complex ndarrayThe truncated or zero-padded input, transformed along the axis indicated by axis, or the last one if
axis is not specified.
fftfreq
Frequency bins for given FFT parameters.
Notes
FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, by using symmetries
in the calculated terms. The symmetry is highest when n is a power of 2 , and the transform is therefore most efficient for
these sizes.
[Link]
[Link] ( n , d=1.0 )
Return the Discrete Fourier Transform sample frequencies.
The returned float array f contains the frequency bin centers in cycles per unit of the sample spacing (with zero at the start).
For instance, if the sample spacing is in seconds, then the frequency unit is cycles/second.
Given a window length n and a sample spacing d :
f = [0,1,..., n/2-1, -n/2,..., -1] / (d*n) if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd
Parameters:
n : int
Window length.
d : scalar, optional
Sample spacing (inverse of the sampling rate). Defaults to 1.
.
Returns:
f : ndarray
Array of length n containing the sample frequencies.
Multilevel decomposition using wavedec
[Link](data, wavelet, Level=None)
Multilevel 1D Discrete Wavelet Transform of data.
Parameters:
data: array_like
Input data
wavelet : Wavelet object or name string
Wavelet to use
level : int, optional
Decomposition level (must be >= 0). If level is None (default) then it will be calculated using the dwt_max_level function.
Returns:
[cA_n, cD_n, cD_n-1, ..., cD2, cD1] : list
Ordered list of coefficients arrays where n denotes the level of decomposition. The first element ( CA n ) of the result is
approximation coefficients array and the following elements (CD_n - CD_1 ) are details coefficients arrays.
Maximum decomposition level dwt_max_level, dwtn_max_level
pywt.dwt_max_level(data_len, filter_len)
Compute the maximum useful level of decomposition.
Parameters:
data_len : int
Input data length.
filter_len : int, str or Wavelet
The wavelet filter length. Alternatively, the name of a discrete wavelet or a Wavelet object can be specified.
Returns:
max_level : int
Maximum level.
The rational for the choice of levels is the maximum level where at least one coefficient in the output is uncorrupted by
edge effects caused by signal extension. Put another way, decomposition stops when the signal becomes shorter than
the FIR filter length for a given wavelet. This corresponds to:
( data_len
max_level = ⌊ log 2
filter_len −1 )
⌋