Wavelet shrinkage

A signal in noise, taken apart by a discrete wavelet transform, its small detail coefficients thresholded away, and put back together. The transform is Math::Wavelet, compiled with the rest of this page from Raku into the script your browser runs; only the scalogram below is computed on the server.

signalsignal + noisereconstructed from the kept coefficients

Coefficients, coarse to fine

cA522, never thresholded
cD57 / 22 kept
cD47 / 38 kept
cD36 / 70 kept
cD28 / 133 kept
cD13 / 259 kept

Blue stems are the coefficients that survive the threshold, grey ones are set to zero, and the dashed lines are ± the threshold. The coarse approximation cA is never thresholded. Hard thresholding keeps a survivor as it is, soft shrinks it toward zero by the threshold, and garrote does something in between.

ψ, the wavelet
φ, its scaling function

Daubechies, 8 taps, 4 vanishing moments

The same, in Raku

use Math::Wavelet;

# @x: the doppler signal, 512 samples, plus Gaussian noise of σ 0.2 (seed 1)
# [cA5, cD5, …, cD1]
my @c = wavedec(@x, 'db4', :level(5));

# the noise, estimated from the finest details, and VisuShrink's threshold × 1.00
my $sigma = noise-sigma(@c[*-1]);
my $thr = 1.00 * $sigma * sqrt(2 * log(+@x));

# every detail level thresholded, the approximation left alone
@c[1..*] = @c[1..*].map({ threshold($_, $thr, :mode<hard>) });
my @denoised = waverec(@c, 'db4').head(+@x);

The continuous transform

Each row is one scale, and time runs left to right; brighter is a larger |W|. Moving a wavelet along the signal at every scale is what the transform does, so a short burst shows up as a bright column, and a steady tone as a bright row.

scalogram
top row: a period of 4 samples, scale 2.00bottom row: 80 samples, scale 40.0

Complex Morlet wavelets — complex, so |W| is the envelope; centre frequency 1.0000

ψ on [-8, 8]; solid is its real part, dashed its imaginary part

The same, in Raku

use Math::Wavelet;

# @x: the noisy signal, 512 samples, averaged in pairs: a sketch at half the rate
my @h = @x.rotor(2).map(*.sum / 2);

# 24 frequencies in cycles per pair, from a period of 4 samples down to one of 80
my @f = (^24).map({ 0.5 * (2/80 / 0.5) ** ($_ / 23) });

# each as a scale, through the wavelet's centre frequency at scale 1
my $cf = central-frequency('cmor1.5-1.0', :precision(10));
my ($W, $freq) = cwt(@h, @f.map({ $cf / $_ }), 'cmor1.5-1.0', :method<fft>, :sampling-period(2));

# |W|, one row per scale: the picture above
my @image = $W.map(*.map(*.abs));