GeneralizedSampling package¶

Installation¶

The GeneralizedSampling package is available through the Julia package system by running Pkg.add("GeneralizedSampling").

Background¶

When sampling a function/signal, one often has to use a basis that is inferior for reconstruction, i.e., a basis requires a large number of elements to give a decent approximation. Think e.g. of the Gibbs phenomenon that occurs when using the Fourier basis for reconstructing discontinuous functions.

Generalized sampling is a general technique for transforming samples of a function w.r.t. one basis into samples w.r.t. another basis, i.e., essentially performing a change of basis.

The GeneralizedSampling package currently implements transformations from the Fourier to the scaling function basis on \(L^2([-1/2,1/2]^d)\).

Basic Usage¶

The primary element in GeneralizedSampling is a change of basis type (CoB) computed from the sample locations, the name of wavelet used for reconstruction and the scale J of the wavelet space:

T = freq2wave(samples, wavename, J)

The T object behaves in many ways like the ordinary matrix it resemples. In particular, multiplication, multiplication with its adjoint and the backslash operator for least squares solutions work as expected. So if the Fourier transform Ghat of a function G is sampled in the locations used for T, G‘s representation w in the wavelet domain of scale J is approximated by

f = Ghat(samples)
w = T \ f

To evaluate w in the wavelet basis, the IntervalWavelets package can be used:

using IntervalWavelets
y = weval(w, wavename, R)

The third argument R controls the number of points in which each scaling function is evaluated – se the documentation for IntervalWavelets for examples.

A simple example is the Gaussian density, since its Fourier transform is again a Gaussian density. Consider an approximation with \(2^5 = 32\) Daubechies 4 scaling functions from 64 Fourier measurements:

using GeneralizedSampling
xi = GeneralizedSampling.grid(64, 0.5);
T = Freq2Wave(xi, "db4", 5)
f = exp(-xi.^2);
w = T \ f;
using IntervalWavelets
x, y = weval(real(w), "db4", 10);

Using the Plots package this can be plotted with the following commands:

using Plots
plot(x, sqrt(pi)*exp(-(pi*x).^2), label="true")
plot!(x, y, label="approximation")

The undesired behavior in the ends can be resolved by increasing the number of Fourier measurements or using a different wavelet for reconstruction.

Complete examples are found in the examples folder. The file truncated_cosine.jl approximates a truncated cosine in the Haar basis and can be called using

tc_path = joinpath(Pkg.dir("GeneralizedSampling"), "examples", "truncated_cosine.jl")
include(tc_path)

Note:

• The theory of generalized sampling promises that the change of basis matrix \(T\) is numerically stable under a number of assumptions:

  • The number of samples must be sufficiently high compared to \(J\).
  • The samples must be in a domain symmetric around the origin.
  • For samples on a grid the distance to neighboring sample points must be less than the inverse length of the reconstruction interval.
  • Non-uniform samples must have a sufficiently high bandwidth and sufficiently low density – see [AGH14] and [AGH15] for further details.

• The condition number of \(T\) is not directly available. To compute the condition number the change of basis matrix has to be computed explicitly with collect(T).

• The change of basis matrix may very well be too large to compute explicitly; check size(T) before collecting.

Wavelets¶

Currently GeneralizedSampling supports reconstruction in Daubechies wavelets/scaling functions. As the reconstruction happens on \([-1/2,1/2]\) the functions near the boundaries needs to be modified – which can happen in multiple ways. We have chosen the boundary wavelets from [CDV93], which has the same number of vanishing moments as the internal/non-boundary wavelets.

The allowed wavelets are named “haar”, “db1”, “db2”, ..., “db8”.

Change of basis¶

freq2wave(samples, wavename, J, B, ...)¶

Compute a frequency-to-wavelet change of basis object.

• samples are the sampling locations. For 1D samples it a vector and for 2D samples it is matrix with 2 columns.
• wavename is a string as described in in the Wavelets section.
• J is the scale of the reconstruction wavelets (the \(V_J\) space in multiresolution terminology). Note that \(2^J\) has to be larger than the length of the wavelet’s support.
• If samples are not a uniform grid, a bandwidth B has to be supplied that is larger than maxabs(samples). Note that if B is too large the density of the samples may also be too large, which degenerates the condition number of T
• The optional arguments are passed to the functions computing Fourier transforms of the wavelet (if needed).

As mentioned in Getting Started, the benefit of generalized sampling is that the computations are numerically stable. However, some assumptions must be fulfilled:

