Math::Wavelet

Wavelet transforms in pure Raku: the single-level and multilevel DWT, the stationary transform, the 2-D transform, thresholding and wavelet denoising, over the 106 discrete wavelets PyWavelets has; and the continuous transform, over its 21 continuous ones — the Gaussian derivatives, the Mexican hat, Morlet, and their complex relatives. It has no dependencies and no native library, and agrees with PyWavelets to ten significant digits.

Version 0.0.2. The interface below is implemented and tested on both engines; what it leaves out is under Scope.

use Math::Wavelet;

my @x = (^64).map({ sin($_ / 4) + ($_ %% 3 ?? 0.1 !! -0.05) });

my ($cA, $cD) = dwt(@x, 'db4');                 # one level
my @y         = idwt($cA, $cD, 'db4');          # and back

my @c = wavedec(@x, 'sym5', :level(3));         # [cA3, cD3, cD2, cD1]
my @r = waverec(@c, 'sym5');
my @t = threshold(@c, 0.1, :mode<hard>);

my @clean = denoise(@x, 'db4');                 # VisuShrink, soft

my @image = (^8).map(-> $r { (^8).map({ $r * $_ }).Array });
my ($a, ($h, $v, $d)) = dwt2(@image, 'haar');   # an image is an array of rows

my ($W, $f) = cwt(@x, [1, 2, 4, 8], 'mexh');    # a row of coefficients per scale

rakudo -Ilib examples/denoise.raku
rakupp -Ilib examples/denoise.raku

Install

Either installer takes it, into the same ~/.raku store:

rakupp install Math::Wavelet
zef install Math::Wavelet

The modules

use Math::Wavelet exports every routine below. The modules behind it can also be used one at a time:

Wavelets

A discrete wavelet is given by name, or as a Math::Wavelet::Filter:

my $w = wavelet('sym4');
say $w.family;                     # Symlets
say $w.dec-len;                    # 8
say $w.vanishing-moments-psi;      # 4
say wavelist('coif').elems;        # 17
say wavelist.elems;                # 127: wavelist(:kind<discrete>) is 106

wavelist takes PyWavelets' :kind, all (the default), discrete or continuous, and wavelet answers a Math::Wavelet::Continuous for a continuous name. The continuous wavelets are under The continuous transform.

A wavelet of your own is either the orthogonal filter bank that one lowpass filter defines, or all four filters given explicitly:

my $own  = Math::Wavelet::Filter.from-lowpass(@dec-lo, :name<mine>);
my $bank = Math::Wavelet::Filter.new(:@dec-lo, :@dec-hi, :@rec-lo, :@rec-hi);

The coefficients are PyWavelets' own doubles, bit for bit. dmey is a 62-tap truncation of the Meyer filter, so it reconstructs only to within about half a percent. The same is true in PyWavelets.

The one-dimensional transform

Every routine takes a name or a Math::Wavelet::Filter for $w, and returns Arrays of Num. Outside periodization, a level of n samples produces ⌊(n + F − 1) / 2⌋ coefficients in each half for an F-tap filter, and the inverse gives back 2·len − F + 2 samples. waverec drops the extra trailing approximation coefficient that an odd length leaves, the same way PyWavelets does.

Signal extension

What the transform assumes lies beyond each end of the signal. The default is symmetric.

periodization is the mode in which an orthogonal wavelet is an orthogonal transform: the coefficients carry exactly the signal's energy and there are exactly as many of them as samples. Every other mode adds F − 2 coefficients per level for the boundary, or F − 1 when the length is odd.

The stationary transform

The undecimated transform keeps every coefficient at every level, so it is shift-invariant, at the cost of n coefficients per band per level.

:trim-approx returns [cA_n, cD_n, …, cD_1] instead. :norm scales the filters by 1/√2, which makes an orthogonal wavelet's transform preserve energy. The length must be divisible by 2^(level + start-level).

Two dimensions

An image is an array of rows. dwt2 filters down the columns first, then along the rows. That is PyWavelets' order, and it gives PyWavelets' band names:

Thresholding and denoising

The threshold modes are PyWavelets':

denoise is VisuShrink. It estimates σ from the finest detail level, thresholds every detail level at σ·√(2 ln n) (or at :threshold, if given), leaves the approximation alone, and reconstructs n samples. :mode is the threshold mode, and :extension the signal-extension mode.

Examples

• examples/denoise.raku — a sine with a step, buried in noise, recovered with three wavelets, with σ estimated and the error before and after.
• examples/compress.raku — keeps the largest tenth of the coefficients and measures how much of the signal six wavelets keep with them.
• examples/image.raku — a two-level 2-D decomposition of a synthetic image with a horizontal and a vertical edge: where its energy goes, and which band sees which edge.
• examples/scalogram.raku — a chirp followed along the ridge of its scalogram: at each moment, the frequency where |W| peaks, against the one the chirp has.

The continuous transform

$scales is one positive number or a list of them, not necessarily integers. A coefficient is a Num, or a Complex for a complex wavelet. The frequencies are in cycles per $sampling-period, which defaults to one sample.

The algorithm is PyWavelets': ψ is sampled at 2¹⁰ points and integrated once, the integral is resampled at each scale, convolved with the signal, and differenced. :method<conv> (the default) convolves directly, which costs n·m for a kernel of m taps, and m grows with the scale. :method<fft> convolves through a radix-2 FFT. The signal's spectrum is computed once per transform size, and a complex kernel costs no more than a real one. The two methods agree to rounding. For 40 scales over 512 samples of mexh, conv takes 3.2 s on Raku++ and 1.5 s on Rakudo, and fft takes 1.1 s on Raku++ and 1.3 s on Rakudo.

Continuous wavelets

B is the bandwidth and C the centre frequency, and M is the spline order, an integer of at least 1: cmor1.5-1.0, fbsp2-1-0.5. A bare shan, fbsp or cmor takes PyWavelets' defaults, which PyWavelets itself now warns about. The explicit form is the one to write.

my $w = wavelet('cmor1.5-1.0');
say $w.family;                     # Complex Morlet wavelets
say $w.complex-cwt;                # True
say $w.bandwidth-frequency;        # 1.5
my ($psi, $x) = $w.wavefun(:level(8));   # 256 samples over [-8, 8]
say $w.psi(0);                     # 0.4606+0i

The discrete transforms refuse a continuous wavelet, and cwt refuses a discrete one, saying so.

Scope

Left out of 0.0.1 on purpose:

The stationary transform inverts only from start level 0, as in PyWavelets.

Compatibility

Neither version is an established floor; no older engine has been tried. The expected values in t/vectors come from PyWavelets 1.8.0, through tools/gen-vectors.py, and the filter table from the same version, through tools/gen-table.py. continuous.vec and cwt.vec came later, from the same PyWavelets on NumPy 2.5.4. That build regenerates the older files with a few last-digit differences, such as 0 against -2.4e-16, so those files were kept as they were.

Author

Andrew Shitov (zef:ash).

Licence

Artistic-2.0.

Why the module is shaped this way, and what running it under two engines turned up, is in notes/Math-Wavelet.md.