API reference

dwt(x, b, l; convention = :aligned)

Return a transformed copy of x. The copying form of dwt!.

dwt!(x, b, l; convention = :aligned)
dwt!(x, w, b, l; convention = :aligned)

Transform x in place with basis b, and return x.

l sets the number of levels. An Int applies the same number along every axis; an NTuple{N,Int} or Vector{Int} sets the level along each axis separately. A dimension of length n supports l levels of perfect reconstruction when 2^l divides n, and is transformed as far as it divides beyond that.

The four-argument form takes w, a scratch array with the same size and element type as x, and avoids the allocation that the three-argument form makes.

convention selects the phase of the detail coefficients. :aligned (the default) applies the scaling and wavelet filters to the same input window. :wavelets reproduces the coefficient layout of Wavelets.jl. The scaling coefficients are identical in both conventions.

x may be any strided AbstractArray, including a subband view from rtree_view, or a GPU array.

Examples

using NDWaveletTransforms

x = rand(128, 128)
dwt!(copy(x), WT_D4, 2)         # two levels along both axes
dwt!(copy(x), WT_D4, (1, 3))    # one level along axis 1, three along axis 2
idwt!(dwt(x, WT_D4, 2), WT_D4, 2) ≈ x

See also dwt, idwt!, wpt!, and the Phase conventions page.

idwt(x, b, l; convention = :aligned)

Return an inverse-transformed copy of x. The copying form of idwt!.

idwt!(x, b, l; convention = :aligned)
idwt!(x, w, b, l; convention = :aligned)

Invert a wavelet transform in place, and return x. The arguments are those of dwt!; convention must match the one used for the forward transform.

wpt(x, b, l; convention = :aligned)

Return a packet-transformed copy of x. The copying form of wpt!.

wpt!(x, b, l; convention = :aligned)

Wavelet packet transform of x in place, and return x. Equivalent to dwt!(x, b, l; wpt = true): the recursion is applied to every subband rather than only the approximation band.

iwpt(x, b, l; convention = :aligned)

Return an inverse packet-transformed copy of x.

iwpt!(x, b, l; convention = :aligned)

Invert a wavelet packet transform in place, and return x.

nsdwt(x, b, l; convention = :aligned)

Return a nonstandard-transformed copy of x.

nsdwt!(x, b, l; convention = :aligned)

Nonstandard discrete wavelet transform of x in place, and return x. The nonstandard transform applies every level along one axis before moving to the next, while dwt! advances one level along every axis at a time. The two agree for a single level and differ deeper. The arguments are those of dwt!.

nsidwt(x, b, l; convention = :aligned)

Return an inverse nonstandard-transformed copy of x.

nsidwt!(x, b, l; convention = :aligned)

Invert a nonstandard wavelet transform in place, and return x.

nswpt(x, b, l; convention = :aligned)

Return a nonstandard packet-transformed copy of x.

nswpt!(x, b, l; convention = :aligned)

Nonstandard wavelet packet transform of x in place, and return x.

nsiwpt(x, b, l; convention = :aligned)

Return an inverse nonstandard packet-transformed copy of x.

nsiwpt!(x, b, l; convention = :aligned)

Invert a nonstandard wavelet packet transform in place, and return x.

rtree_view(x, band)

Return a view of one subband of x, the coefficient array of a wavelet transform.

band is an Int index into rtree_views, or a Symbol or String naming the band. Names read one letter per axis, l for low and h for high, so :ll is the 2-D approximation band, :lh is low along axis 1 and high along axis 2, and :hh is high along both. Indices run from 1 for all-low to 2^ndims(x) for all-high.

A name must have one letter per axis. _ leaves that axis whole, so :l_ is the low band along axis 1 across every axis-2 coefficient, and :_l is the low band along axis 2 across every axis-1 coefficient. This reaches the bands of an asymmetric transform: after dwt!(x, WT_HAAR, (1, 2)), @rtview x[:ll, :_l] is the approximation band one level deep along axis 1 and two along axis 2.

Chaining navigates a multi-level transform. After dwt!(x, WT_D4, 2), rtree_view(rtree_view(x, :ll), :ll), or @rtview x[:ll, :ll], is the coarsest approximation band. The result is a SubArray, so writing to it writes to x.

rtree_views(x)

Return every first-level subband of x as a tuple of views, in the order used by rtree_view. Band 1 is the approximation band and band 2^ndims(x) is high along every axis.

Each view is half the length of x along every axis, so on a multi-level coefficient array rtree_views returns the bands of the finest level.

@rtview x[band]

Rewrite subband indexing into calls to rtree_view. The result is a SubArray, so the macro works on either side of an assignment or a broadcast:

@rtview a[:ll] .= 0       # zero the approximation band
b = @rtview a[:hl, :hh]   # HH band of the HL band
@rtview a[:ll, :_l]       # axis 2 twice down, axis 1 once

Multiple bands chain from the outside in, so _ keeps whatever range an earlier name established. Only references whose indices are all band names (Symbol, String, or QuoteNode) are rewritten; a[1] and a[1:4] are left alone.

