Skip to content

cds2.wavelets

Haar wavelet transforms and wavelet denoising.

Pure NumPy Haar wavelet transforms: decomposition, reconstruction, denoising.

DWTResult dataclass

Multi-level Haar decomposition of a finite series.

Source code in src/cds2/wavelets.py 
69
70
71
72
73
74
75
@dataclass(frozen=True)
class DWTResult:
    """Multi-level Haar decomposition of a finite series."""

    approximations: list[FloatArray]
    details: list[FloatArray]
    original_length: int

DenoiseResult dataclass

Outcome of wavelet shrinkage denoising.

Source code in src/cds2/wavelets.py 
109
110
111
112
113
114
115
@dataclass(frozen=True)
class DenoiseResult:
    """Outcome of wavelet shrinkage denoising."""

    data: FloatArray
    thresholded_counts: int
    levels: int

haar_dwt(signal, normalize=True)

Single-level Haar decomposition into pairwise averages and differences.

Source code in src/cds2/wavelets.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
def haar_dwt(
    signal: Sequence[float] | FloatArray, normalize: bool = True
) -> tuple[FloatArray, FloatArray]:
    """Single-level Haar decomposition into pairwise averages and differences."""
    series = _as_signal(signal)
    if series.size % 2 != 0:
        msg = "signal length must be even"
        raise ValueError(msg)
    pairs = series.reshape(-1, 2)
    scale = np.sqrt(2.0) if normalize else 2.0
    approx = np.asarray((pairs[:, 0] + pairs[:, 1]) / scale)
    detail = np.asarray((pairs[:, 0] - pairs[:, 1]) / scale)
    return approx, detail

haar_idwt(approx, detail, normalize=True)

Single-level Haar inverse transform restoring the paired samples.

Source code in src/cds2/wavelets.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def haar_idwt(
    approx: Sequence[float] | FloatArray,
    detail: Sequence[float] | FloatArray,
    normalize: bool = True,
) -> FloatArray:
    """Single-level Haar inverse transform restoring the paired samples."""
    averages = np.asarray(approx, dtype=float)
    differences = np.asarray(detail, dtype=float)
    if averages.ndim != 1 or differences.ndim != 1 or averages.size < 1 or differences.size < 1:
        msg = "approximation and detail must be non-empty 1-D sequences"
        raise ValueError(msg)
    if averages.shape != differences.shape:
        msg = "approximation and detail must have the same length"
        raise ValueError(msg)
    scale = np.sqrt(2.0) if normalize else 1.0
    restored = np.empty(averages.size * 2, dtype=float)
    restored[0::2] = (averages + differences) / scale
    restored[1::2] = (averages - differences) / scale
    return restored

dwt_levels(signal, levels, normalize=True)

Multi-level Haar decomposition applying haar_dwt to each approximation.

Source code in src/cds2/wavelets.py 
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def dwt_levels(
    signal: Sequence[float] | FloatArray, levels: int, normalize: bool = True
) -> DWTResult:
    """Multi-level Haar decomposition applying ``haar_dwt`` to each approximation."""
    if levels < 1:
        msg = "levels must be at least 1"
        raise ValueError(msg)
    series = _as_signal(signal)
    if series.size < 2**levels:
        msg = "signal too short"
        raise ValueError(msg)
    approximations: list[FloatArray] = []
    details: list[FloatArray] = []
    current = series
    for _ in range(levels):
        current, detail = haar_dwt(current, normalize=normalize)
        approximations.append(current)
        details.append(detail)
    return DWTResult(
        approximations=approximations, details=details, original_length=int(series.size)
    )

idwt_levels(result, normalize=True)

Reconstruct the original signal from a :class:DWTResult.

Source code in src/cds2/wavelets.py 
101
102
103
104
105
106
def idwt_levels(result: DWTResult, normalize: bool = True) -> FloatArray:
    """Reconstruct the original signal from a :class:`DWTResult`."""
    current = result.approximations[-1]
    for detail in reversed(result.details):
        current = haar_idwt(current, detail, normalize=normalize)
    return current

wavelet_denoise(signal, threshold_factor=3.0, levels=2, normalize=True)

Zero detail coefficients below threshold_factor standard deviations per level.

Source code in src/cds2/wavelets.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def wavelet_denoise(
    signal: Sequence[float] | FloatArray,
    threshold_factor: float = 3.0,
    levels: int = 2,
    normalize: bool = True,
) -> DenoiseResult:
    """Zero detail coefficients below ``threshold_factor`` standard deviations per level."""
    decomposition = dwt_levels(signal, levels, normalize=normalize)
    filtered: list[FloatArray] = []
    zeroed = 0
    for detail in decomposition.details:
        sigma = float(np.median(np.abs(detail)) / 0.6745)
        cutoff = threshold_factor * sigma
        keep = np.abs(detail) >= cutoff
        zeroed += int((~keep).sum())
        filtered.append(np.asarray(np.where(keep, detail, 0.0)))
    reconstructed = DWTResult(
        approximations=decomposition.approximations,
        details=filtered,
        original_length=decomposition.original_length,
    )
    return DenoiseResult(
        data=idwt_levels(reconstructed, normalize=normalize),
        thresholded_counts=zeroed,
        levels=levels,
    )