The WaveD Transform in R : Performs Fast Translation-Invariant Wavelet Deconvolution

This paper provides an introduction to a software package called waved making avail-able all code necessary for reproducing the ﬁgures in the recently published articles on the WaveD transform for wavelet deconvolution of noisy signals. The forward WaveD transforms and their inverses can be computed using any wavelet from the Meyer family. The WaveD coeﬃcients can be depicted according to time and resolution in several ways for data analysis. The algorithm which implements the translation invariant WaveD transform takes full advantage of the fast Fourier transform (FFT) and runs in O ( n (log n ) 2 ) steps only. The waved package includes functions to perform thresholding and ﬁne resolution tuning according to methods in the literature as well as newly designed visual and statistical tools for assessing WaveD ﬁts. We give a waved tutorial session and review benchmark examples of noisy convolutions to illustrate the non-linear adaptive properties of wavelet deconvolution.


Introduction
In this paper we present the WaveD transform in R and illustrate some statistical applications of the WaveD transform to the deconvolution of noisy signals. The aim of deconvolution is to recover an unknown function f from a noisy observation of g * f (t) = T g(t − u)f (u)du, where the convolution kernel g is observed with noise (taking = ε below) or without noise (taking = 0 below), g (t) = g(t) + ζ(t), t ∈ T = [0, 1], and ξ, ζ are independent white noises and 0 < ε, < 1 are noise levels. Both f and g are supposed to be periodic on T and g * f (t) denotes the circular convolution. In the finite sample implementation of model (1) at points t i = i/n, i = 1, ..., n, we let where σ is the noise standard deviation and n is the sample size. Let y be a vector with elements (y 1 , ..., y n ) where y i = Y (t i ), i = 1, . . . , n. An illustration of model (1) is given in Figure 3 using the test functions of Figure 1. As for model (2) we consider the two cases: (a) = 0 in which case g (t) = g(t) (known kernel); (b) = ε = σ/ √ n (noisy kernel). We denote g a vector (g 1 , ..., g n ) with elements g i = g (t i ), i = 1, . . . , n. An illustration of model (2) in the Fourier domain is given in Figure 8. In its simplest form, the WaveD transform as discussed in this paper requires only y and g as input, other arguments are optional.
• WaveD is capable of representing functions with discontinuities or with non-homogeneous time and frequency behaviour.
• The translation invariant version of WaveD improves upon the numerical performance of ordinary WaveD by cycle spinning over all circulant shifts.
• WaveD is an non-iterative deconvolution technique.

What's new in the waved package?
Earlier versions of the WaveD method have been implemented through various small Matlab packages, corresponding to various existing WaveD transforms. For example one package uses the algorithm of Kolaczyk (1994) to compute the ordinary Meyer wavelet transform. Another uses the algorithm of Donoho and Raimondo (2004) to compute the translation invariant Meyer transform. These small WaveD packages are not self-contained and are intented for use with WaveLab (Buckheit, Donoho, Johnstone, and Sargle 1995).
This paper describes a unified setting where all the WaveD transforms are implemented in the software environment R (R Development Core Team 2007) via a contributed package named waved (Raimondo and Stewart 2006). The aims of the waved package are: 1. To make available, in one self-contained package all code necessary to compute the various WaveD transforms with optimal data-driven tuning for wavelet deconvolution.
2. To take full advantage of the object-oriented R environment: the (top) function, called WaveD, produces objects of class wvd. The wvd class of objects are R lists containing the various WaveD transforms as well as all the WaveD estimate characteristics such as threshold, fine resolution level, degree of Meyer wavelet and so on.
3. To introduce visual and statistical tools to assess the validity and the quality of a WaveD fit. Special features of the waved package include a summary and a plot function specifically designed for objects of class wvd.
4. To allow a user to reproduce illustrative figures and analyses from the literature.
Finally, we discuss how the waved package differs from existing R packages for wavelet analysis. Existing wavelet R packages include: the wavethresh package of Nason, Kovac, and Mächler (2006): a software to perform wavelet statistics and transforms; the waveslim package of Whitcher (2006): basic wavelet routines for one, two and three dimensional signal analysis; the wavelets package of Aldrich (2007): a package of functions for computing wavelet filters. These packages offer a wide range of compactly supported wavelet transforms, typically Daubechies wavelets, for direct data analysis. On the other hand the waved package is designed for indirect data analysis (such as noisy-convolution) and uses band-limited wavelets, typically Meyer wavelets.