WTBasis

Abstract supertype for wavelet bases. The only concrete basis is WTOrthogonalBasis.

WTOrthogonalBasis(φ, ψ)
WTOrthogonalBasis(; φ = nothing, ψ = nothing)

An orthogonal wavelet basis, stored as a scaling filter φ and a wavelet filter ψ of equal-length SVectors. Both filters are normalised on construction.

The keyword form needs only one filter: the other is derived as its orthogonal complement with complement. The positional form takes both.

Examples

using StaticArrays

const MY_HAAR   = WTOrthogonalBasis(φ = SA[1.0, 1.0])
const MY_HAAR_ψ = WTOrthogonalBasis(ψ = SA[1.0, -1.0])

See the Bases page for the built-in bases.

complement(ψ)

Return the orthogonal complement of the wavelet filter ψ, by the reverse-and-multiply trick.

threshold!(x, λ; mode = :soft)

Shrink every entry of x toward zero in place, and return x. x is usually a coefficient array from dwt. mode selects the rule:

• :soft (the default) subtracts λ from the magnitude, sign(x) * max(|x| - λ, 0);
• :hard keeps x where |x| > λ and zeroes it otherwise;
• :garrote applies x * (1 - λ^2 / |x|^2) where |x| > λ.

Examples

c = dwt(s, WT_D4, 4)
threshold!(c, 0.1; mode = :soft)
s_denoised = idwt!(c, WT_D4, 4)

keeplargest!(x, k)

Zero all but the k largest-magnitude entries of x, in place, and return x. Ties at the cut are broken by position, so exactly k entries are retained.

Examples

c = dwt(img, WT_D4, 3)
keeplargest!(c, round(Int, 0.05 * length(c)))
compressed = idwt!(c, WT_D4, 3)

noisiness(x)

Estimate the standard deviation of additive Gaussian noise from an array of detail coefficients, as median(abs, x) / 0.6745.

Pass the finest-scale detail coefficients, which carry almost no signal. denoise picks that band for you.

rtenergy(c)

Energy sum(abs2, ...) of every first-level subband of the coefficient array c, as a NamedTuple keyed by band name: :ll, :lh, :hl, :hh in two dimensions, :l and :h in one.

The bands are those of rtree_views, so recursing into :ll gives the next level.

Examples

c = dwt(img, WT_D4, 2)
e = rtenergy(c)
e.ll / sum(abs2, c)     # share of the energy in the approximation band

sparsity(x)

Hoyer's sparsity of x, between 0 for a uniform array and 1 for a single nonzero entry. Apply it to a coefficient array to compare how well two bases concentrate the same signal.

denoise(x, b, l; rule = :universal, mode = :soft, sigma = nothing, cycles = 0)

Denoise x by thresholding its wavelet coefficients and inverting.

rule selects the threshold: :universal is the VisuShrink threshold σ sqrt(2 log n), :sure is Stein's unbiased risk estimate, and a number is used as the threshold directly. The noise σ is estimated from the finest-scale detail coefficients with noisiness unless the sigma keyword supplies it. mode is passed to threshold!.

cycles > 0 runs the whole denoiser under that many circular shifts with cyclespinning!, which suppresses boundary artefacts.

Examples

s_denoised = denoise(s, WT_D4, 4)
s_smooth   = denoise(s, WT_D4, 4; cycles = 8)

compress(x, b, l; keep = 0.05)

Compress x by keeping the largest keep fraction of its wavelet coefficients and inverting. Returns a NamedTuple with the reconstruction x, the number of coefficients kept, and the fraction of the energy energy retained.

Examples

result = compress(img, WT_D4, 3; keep = 0.05)
result.x        # the reconstruction
result.energy   # fraction of the energy kept

cascade(basis, iterations) -> (φ, ψ)

Sample the scaling function φ and the wavelet function ψ of basis with the cascade algorithm.

The cascade starts from the box function and applies iterations subdivisions of the refinement equation φ(t) = √2 Σ_k h[k] φ(2t - k), where h is the scaling filter. Each step doubles the resolution, so the vectors are sampled at k / 2^iterations; the nonzero samples cover the support [0, N - 1] of the filters, where N is the number of taps. The wavelet function is one high-pass subdivision of the previous scaling iterate, so both vectors have the same length.

iterations must be at least 1. Six to eight iterations are enough to plot most bases.

Examples

φ, ψ = cascade(WT_D4, 8)
t = (0:length(φ) - 1) ./ 2^8

See the Algorithms and Bases pages.

cyclespin!(x, n)

Circularly shift x in place by n samples, wrapping the tail to the front. A matrix shifts both axes by the same amount.

cyclespinning!(f, x, n = 4; start = 16)

Average f over n circular shifts of x, in place, and return x. Each shift is undone before the results are averaged, so f always sees a shifted copy of the original: this is the Coifman-Donoho cycle-spinning estimator. The shifts are the primes prime(start + 1) upward. Averaging cancels the edge artefacts that the periodic transform leaves behind, at a cost linear in n.