base/WAVELET.md
The smile.wavelet package provides a fast, in-place Discrete Wavelet Transform (DWT) and its inverse for 1-D and 2-D signals, together with wavelet shrinkage for automatic signal denoising.
Key properties of the implementation:
| Property | Value |
|---|---|
| Transform type | Orthogonal (energy-preserving) |
| Algorithm | Mallat pyramid algorithm |
| In-place | Yes — the signal array is modified in-place |
| Complexity | O(N) per level, O(N) total |
| Signal length | Power of 2 |
import smile.wavelet.*;
// 1. Pick a wavelet
Wavelet wavelet = new DaubechiesWavelet(8);
// 2. Forward transform (in-place)
double[] signal = { 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0, 0.0 };
wavelet.transform(signal);
// signal[] now holds the wavelet coefficients
// 3. Inverse transform (in-place) — recover the original signal
wavelet.inverse(signal);
// 4. Denoise a noisy signal in one line
double[] noisy = { /* your data */ };
WaveletShrinkage.denoise(noisy, new HaarWavelet());
All wavelets extend the base Wavelet class:
Wavelet
├── HaarWavelet
├── D4Wavelet
├── DaubechiesWavelet (n = 4, 6, 8, 10, 12, 14, 16, 18, 20)
├── CoifletWavelet (n = 6, 12, 18, 24, 30)
├── SymletWavelet (n = 8, 10, 12, 14, 16, 18, 20)
└── BestLocalizedWavelet (n = 14, 18, 20)
The simplest possible wavelet — a step function. Equivalent to D2 (Daubechies with 2 coefficients). Computes pair-wise averages and differences, producing a fast transform with the lowest computational cost.
Wavelet w = new HaarWavelet();
Best for: piecewise-constant signals, quick prototyping, signal with abrupt discontinuities.
Daubechies D4 (4-coefficient) wavelet using a different centering convention from DaubechiesWavelet(4). The forward and inverse steps are computed explicitly (no generic filter loop) for maximum speed.
Wavelet w = new D4Wavelet();
Best for: signals with linear trends. One extra vanishing moment compared to Haar.
The standard Daubechies family with maximal vanishing moments for a given support width n. The number of vanishing moments equals n/2.
Wavelet w4 = new DaubechiesWavelet(4); // 2 vanishing moments
Wavelet w6 = new DaubechiesWavelet(6); // 3 vanishing moments
Wavelet w8 = new DaubechiesWavelet(8); // 4 vanishing moments
// ... up to n=20 (10 vanishing moments)
Supported orders: 4, 6, 8, 10, 12, 14, 16, 18, 20.
Best for: general-purpose use; more vanishing moments → better polynomial approximation at the cost of a wider support.
Coiflets are designed so that both the wavelet function and the scaling function have vanishing moments. This makes the scaling coefficients better approximations of signal samples, which is useful in numerical analysis.
Wavelet c6 = new CoifletWavelet(6); // 2 vanishing moments each
Wavelet c12 = new CoifletWavelet(12); // 4 vanishing moments each
// ... up to n=30
Supported orders: 6, 12, 18, 24, 30 (multiples of 6).
Best for: signals where the coarsest approximation coefficients should closely match the original sample values (e.g. numerical integration, quadrature).
Symlets (S8–S20) are near-symmetric modifications of the Daubechies family. They have the same number of vanishing moments for a given support, but the filter energy is more symmetrically distributed, reducing phase distortion.
Wavelet s8 = new SymletWavelet(8);
Wavelet s12 = new SymletWavelet(12);
// ... up to n=20
Supported orders: 8, 10, 12, 14, 16, 18, 20 (even numbers only).
Best for: audio processing and signals where linear phase response matters.
Best Localized (BL) wavelets minimize time-frequency uncertainty while maintaining compact support. They are better localized in time than standard Daubechies wavelets of the same order.
Wavelet bl14 = new BestLocalizedWavelet(14);
Wavelet bl18 = new BestLocalizedWavelet(18);
Wavelet bl20 = new BestLocalizedWavelet(20);
Supported orders: 14, 18, 20.
Best for: transient detection in signals where both time and frequency localization matter.
Any orthogonal filter bank coefficients can be wrapped in the base Wavelet class directly:
// Example: Daubechies D4 coefficients
double[] h = {
0.4829629131445341, 0.8365163037378079,
0.2241438680420134, -0.1294095225512604
};
Wavelet myWavelet = new Wavelet(h);
The constructor automatically derives the high-pass (wavelet) filter g from the low-pass (scaling) filter h via the quadrature mirror filter (QMF) relationship:
g[n] = (-1)^n * h[N-1-n]
wavelet.transform(double[] a) applies the full multi-level DWT in-place. Starting from the full signal length N, it repeatedly applies the one-level forward filter, halving the active length at each step, until the active length drops below ncof (the number of filter coefficients).
double[] a = { 1, 2, 3, 4, 5, 6, 7, 8 }; // length must be a power of 2
new DaubechiesWavelet(4).transform(a);
// a[] now contains the multi-resolution wavelet coefficients
The resulting array is ordered coarsest-to-finest:
[ s1 | d1 | d2 .... | d3 ............ ]
^ ^ ^--- level 2 detail ^--- level 3 detail (finest)
| |
| +--- level 1 detail (coarsest)
+--- level 0 scaling (global mean)
wavelet.inverse(double[] a) reconstructs the original signal from its wavelet coefficients in-place. It is the exact numerical inverse of transform():
double[] original = { 1, 2, 3, 4, 5, 6, 7, 8 };
double[] a = original.clone();
Wavelet w = new DaubechiesWavelet(8);
w.transform(a);
w.inverse(a);
// a[] ≈ original[] (to within floating-point rounding, ≈ 1e-14)
| Requirement | Detail |
|---|---|
| Length | Must be a power of 2 (2, 4, 8, 16, …) |
| Minimum length | Must be ≥ the number of filter coefficients ncof |
| In-place | The input array is modified; clone first if the original is needed |
Both transform() and inverse() throw IllegalArgumentException if the input length is not a power of 2, or is shorter than ncof.
// Zero-pad to the next power of 2 for non-power-of-2 signals:
int n = signal.length;
int padded = Integer.highestOneBit(n - 1) << 1; // next power of 2
double[] a = Arrays.copyOf(signal, padded); // right-pad with zeros
wavelet.transform(a);
After transform(), index 0 holds the global scaling coefficient (proportional to the signal mean). This coefficient is never thresholded by WaveletShrinkage.
WaveletShrinkage is a utility interface with only static methods. No instantiation is needed.
The threshold is computed automatically using Donoho & Johnstone's universal threshold:
λ = σ · √(2 · log(N))
where:
N is the signal length,σ is estimated from the finest-scale detail coefficients (the upper half of the transformed array) using the scaled Median Absolute Deviation (MAD):σ̂ = MAD(finest detail coefficients) / 0.6745
This estimator is robust to signal content because the finest-scale detail coefficients are dominated by noise (not signal) for smooth signals.
| Mode | Formula | Effect |
|---|---|---|
| Hard (default) | `c → c if | c |
| Soft | `c → sign(c)·max( | c |
Hard thresholding tends to produce better MSE for sparse signals; soft thresholding gives smoother reconstructions and is preferred when the signal is expected to be continuous.
import smile.wavelet.*;
double[] signal = /* ... your noisy 1-D signal, length must be a power of 2 ... */;
// Hard thresholding (default)
WaveletShrinkage.denoise(signal, new DaubechiesWavelet(8));
// Explicit hard thresholding
WaveletShrinkage.denoise(signal, new DaubechiesWavelet(8), false);
// Soft thresholding
WaveletShrinkage.denoise(signal, new DaubechiesWavelet(8), true);
Note: denoise() modifies the array in-place. Clone it first if the original is needed.
The global scaling coefficient signal[0] is always preserved — it encodes the mean of the signal and is never zeroed by thresholding.
The 2-D extension applies the 1-D denoising independently to each row, then to each column (separable 2-D wavelet transform):
double[][] image = /* rows × cols, both dimensions must be powers of 2 */;
// Hard thresholding (default)
WaveletShrinkage.denoise2D(image, new HaarWavelet());
// Soft thresholding
WaveletShrinkage.denoise2D(image, new HaarWavelet(), true);
| Use case | Recommended wavelet |
|---|---|
| General denoising / prototyping | DaubechiesWavelet(8) |
| Piecewise-constant or discontinuous signals | HaarWavelet |
| Signals with smooth, polynomial-like trends | DaubechiesWavelet(12) or DaubechiesWavelet(20) |
| Audio / minimal phase distortion | SymletWavelet(8) – SymletWavelet(20) |
| Numerical integration / quadrature | CoifletWavelet(12) |
| Transient / event detection | BestLocalizedWavelet(14) |
| Speed-critical applications | HaarWavelet or D4Wavelet |
As a rule of thumb: more filter coefficients → more vanishing moments → better suppression of smooth signal content in detail coefficients → more accurate denoising for smooth signals, at the cost of a wider support and slightly higher boundary effects.
The DWT requires that the signal length be a power of 2. For signals of arbitrary length, zero-pad on the right:
import java.util.Arrays;
double[] rawSignal = /* length N, not necessarily a power of 2 */;
// Find the next power of 2
int n = rawSignal.length;
int paddedLength = n == 1 ? 1 : Integer.highestOneBit(n - 1) << 1;
double[] padded = Arrays.copyOf(rawSignal, paddedLength); // pads with 0.0
Wavelet w = new DaubechiesWavelet(8);
WaveletShrinkage.denoise(padded, w);
// Trim back to original length
double[] denoised = Arrays.copyOf(padded, n);
Tip: Symmetric (periodic) padding often gives better results than zero padding at signal boundaries, especially when the signal does not start and end near zero. For simple use cases, zero padding is sufficient.
The DWT is computed via the Mallat pyramid (also called sub-band coding). At each level, the signal of length n is decomposed by:
h (downsampled by 2) → approximation coefficients of length n/2.g (downsampled by 2) → detail coefficients of length n/2.Both convolutions use circular (periodic) boundary conditions, so the output has exactly the same length as the input.
The high-pass filter is derived from the low-pass filter via the quadrature mirror relationship:
g[k] = (-1)^k · h[N−1−k]
This guarantees that the transform is orthogonal (energy-preserving, Parseval's theorem holds).
For any orthogonal wavelet, the DWT preserves signal energy (Parseval's theorem):
Σ a[i]² = Σ w[i]²
where a is the original signal and w are the wavelet coefficients.
For a length-N signal with i.i.d. Gaussian noise of standard deviation σ, the universal threshold
λ = σ · √(2 log N)
guarantees (with high probability as N → ∞) that all pure-noise coefficients are zeroed while at least one signal coefficient survives, minimising the minimax risk over a broad class of signals (Donoho & Johnstone, 1994).
import smile.wavelet.*;
public class RoundTripExample {
public static void main(String[] args) {
double[] signal = {1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 0};
double[] original = signal.clone();
Wavelet w = new DaubechiesWavelet(8);
w.transform(signal);
System.out.println("Wavelet coefficients:");
for (int i = 0; i < signal.length; i++) {
System.out.printf(" w[%2d] = %10.6f%n", i, signal[i]);
}
w.inverse(signal);
System.out.println("\nReconstruction error:");
double maxErr = 0;
for (int i = 0; i < signal.length; i++) {
maxErr = Math.max(maxErr, Math.abs(signal[i] - original[i]));
}
System.out.printf(" max |error| = %.2e%n", maxErr); // expect ~1e-15
}
}
import smile.wavelet.*;
import java.util.Arrays;
public class DenoiseExample {
public static void main(String[] args) {
int n = 128;
double[] clean = new double[n];
double[] noisy = new double[n];
// Build a clean signal + broadband noise
java.util.Random rng = new java.util.Random(42);
for (int i = 0; i < n; i++) {
clean[i] = Math.sin(2 * Math.PI * i / 32.0);
noisy[i] = clean[i] + 0.4 * rng.nextGaussian();
}
// Denoise with Daubechies D8, hard thresholding
double[] denoised = noisy.clone();
WaveletShrinkage.denoise(denoised, new DaubechiesWavelet(8), false);
// Report MSE improvement
System.out.printf("Noisy MSE = %.6f%n", mse(clean, noisy));
System.out.printf("Denoised MSE = %.6f%n", mse(clean, denoised));
}
static double mse(double[] a, double[] b) {
double sum = 0;
for (int i = 0; i < a.length; i++) { double d = a[i]-b[i]; sum += d*d; }
return sum / a.length;
}
}
import smile.wavelet.*;
public class SoftVsHardExample {
public static void main(String[] args) {
double[] signal = buildNoisySignal(256);
double[] hard = signal.clone();
WaveletShrinkage.denoise(hard, new SymletWavelet(8), false);
double[] soft = signal.clone();
WaveletShrinkage.denoise(soft, new SymletWavelet(8), true);
System.out.println("Hard-thresholded signal (first 8): "
+ java.util.Arrays.toString(java.util.Arrays.copyOf(hard, 8)));
System.out.println("Soft-thresholded signal (first 8): "
+ java.util.Arrays.toString(java.util.Arrays.copyOf(soft, 8)));
}
static double[] buildNoisySignal(int n) {
double[] t = new double[n];
java.util.Random rng = new java.util.Random(0);
for (int i = 0; i < n; i++)
t[i] = Math.cos(2 * Math.PI * i / 64.0) + 0.3 * rng.nextGaussian();
return t;
}
}
import smile.wavelet.*;
public class Denoise2DExample {
public static void main(String[] args) {
int n = 64; // must be a power of 2
double[][] image = new double[n][n];
// Synthetic image: smooth 2-D cosine + noise
java.util.Random rng = new java.util.Random(1);
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
image[i][j] = Math.cos(2 * Math.PI * i / 16.0)
* Math.cos(2 * Math.PI * j / 16.0)
+ 0.3 * rng.nextGaussian();
WaveletShrinkage.denoise2D(image, new DaubechiesWavelet(8), false);
System.out.println("Centre pixel after denoising: " + image[n/2][n/2]);
}
}
import smile.wavelet.*;
public class CustomWaveletExample {
// Daubechies D6 coefficients
private static final double[] D6 = {
0.3326705529500825, 0.8068915093110924, 0.4598775021184914,
-0.1350110200102546, -0.0854412738820267, 0.0352262918857095
};
public static void main(String[] args) {
Wavelet w = new Wavelet(D6);
double[] a = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
15,14,13,12,11,10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0};
double[] original = a.clone();
w.transform(a);
w.inverse(a);
double maxErr = 0;
for (int i = 0; i < a.length; i++)
maxErr = Math.max(maxErr, Math.abs(a[i] - original[i]));
System.out.printf("Custom D6 roundtrip error: %.2e%n", maxErr);
}
}
SMILE — © 2010-2026 Haifeng Li. GNU GPL licensed.