Paper organisation
In section 2 we give a brief introduction to the WaveD transform based on the Fourier transform. Section 3 is concerned with setting-up the waved software and its demo. We also present the WaveD function in R and introduce objects of class wvd. In Section 4, we discuss some more advanced features of the WaveD function in R, this includes statistical applications, fine tuning of the parameters and WaveD fit assessment. Section 5 contains a list of waved main functions.

Fourier transforms
Convolution products are naturally represented in the Fourier domain. In the periodic setting, we can write the model (1) in terms of Fourier coefficients, where, with e (t) = e 2πi t and f, g = T fḡ, f = f, e , g = g, e and ξ = ξ, e are i.i.d. standard (complex-valued) normal random variables. As for the model (2) we have where z are i.i.d. standard (complex-valued) Gaussian r.v.'s independent of ξ , and noise level 0 < < 1. This model includes cases where the eigen-values (g ) are not fully known but are observed with noise as illustrated in Figure 8.
Illustration of the WaveD paradigm in the time, Fourier and wavelet domain,

Adaptive denoising via non-linear WaveD transform
In the case of noisy data (1), (2) we note that provides an unbiased estimator of (β κ ) κ . The waved software uses statistical techniques to perform wavelet regression and smoothing. The main idea is to remove small wavelet coefficients (noise) and keep large wavelet coefficients (signal). Optimal and data driven choices of WaveD tuning parameters are further discussed in section 4; here we shall only present the WaveD method in broad terms using a generic threshold function where λ is a threshold parameter. The WaveD estimator  is which with a slight abuse of terminology we also call the WaveD transform. We summarise the main steps in the diagram below

The translation invariant WaveD transform
Numerical (and computational) properties of the WaveD transform are improved using cycle spinning (Donoho and Raimondo 2004). For any h > 0, we denote T h f (x) = f (x + h) the shift operator and → FWaveD , ← IWaveD the Forward WaveD transform and its inverse (5). For an arbitrary time shift h we define one cycle-spin of the WaveD transform as Let H n = {1/n, 2/n, ..., 1 − 1/n, 1} be the set of all possible circulant shifts. The translation invariant WaveD transform is 3. The WaveD transform in R

Software access
The waved software (Raimondo and Stewart 2006) is provided as an R package obtainable from the Comprehensive R Archive Network at http://CRAN.R-project.org/. Installation instructions are provided there also.

Getting help
Once the waved package has been installed detailed help pages for basic functions may be obtained within R using the help() function. For example help(WaveD) gives the help page of the main waved function. Note that waved refers to the R package whereas WaveD is the main function which performs wavelet deconvolution. See section 5 for a list of basic waved functions.

The waved demo
From now on we assume that the waved package has been attached. Typing demo(waved) provides a series of examples which illustrate various applications of the WaveD transform.
To simulate data according to (1) and (A ε ) one needs to specify: (a) a target function f ; (b) a convolution kernel g; (c) a sample size n; (d) a standard deviation σ. The simplest way to get started is to use the waved package demo. Just type demo(waved) and answer questions at the prompt. Alternatively, the function waved.example() can be used (recommended) to generate the data sets and figures in this paper by setting its two arguments to TRUE (default).
In the default setting the convolution kernel g is defined using the density of a Gamma distribution with shape and scale parameters set to 0.5 and 0.25 respectively. In this setting, the eigen-values (g ) satisfy |g | ∼ |l| −ν , with ν = 0.5. The parameter ν which drives the decay of the eigen-values is often referred as the Degree of Ill-Posedness (DIP) of the convolution problem ).

Setting up your examples
Once you become more familiar with the waved package you may want to generate your own data by modifying the default parameters of the demo. The function waved.example() can be used to generate simulated examples with different model parameters: sample size, noise level, Degree of Ill-Posedness, seed and so on. This is done by setting its first argument to FALSE and answering questions at the prompt. For example to set the sample size to n = 4096, R> my.own.simulation <-waved.example(FALSE) Please enter the sample size (must be a power of 2)

