Port origin: Algorithms ported from Madagascar (© 2004–2008 University of Texas at Austin) into the Seismic Unix (SU) framework.
Primary application: Ground-roll (surface-wave) attenuation in 2-D seismic reflection data.
suseislet is a C-language signal-processing toolkit that extends Seismic Unix with three new command-line programs:
| Program | Role |
|---|---|
sudip |
Estimates the local dip field (instantaneous slope of seismic events) using plane-wave destruction |
suseislet |
Applies the forward or inverse seislet transform guided by the estimated dip field |
suthreshold |
Thresholds seislet-domain coefficients to separate signal from noise |
Together they implement a complete transform-threshold-reconstruct workflow for coherent noise attenuation — analogous to wavelet shrinkage, but adapted to the geometry of seismic data.
In seismic exploration, ground roll (Rayleigh-wave surface noise) contaminates reflection records. It is low-frequency, high-amplitude, and coherent — making it hard to remove with simple frequency filters without damaging the target reflections. The seislet transform exploits the structural geometry of seismic events (their local dip) to separate ground roll from reflections in a sparse coefficient domain.
- Dip-adaptive transform — the seislet basis functions follow the local slope of seismic events, producing highly sparse representations of coherent wavefields
- Three lifting-scheme wavelet types:
h— Haar (simplest, O(N) operations, low quality)l— Linear (better quality, uses two neighboring traces for prediction)b— Biorthogonal CDF 9/7 (highest quality, four lifting steps + scaling)
- Plane-wave prediction engine — trace-to-trace prediction is performed by solving a regularised linear system using the plane-wave destruction (PWD) all-pass filter
- Three thresholding modes — mute-smooth (ground-roll removal), keep-smooth (ground-roll extraction), and IQR shrinkage (statistical soft thresholding)
- Power-of-2 panel processing — the lifting scheme operates on panels of
n2 = 32, 64, 128, …traces, consistent with classical wavelet multiresolution analysis - Full invertibility — the forward and inverse transforms are exact inverses when no thresholding is applied (perfect reconstruction)
- Seismic Unix native — reads and writes standard SU (SEGY-like) binary trace files; integrates seamlessly into existing SU pipelines
- Ported from Madagascar — algorithm code is mathematically identical to the peer-reviewed Madagascar reference implementation; only the memory/IO layer was adapted
Classical wavelets decompose a signal by repeatedly splitting it into smooth (scaling) and detail (wavelet) sub-bands via convolution and downsampling. The lifting scheme (Sweldens, 1996) achieves the same result using only in-place additions and multiplications — no convolutions — through two alternating steps:
Given: even-indexed samples s, odd-indexed samples d
PREDICT step: d ← d − P(s) (detail = residual of prediction)
UPDATE step: s ← s + U(d) (smooth = corrected average)
The seislet transform replaces the classical shift-based prediction operator P with a plane-wave prediction operator that shifts traces along the local dip direction. This makes the detail coefficients near zero wherever the data follow the estimated dip — producing an extremely sparse representation.
A plane wave in a seismic panel satisfies:
u(t, x + Δx) ≈ u(t − p·Δx, x)
where p (samples per trace, or dt/dx · tan θ) is the local dip (slope). The plane-wave destruction filter (Fomel, 2002) is defined such that:
D(p) · u ≈ 0 when u is a plane wave with slope p
D(p) is implemented as a banded all-pass filter with coefficients computed by a binomial approximation (Wilson, 1992). For a filter of order nw, the coefficient vector a[k] for slope p satisfies:
passfilter(p):
a[k] = b[k] · ∏_{j<n−k} (n − j − p) · ∏_{j≥n−k} (p + j + 1)
where b[k] are pre-computed binomial base weights. This all-pass filter achieves a fractional-sample shift of exactly p in the z-domain, enabling sub-sample accuracy in dip tracking.
sudip estimates p(t, x) iteratively. At each outer iteration iter:
- Apply the current plane-wave filter:
u2 = D(p) · u(destruction residual) - Apply its slope-derivative:
u1 = ∂D/∂p · u(sensitivity) - Normalise both by the RMS of
u1 - Compute a smooth dip update
dpby solving a regularised least-squares division:
dp = argmin_q ‖u1·q − u2‖² + ε²·L²q
where L² is a 2D triangle smoother (controlled by rect1, rect2). The regularization suppresses spurious oscillations in the dip field.
- Update:
p ← p + dp(additive) orp ← p · exp(dp/p)(multiplicative, preserves sign)
For a panel of n2 traces and n1 time samples, the transform iterates over dyadic scale levels j = 1, 2, 4, …, n2/2. At each level, even-indexed trace groups are the smooth component and odd-indexed groups are the detail.
PREDICT: d[i+j] ← d[i+j] − predict_forward(s[i])
UPDATE: s[i] ← s[i] + predict_backward(d[i+j]) / 2
predict_forward(s[i]) shifts trace s[i] forward by j traces along the estimated dip field using j sequential calls to predict_step.
Uses two neighbors for prediction (one on each side), giving a second-order accurate prediction:
PREDICT: d[i+j] ← d[i+j] − (predict_forward(s[i]) + predict_backward(s[i+2j])) / 2
UPDATE: s[i] ← s[i] + predict_backward(d[i+j]) / 2 (boundary)
s[i] ← s[i] + (predict_forward(d[i+j]) + predict_backward(d[i−j])) / 4 (interior)
This is the Cohen–Daubechies–Feauveau 9/7 wavelet (used in JPEG 2000), adapted for seismic dip. It has four lifting steps and a normalization step with fixed constants:
| Step | Constant α |
|---|---|
| Predict 1 (α₁) | +1.586134342 |
| Update 1 (α₂) | −0.05298011854 |
| Predict 2 (α₃) | −0.8829110762 |
| Update 2 (α₄) | +0.4435068522 |
| Scale (K) | 1.230174105 |
Each predict/update step replaces the fixed shift of a standard CDF 9/7 wavelet with the plane-wave prediction operator, making the filter geometry-adaptive.
After the forward transform, the output panel is organized as:
Trace index 0: SMOOTH component (lowest frequency / coarsest scale)
Trace index 1: DETAIL at scale j=n2/2
Trace index 2: DETAIL at scale j=n2/4
Trace index 3: DETAIL at scale j=n2/4
...
Trace index n2−1: DETAIL at scale j=1 (finest scale)
Ground roll, being a coherent low-wavenumber event, concentrates energy into trace 0 (the smooth component). Reflections, being nearly orthogonal to the dip basis, distribute into the detail traces.
For statistical thresholding, the interquartile range (IQR) method is used:
Q1 = 25th percentile of trace amplitudes
Q3 = 75th percentile of trace amplitudes
IQR = Q3 − Q1
dmin = Q1 − scale × IQR
dmax = Q3 + scale × IQR
threshold(x):
if x < dmin: x ← x − dmin
if x > dmax: x ← x − dmax
else: x ← 0
This is a soft-thresholding variant: large outliers are shrunk toward the fence, and typical amplitudes are zeroed. The scale parameter (default 1.5) controls how aggressively the thresholding operates.
# Install dependency (Ubuntu/Debian)
sudo apt install libtirpc-dev
# Edit Makefile: set CWPROOT to your SeisUnix installation path
# e.g., CWPROOT = /usr/local/cwp
make # compile all three programs
make install # copy to $CWPROOT/bin (optional)# Step 1: Estimate local dip field from the raw shot gather
sudip n2=64 rect1=5 rect2=5 niter=5 < shot.su > dip.su
# Step 2: Forward seislet transform (data → seislet domain)
suseislet n2=64 dip=dip.su type=b < shot.su > seislet.su
# Step 3: Threshold — zero the smooth component (ground roll)
suthreshold mode=0 < seislet.su > thresh.su
# Step 4: Inverse seislet transform (seislet domain → clean data)
suseislet n2=64 dip=dip.su inv=1 type=b < thresh.su > clean.su# Keep only the smooth component (ground-roll model)
suthreshold mode=1 < seislet.su > gr_seislet.su
# Reconstruct the ground-roll model in data space
suseislet n2=64 dip=dip.su inv=1 < gr_seislet.su > groundroll.su
# Subtract ground-roll model from original data
suop2 shot.su groundroll.su op=diff > clean.su# Forward then inverse should return the original data
suseislet n2=64 dip=dip.su < shot.su > seislet.su
suseislet n2=64 dip=dip.su inv=1 < seislet.su > reconstructed.su
suop2 shot.su reconstructed.su op=diff > residual.su
# residual.su should be near-zero (machine precision)# Display seislet coefficients (trace 0 = smooth, rest = details)
suximage < seislet.su title="Seislet Domain" &# Haar wavelet (fastest, lowest quality)
suseislet n2=64 dip=dip.su type=h < shot.su > s_haar.su
# Linear wavelet (moderate quality)
suseislet n2=64 dip=dip.su type=l < shot.su > s_linear.su
# Biorthogonal CDF 9/7 (default, best quality)
suseislet n2=64 dip=dip.su type=b < shot.su > s_bio.su# Zero the two coarsest-scale components (trace 0 and 1)
suthreshold mode=0 nsmooth=2 < seislet.su > thresh.suThe core transform library. All state is module-level (not thread-safe).
Initializes the transform for a panel of n1 samples × n2 traces.
| Parameter | Type | Description |
|---|---|---|
n1 |
int |
Samples per trace |
n2 |
int |
Number of traces (must be power of 2) |
inv1 |
bool |
Always pass true (selects Predict-then-Update lifting order for correct reconstruction) |
unit1 |
bool |
If true, apply orthonormal scaling weights w[i] = sqrt(level) |
eps |
float |
Regularization for plane-wave prediction (typical: 0.01) |
order |
int |
All-pass filter order: 1 = first-order, 2 = second-order |
type |
char |
Wavelet type: 'h' Haar, 'l' linear, 'b' biorthogonal |
Registers the dip field float **dip[n2][n1] with the transform engine. Must be called after seislet_init and before seislet_lop.
Applies the transform as a linear operator.
adj |
Direction |
|---|---|
true |
Forward (analysis): reads y (data), writes x (seislet coefficients) |
false |
Inverse (synthesis): reads x (coefficients), writes y (data) |
add=false zeroes the output before accumulation (standard usage).
Frees all internally allocated memory and calls predict_close().
Alternative operator forms that internally flip the inv flag based on adj. Used in iterative inversion contexts.
Implements the core predict-and-solve engine used by the lifting steps.
| Parameter | Description |
|---|---|
nx |
Samples per trace (n1) |
ny |
Number of traces (n2) |
e |
Regularization parameter squared (eps²) |
nw |
All-pass filter half-order (1 or 2) |
k |
Starting trace index for predict_lop |
two |
If true, initialise a second PWD filter (w2) |
The fundamental operation: given a dip field pp[n1], solves the regularised system (R + AᵀA) · trace = rhs where A is the plane-wave destruction matrix. This computes the prediction of one trace from its neighbor along the dip direction.
forw=true→ shift trace forward (left-to-right across panel)forw=false→ shift trace backward (right-to-left)
Applies the prediction as a full-panel linear operator (used in iterative inversion).
Builds the banded matrix representation of the plane-wave destruction operator.
Allocates a pwd struct containing filter coefficient arrays a[na][n1] where na = 2*nw + 1. Calls apfilt_init(nw) internally.
Fills the diagonal (diag[n1]) and off-diagonal bands (offd[nb][n1]) of the banded system matrix by computing AᵀA contributions from the all-pass filter coefficients at each sample position.
Multiplies an input trace by the PWD banded matrix. Used in predict_step after pwd_define to apply the operator.
The lowest-level numerical kernel. No dependencies on other project modules.
Precomputes the binomial base coefficients b[k] for a filter of order 2*nw. Uses double precision internally for numerical accuracy.
Computes the 2*nw + 1 all-pass filter coefficients for fractional shift p (the local dip in samples/trace). These coefficients define the plane wave operator at a single sample location.
Computes the derivative of the filter coefficients with respect to slope p. Used in dip2.c (step 2 of dip estimation) to form the sensitivity matrix.
| Parameter | Description |
|---|---|
niter |
Inner CG iterations for the smooth division solver |
nx, ny |
Panel dimensions (n1, n2) |
fx, fy |
Smoothing radii for the triangle smoother (rect1, rect2) |
sign |
If true, use multiplicative dip update |
gauss |
Gaussian smoother flag (not implemented; pass false) |
verb |
Print convergence info each iteration |
Runs niter outer iterations of dip update. Modifies p[n2][n1] in-place. mask can be NULL or a boolean array to fix dip values at certain samples.
Allocates workspace of size n1 and stores the IQR fence multiplier scale.
Computes Q1 and Q3 from the data, then sets dmin = Q1 − scale·IQR and dmax = Q3 + scale·IQR.
Applies the soft threshold in-place: values outside [dmin, dmax] are shrunk toward the fence; values inside are zeroed.
Provides two portable replacements for RSF macros/functions:
SF_NINT(x)— rounds a float to the nearest integer (used inallp2.cfor aliasing-stretch index)adjnull(adj, add, nx, ny, x, y)— zeros the output array before accumulation; replacessf_adjnullfrom Madagascar's RSF library
suseislet (main)
└── seislet.c
└── predict.c
├── pwd.c
│ └── apfilt.c ← no further deps
└── suband.c ← banded matrix solver
sudip (main)
└── dip2.c
├── allp2.c
│ └── apfilt.c
└── div2.c
└── sudivutil.c ← smoothing + CG solver
suthreshold (main)
└── threshold.c ← no further deps
Important
n2 must be a power of 2. The lifting scheme requires dyadic subsampling. Use n2 = 32, 64, 128, 256, …. The code validates this at startup and aborts with an error if violated.
Warning
The inv1 flag in seislet_init should always be true. The forward (analysis) transform internally uses the Predict-then-Update (PU) lifting order when inv1=true. The inverse (synthesis) transform always undoes the PU scheme. Passing inv1=false for the forward pass creates a lifting-order mismatch and produces incorrect reconstruction. This is a known quirk of the ported code and is documented in suseislet.c.
Note
Dip units: The dip field p is in samples per trace (i.e., dt/dx · tan θ). It is a dimensionless slope in the discrete (sample, trace) coordinate system. Typical values for seismic reflection data lie in [-5, 5]; ground roll can have steeper dips.
Tip
Start with type=b (biorthogonal). For most real datasets, the CDF 9/7 biorthogonal wavelet gives the best separation between ground roll and reflections. Switch to type=h (Haar) only for diagnostic purposes or when computational speed is critical.
Note
Zero-padding: If the input has fewer than n2 traces, both sudip and suseislet zero-pad the panel and print a warning. The output always contains exactly n2 traces.
- Fomel, S. (2002). Applications of plane-wave destruction filters. Geophysics, 67(6), 1946–1960.
- Fomel, S. (2006). Towards the seislet transform. SEG Technical Program Expanded Abstracts, 2847–2851.
- Sweldens, W. (1996). The lifting scheme: A custom-design construction of biorthogonal wavelets. Applied and Computational Harmonic Analysis, 3(2), 186–200.
- Daubechies, I. & Sweldens, W. (1998). Factoring wavelet transforms into lifting steps. Journal of Fourier Analysis and Applications, 4(3), 247–269.
- Wilson, R. (1992). Minimax-optimal all-pass fractional delay filters. IEEE Transactions on Signal Processing.
- Madagascar software package — reproducibility.org. Original
sfseislet,sfdip, and related programs that this project ports to Seismic Unix.