• There should be more samples than reconstruction coefficients. The ratio between samples and the number of reconstruction coefficients that ensures a numerically stable matrix is called the stable sampling rate. For uniform samples the stable sampling rate is well described – see References. For non-uniform samples the stable sampling rate also depends on the density of the samples, which is defined as the minimum radius that gives a covering of the bandwidth area with equal sized circles centered at the sampling points.
• The samples should be distributed around the origin, i.e., only positive samples does not work.

Example:

samples = grid(2^7, 0.5)
T = freq2wave(samples, "db2", 6)

collect(T)¶

Compute the explicit change of basis matrix from T. Note that this is time consuming and possible impossible to hold in memory for large sampling/reconstrution sets.

\

Overload of the usual backslash operator: x = T \ b computes the least squares solution of \(Tx \approx b\).

*, '*

Overload of multiplication with T and T', the adjoint of T.

size(T)¶

Get the size of T as a tuple \((M,N)\).

size(T, d)

Get the size along dimension d of T.

wsize(T)¶

Get the size of the reconstructed wavelets as a tuple. In 1D the result is \((N,)\) and in 2D the result is \((N,N)\).

isuniform(T)¶

Returns true if the samples used for T are uniform and false otherwise.

hasboundary(T)¶

Returns true if the wavelet used for reconstruction in T has special functions near boundaries and false otherwise.

van_moment(T)¶

Get the number of vanishing moments of the wavelet used for reconstruction in T.

Types¶

The abstract change of basis supertype is denoted CoB.

The specific change of basis types implemented are from Fourier to wavelet bases. They are collectively denoted Freq2Wave and are a subtype of CoB:

Freq2Wave <: CoB

The computations for wavelets with boundary correction are more involved than for those without and therefore two subtypes of Freq2Wave are introduced for both 1D and 2D:

Freq2NoBoundaryWave1D <: Freq2Wave
Freq2BoundaryWave2D <: Freq2Wave
Freq2NoBoundaryWave1D <: Freq2Wave
Freq2BoundaryWave2D <: Freq2Wave

Fourier transforms¶

Fourier transforms of the scaling functions are available. The high level interface is

FourScalingFunc(xi, wavename, J, k; ...)¶

Evaluate the Fourier transform of wavename at xi.

• xi is either a real number of an array of real numbers.
• wavename is a string as described in in the Wavelets section.
• Optional J is the scale of the scaling function, which by default is 0.
• Optional k is the translation of the scaling function, which by default is 0.

The remaining arguments relate to the iterative computations of the Fourier transforms and are usually not needed. Check the inside documentation for more info.

As an example, the following command computes the Fourier transform of the Daubechies 2 scaling function and plots the real and imaginary part using Winston:

x = linspace(-5, 5, 1000)
y = FourDaubScaling(x, "db2")
using Winston
plot(x, real(y), x, imag(y))

FourScalingFunc(xi, wavename, side, J, k; ...)

As above, but for the boundary scaling functions. side is either 'L' or 'R'.

Note that these Fourier transforms are for the scaling functions that in the time domain are translated to fir the reconstruction interval \([-1/2,1/2]\), i.e., their Fourier transforms are phase shifted.

The lower level functions are available for each type of scaling function, but not documented here. Check the documentation in Julia with the usual ?function where function is FourHaarScaling or FourDaubScaling.

Miscellaneous¶

Functions that are used for internal documentation are not documented here; they all have documentation available from within Julia.

To generate sampling locations from a uniformly spaced grid there are functions in 1D and 2D.

grid(M, D)¶

Return a vector of M locations evenly distributed around the origin with distance D. By default, D = 1.

grid((M, N), D)

Return a matrix with 2 columns containing the x- and y-values of a uniformly distributed grid of locations around the origin with distance D. There are M different locations in the 1st dimension and N different locations in the 2nd dimension.

isuniform(points)

Returns true if points are located on a uniform grid such as the output from grid and false otherwise.

For a configuration of sampling locations xi the density correcting weights and its density are available as

weights(xi, K)¶

density(xi, K)¶

The bandwidth K is explained in Change of basis and must be at least maxabs(xi).

When dealing with wavelets with boundary corrections, computations differs for the internal and boundary parts. To this end, the split function is available to help divide a vector or matrix of coefficients into the parts related to internal/boundary functions.

split(x, B)¶

Returns three vector slices of the B leftmost, the internal and the B rightmost entries of x, respectively.

split(A, B)

Returns slices of the outer parts of A and its internal parts. The outer parts are each of the four \(B\times B\) corners and each of the four non-corner sides (with one dimension equal to B).