R> 4096
and so on to keep or change the other model parameters. To recover the data sets used in this paper just set the first argument of waved.example() to TRUE; the second argument refers to graphics display, the default setting is TRUE, hence Alternatively, the figures can be produced by calling in waved data names R> plot(t, lidar, type="l") R> plot(t, lidar.blur, type="l") R> plot(t, lidar.noisy, type="l")

The WaveD function and wvd objects
The function WaveD creates R objects of class wvd. The wvd class objects are lists which contain the various WaveD transforms as well as all the WaveD estimate characteristics such as threshold, resolution level, degree of the Meyer wavelet and so on. Statistical properties of objects of class wvd are discussed in Section 4. The summary and plot functions for objects of class wvd are discussed in Section 4.5. In its simplest version the WaveD function requires two input arguments: y a vector with elements (y 1 , ..., y n ) where y i = Y (t i ), i = 1, . . . , n, see (1), and g a vector (g 1 , ..., g n ) with elements g i = g (t i ), i = 1, . . . , n, see (2). Optional arguments to WaveD include: F the finest resolution level j used in the expansion (7) as well as the threshold value λ at (6). The parameter F may take any value within the range L, ..., (log 2 (n) − 1) where L is a low resolution level (default L = 3). In our examples n = 2048 so that F may take any value within the range 3,...,10. For illustration purposes, R> lidar.wvd <-WaveD(lidar.blur, g, F=6, thr=0) R> multires(lidar.wvd$w, lo=3, hi=6) R> lidar.noisy.wvd <-WaveD(lidar.noisy, g, F=6, thr=0) R> multires(lidar.noisy.wvd$w, lo=3, hi=6) these commands produce the WaveD transform of the blurred lidar data of Figure 2 as well as the WaveD transform of the noisy-blurred lidar data of Figure 3. In both cases using F=6 as the finest resolution level and threshold thr=0 (no thresholding). The corresponding plots can be seen in Figure 4.
The forward WaveD transform can be obtained by typing R> lidar.w <-FWaveD(lidar.blur, g, F=6) or lidar.wvd$w which returns the same output as lidar.w as defined above. The vector lidar.w is a vector of wavelet coefficients stored from the lowest resolution level to the highest resolution level. The function dyad(j) may be used to access wavelet coefficients at a given resolution level R> lidar.wavelet.coef.at.level.5 <-lidar.w[dyad(5)] A useful property of wavelet coefficients is that they are large (in absolute value) near discontinuities, see e.g. the top RHS plot of Figure 4. Another feature of wavelet coefficients is that they become more and more sensitive to noise as the resolution level increases. See e.g. the top LHS plots of Figure 4. In Figure 4 (top plots), the function multires() is used to depict wavelet coefficients according to time and frequency. More details about the multires() function and the data structure of the WaveD transforms are given in Section 5.
The inverse WaveD transform can be obtained by typing R> inverse.waved <-IWaveD(lidar.w) or lidar.wvd$iw. The vector lidar.wvd$iw returns the inverse WaveD transform (5) computed from lidar.w without any thresholding. Two illustrations of the inverse WaveD transform are depicted on the bottom plots of Figure 4 (with corresponding Forward WaveD transforms depicted on the top plots).
The ordinary WaveD transform is a combination of the FwaveD and IWaveD transforms together with some thresholding options (7). The ordinary WaveD transform is obtained from a wvd object by typing lidar.wvd$ord which, here, returns approximations to lidar as seen in Figure 4. If no thresholding is performed (thr=0) the ordinary WaveD transform returns the same output as the inverse WaveD transform. If a non-zero threshold is used the ordinary WaveD transform returns the inverse WaveD transform after thresholding. For noisy data it is desirable to improve WaveD approximations such as depicted on the RHS bottom plot of Figure 4 by using a non-zero threshold. This is detailed next.

Statistical applications of the WaveD transform
In this section we discuss some more advanced features of the WaveD transform when dealing with noisy data. We use the simulated data of Figure 3 to illustrate how the WaveD function chooses the fine tuning parameters F and thr in a data-driven fashion in agreement with the optimal choices prescribed in the literature (Cavalier and Raimondo 2007).

Choosing a threshold
The threshold value λ at (6) may be thought of as a smoothing parameter since it dictates the amount of smoothing in the estimate, large λ yields smoother estimates and vice-versa. A single threshold value λ may be entered directly in the WaveD function R> plot(t, WaveD(lidar.noisy, g, F=6, thr=0.2)$ord, type="l") R> plot(t, WaveD(lidar.noisy, g, F=6, thr=0.02)$ord, type="l") with corresponding plots depicted in Figure 5. Alternatively a set of level dependent thresholds may be entered as a vector, R> my.thr <-c(0.01, 0.02, 0.03, 0.04) R> lidar.my.thr.wvd <-WaveD(lidar.noisy, g, L=3, F=6, thr=my.thr) will use threshold λ = 0.01 at level j = 3, λ = 0.02 at resolution level j = 4 and so on. The maxiset threshold. If no threshold parameter is specified the WaveD function will use the so-called "maxiset threshold" (11). This level-dependent threshold is derived from the maxiset theory . For example, R> lidar.maxi.wvd <-WaveD(lidar.noisy, g) will use the follwing threshold values R> round(maxithresh(lidar.noisy, g, L=3, F=6), 4) [1] 0.0134 0.0198 0.0298 0.0459 corresponding to a vector thr with entries (λ 3 , λ 4 , ..., λ 6 ) of level dependent thresholds (11). The effect of the maxiset threshold is illustrated in Figure 6. As seen in the RHS plots of Figure 6, the WaveD estimate with the maxiset threshold automatically select significant coefficients to be kept for the reconstruction. This process removes noise (small coefficients) and smoothes the estimate. The un-thresholded and the thresholded Forward WaveD transforms may be obtained from a wvd object as follows, R> unthresholded.w <-lidar.maxi.wvd$w R> multires(unthresholded.w, lo=3, hi=6) R> thresholded.w <-lidar.maxi.wvd$w.thr R> multires(thresholded.w, lo=3, hi=6) which produce the plots of Figure 6. The maxiset threshold is computed as follows • γ: constant which depends on the tail of the noise distribution. For Gaussian noise, the range √ 2 ≤ γ ≤ √ 6 gives good results in practice. The default setting for WaveD is the conservative choice γ = √ 6.
• σ j : level-dependent scaling factor which depends on the convolution kernel.
• c n : sample size-dependent scaling factor reminiscent of the Universal threshold:

Choosing the finest resolution level
The fine resolution level F is related to the highest (Fourier) frequency M allowed in the WaveD estimator 2 F ≈ M. The tuning parameter F stipulates the range of resolution levels where the approximations (7) or (8) are used: Here L is a low resolution parameter (default is L = 3). A numerical value for F within the range L ≤ F ≤ log 2 (n) − 1 may be entered directly in the WaveD function as in the examples of Section 3.5 where we used F = 6. If n = 2048 the maximum value allowed is F=10, R> lidar.Fmax.wvd <-WaveD(lidar.noisy, g, F=10) R> multires(lidar.Fmax.wvd$w.thr) R> plot(t, lidar.Fmax.wvd$ord, type="l") Unlike direct estimation problems (Donoho, Johnstone, Kerkyacharian, and Picard 1995) where it is customary to keep all resolution levels setting F = log 2 (n) − 1, the asymptotic theory for deconvolution  shows that one should stop at a fine resolution level F = j 1 where j 1 depends on the degree of ill-posedness of the convolution kernel (10) The last condition shows that the faster the eigen values go to zero the sooner the wavelet expansion should stop. In pratical terms this means that the maxiset threshold will prevent noise in the estimate up until a high resolution level j 1 which depends on the degree of difficulty of the convolution as well as the noise level. It is important to note that, even after maxiset thresholding, the WaveD estimate based on all resolution levels may, sometimes, contain high noise perturbations. This is illustrated on Figure 7 where there are large noise residuals in the WaveD estimate due to large (unthresholded) coefficients at resolution level F=10.
Data driven fine level selection. To prevent noise perturbation at high resolution level, WaveD is fitted with a function find.j1 which implement the data-driven method of Cavalier and Raimondo (2007) to find the optimal fine resolution level j 1 for noisy deconvolution based on the maxisets threshold. The idea is to keep all (Fourier) frequencies until (the moduli of) the eigen values fall below an appropriate noise level. This is illustrated on Figure 8. Let denote the maximum Fourier frequency allowed in the WaveD formula (5). Then we define the maximum wavelet resolution level aŝ where x is the largest integer below x. As seen on Figure 8 this process works for both noise-free eigen-values ( = 0) and noisy eigen-values. For example, R> print(find.j1(g, scale(lidar.noisy))) [1] 198 6 returns the optimal Fourier frequency M and fine resolution level F. The plotspec90 function produces plots of the fine level selection process (13) R> plotspec(g, scale(lidar.noisy)) R> plotspec(g.noisy, scale(lidar.noisy)) as seen in Figure 8.

Improving the fit using the translation invariant WaveD transform
While thresholding wavelet coefficients reduces the noise and smoothes the WaveD estimate it also introduces Gibbs phenomenon near discontinuities see e.g. RHS bottom plots of Figure 6. Such Gibbs effects can be reduced by cycle spining (Donoho and Raimondo 2004). Both the ordinary (7) and the translation invariant (8) WaveD transform can be obtained from a wvd object, R> plot(t, lidar.maxi.wvd$ord, type="l") R> plot(t, lidar.maxi.wvd$waved, type="l") with corresponding outputs depicted in Figure 10. In any case where some thresholding is performed we recommend using the translation invariant WaveD transform as it reduces visual artifacts.
The MC-option. The algorithm which implements the translation-invariant WaveD transform takes full advantage of the Fast Fourier Transform and requires only O(n(log n) 2 ) steps. This is faster than the algorithm which implements the ordinary WaveD transform.

Thresholding policy
There are many ways to threshold wavelet coefficients and different strategies may be used . The two main thresholding policies studied in the literature are the HARD thresholding policy as in (6) or the SOFT thresholding policy, where λ is a threshold parameter. The statistical theory for WaveD estimation ) is established for the HARD threshold policy (6) which is the default setting in WaveD. However, for data analysis purposes and experimental study we provide a SOFT thresholding option (15) in the WaveD function, Normal Q−Q Plot q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q −4 −2 0 2 4 Figure 13: plot(lidarT.wvd). These residual plots suggest a poor WaveD fit.

Boxplot
R> plot(t, WaveD(lidar.noisy, g, SOFT=FALSE)$ord, type="l") R> plot(t, WaveD(lidar.noisy, g, SOFT=TRUE)$ord, type="l") As seen in Figure 10, SOFT thresholding tends to further smooth the WaveD estimate but the general appearance does not appear as sharp as the translation invariant WaveD estimate of Figure 9.

The summary and plot functions for wvd objects
For convenience we provide a summary and a plot function specifically for objects of class wvd. These functions can be used to assess the quality of WaveD fits. We illustrate this using the doppler.noisy data set (doppler in Gaussian noise, see Figure 11) and the lidar.noisyT data set (lidar in Student-t 2 noise, see Figure 14).
R> doppler.wvd <-WaveD(doppler.noisy, g) R> lidarT.wvd <-WaveD(lidar.noisyT, g) R> plot(doppler.wvd) In addition to providing the tuning parameters F = j 1 , thr = threshold, M = maximum Fourier frequency, γ = maxiset threshold noise constant as well as thresholding policy, the summary function gives some additional statistics such as the percentage of thresholding at a given resolution level as well as the maximum (in absolute value) of the wavelet coefficients at a given resolution level. It also gives the result of a test for normality based on the estimated noise in the data. This can be used to assess the WaveD fit as discussed next.
Estimating noise contribution. In statistical application of wavelet methods it is customary to estimate noise feature such as variance or tail index using the wavelet coefficients of the raw data at the largest resolution level (Raimondo and Tajvidi 2004). Here we shall call the vector of wavelet coefficients at the largest resolution level: noise.proxy. This vector may be obtained from a wvd object: noise.proxy <-lidar.maxi.wvd$noise. The summary and plot functions use the noise.proxy vector to perform some elementary data analysis, compare Figure 12 with Figure 13.
WaveD-fit assessment. The asymptotic theory and fine tuning of the WaveD parameters (Cavalier and Raimondo 2007) is based on the Gaussian white noise model (1) in which the error terms follow a normal distribution. A close inspection to the proof shows that the constant γ used in the maxiset threshold depends on the tail of the noise. For Gaussian noise the value γ = √ 6 gives good results in simulation. However, in other scenarios a larger value may be needed as this would be the case for heavy tailed noise. To assess the appropriateness of the WaveD fit and of the maxiset threshold, the summary function gives the result of a (Shapiro) test for normality based on the estimated noise in the data. For example, R> plot(lidarT.wvd) R> summary(lidarT.wvd) ...Estimated standard deviation= 0.094 Shapiro test for normality, P= 1.490849e-12 we see that for the Student-t 2 noise scenario, the WaveD fit residuals fail the normality-test with a Shapiro-test P -value close to zero. The corresponding WaveD estimate exhibits a large noise residual even after thresholding as seen on the bottom-RHS plot of Figure 14. This combined with the residual plots of Figure 11 suggest a poor WaveD fit.
Improving WaveD-fit in non-Gaussian noise. We suggest some heuristic approaches to improve WaveD fit in non-Gaussian scenarios, R> plot(WaveD(lidar.noisyT, g, SOFT=TRUE)$ord, type="l") R> plot(WaveD(lidar.noisyT, g, SOFT=TRUE, eta=sqrt(8))$ord, type="l") R> plot(WaveD(lidar.noisyT, g, SOFT=FALSE, eta=sqrt(8))$waved, type="l") R> plot(WaveD(lidar.noisyT, g, SOFT=FALSE, eta=sqrt(12))$waved, type="l") (a) Using a soft threshold tends to reduce noise contributions and is more robust against non-normal noise; (b) using an ordinary WaveD estimate with a slightly larger γ tends to reduce noise contributions and is more robust against non-normal noise (in the summary function check max|w| against the threshold and increase γ accordingly); (c) using a TI-WaveD estimate with a slightly larger γ tends to reduce noise contributions but may not remove residuals contributions as effectively as in Ordinary WaveD; (d) using a TI-WaveD estimate with a bigger γ tends to reduce noise contributions and Gibbs phenomena. These four approaches are illustrated on Figure 15 using the lidar.noisyT data. On this occasion the TI-WaveD estimate with γ = 2 √ 3 yields a better estimate.

WaveD estimation with noisy eigen-values
We finish this section by illustrating further adaptive properties of WaveD estimates. Depicted in Figure 16 is a WaveD lidar fit constructed from the noisy-blurred data of Figure 3 and the noisy eigen-values in the RHS plot of Figure 8.
R> lidar.NEV.wvd <-WaveD(lidar.noisy, g.noisy) R> plot(lidar.NEV.wvd) By comparing Figure 16 with Figure 9 we see that the quality of the WaveD approximation is not affected much if one uses noisy eigen values instead of the true eigen values. This is consistent with the asymptotic theory and numerical results of Cavalier and Raimondo (2007).

The WaveD command
The command WaveD(y,g) performs wavelet deconvolution using the data y and the convolution kernel g. If g is not specified WaveD(y) performs a (direct) wavelet transform.
-MC: if Monte Carlo (MC=TRUE) WaveD returns only the TI-WaveD (default=FALSE). Note if MC=TRUE the WaveD output is a simple vector not a list.
-SOFT: if SOFT=TRUE WaveD uses the soft-thresholding policy else hard (default=FALSE).
• Value: in the case that MC=TRUE, WaveD returns a vector consisting of the translation invariant WaveD estimate (8). In the case that MC=FALSE (the default), WaveD returns an object of class wvd, list with following components -j1: estimate of optimal resolution level (14).
-F: fine resolution level used (may be different than j1).
percent: percent of thresholding per resolution level.
noise: noise proxy, wavelet coefficients at the largest resolution level.
p: P-value of the Shapiro normality test based on noise.
residuals: wavelet coefficients that have been removed before fine level F.

Other useful commands
We give a list of other waved commands which can be used independently of the WaveD() function.
In the examples below it is assumed that y is a vector with elements (y 1 , ..., y n ) where y i = Y (t i ), i = 1, . . . , n, as in (1) and that g: a vector (g 1 , ..., g n ) with elements g i = g (t i ), i = 1, . . . , n, as in (2).
• FWaveD(y, g): the command lidar.w=FWaveD(y, g) returns a vector of wavelet coefficients as in WaveD(y, g)$w. This vector has length n, the last n/2 entries are wavelet coefficients at resolution level (J − 1) where J = log 2 (n); the n/4 entries before that are wavelet coefficients at resolution level (J − 2), and so on until level L. In addition, the first 2 L entries are scaling coefficients at coarse resolution level C = L. See the dyad() function below for how to access wavelet coefficients at a given resolution level.