On the stable sampling rate for binary measurements and wavelet reconstruction

TL;DR

This paper demonstrates that binary measurements, modeled with Walsh functions and Hadamard matrices, have a linear stable sampling rate for wavelet reconstruction, making them as efficient as Fourier samples for infinite-dimensional signal recovery.

Contribution

It establishes the linear stable sampling rate for binary measurements in wavelet reconstruction, enabling the use of binary sampling in advanced reconstruction techniques.

Findings

Binary measurements have a linear stable sampling rate.

Binary measurements are as efficient as Fourier samples with wavelets.

The results facilitate the use of binary measurements in generalized sampling methods.

Abstract

This paper is concerned with the problem of reconstructing an infinite-dimensional signal from a limited number of linear measurements. In particular, we show that for binary measurements (modelled with Walsh functions and Hadamard matrices) and wavelet reconstruction the stable sampling rate is linear. This implies that binary measurements are as efficient as Fourier samples when using wavelets as the reconstruction space. Powerful techniques for reconstructions include generalized sampling and its compressed versions, as well as recent methods based on data assimilation. Common to these methods is that the reconstruction quality depends highly on the subspace angle between the sampling and the reconstruction space, which is dictated by the stable sampling rate. As a result of the theory provided in this paper, these methods can now easily use binary measurements and wavelet reconstruction bases.

Equations

$$⟨P_{{\mathcal{S}}_M}{G}_{N,M}(f),r_j⟩=⟨P_{{\mathcal{S}}_M}f,r_j⟩,r_j\in {\mathcal{R}}_N,$$

$$\Vert f-G_{N,M}(f)\Vert \le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\Vert f-P_{{\mathcal{R}}_N}f\Vert ,$$

$$\cos (\omega (U,V)):=\frac{1}{\mu (U,V)}:=\underset{u\in U,\Vert u\Vert =1}{\inf }\Vert P_{V}u\Vert ,$$

$$F_{N,M}(f)=\underset{u\in P_{{\mathcal{S}}_M}f+{\mathcal{S}}_M^{\perp }}{\mathrm{argmin}}\Vert u-P_{{\mathcal{R}}_N}u\Vert ,$$

$$\Vert f-F_{N,M}(f)\Vert \le \mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\mathrm{dist}(f,{\mathcal{R}}_N\oplus ({\mathcal{S}}_M\cap {\mathcal{R}}_N{)}^{\perp }).$$

$$\mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\le \theta ,\theta \in (1,\infty ).$$

$$\Theta (N,\theta )=\min \left\{M\in \mathbb{N}:\mu ({\mathcal{R}}_N,{\mathcal{S}}_M)\le \theta \right\}.$$

$$\sum _{j=k}^r{x}_j:=\sum _{j_1=k_1}^{r_1}\dots \sum _{j_d=k_d}^{r_d}{x}_{j_1,\dots ,j_d},$$

$$f(x)=f(x_1)\otimes \dots \otimes f(x_d)\text{ (d-times)}.$$

$$x=\sum _{i\in \mathbb{Z}}{x}_i{2}^i,$$

$$x\oplus y=\sum _{i\in \mathbb{Z}}(x_i{\oplus }_2{y}_i)2^i,$$

$${\omega }_W=\left(\begin{array}{ccccc}0& \cdots & 0& 1& 1\\ ⋮& \text{udots}& \text{udots}& 1& 0\\ 0& \text{udots}& \text{udots}& \text{udots}& ⋮\\ 1& 1& \text{udots}& & ⋮\\ 1& 0& \cdots & \cdots & 0\end{array}\right).$$

$$\mathrm{wal}(s;x)=(-1{)}^{s^n\cdot {\omega }_W{x}^n}.$$

$${\omega }_{WP}=\left(\begin{array}{ccccc}0& \cdots & 0& 0& 1\\ ⋮& \text{udots}& \text{udots}& 1& 0\\ 0& \text{udots}& \text{udots}& \text{udots}& 0\\ 0& 1& \text{udots}& \text{udots}& ⋮\\ 1& 0& 0& \cdots & 0\end{array}\right).$$

$$\mathrm{Wal}(s,x)=(-1{)}^{s_0{x}_0}\mathrm{wal}(\left[s\right];x)\mathrm{wal}(\left[x\right];s),$$

$$\mathrm{Wal}(-s,x):=-\mathrm{Wal}(s,x)$$

$$\mathrm{Wal}(s,-x):=-\mathrm{Wal}(s,x).$$

$$\mathrm{Wal}(s,x)=\underset{k=1}{\overset{d}{⨂}}\mathrm{Wal}(s_k,x_k).$$

$$W{P}_n=\left\{\sum _{j=A}^B{\alpha }_j\mathrm{Wal}(j,z),{\alpha }_{j_i}\in \mathbb{R},A,B\in {\mathbb{Z}}^d,|B|\le n\right\}.$$

$$\text{stackon}[2pt]f\text{stretchto}\text{scalerel}\ast [\text{widthof}f]\bigwedge {0.5ex}_{\text{belowbaseline}[-2pt]W}(s)=⟨f(\cdot ),\mathrm{Wal}(s,\cdot )⟩={\int }_{[0,1{]}^d}f(x)\mathrm{Wal}(s,x)dx,s\in {\mathbb{R}}^d.$$

$$X_j=\frac{1}{N}\sum _{k=0}^{N-1}{x}_k\mathrm{Wal}(j,\frac{k}{N}).$$

$$X_j=\frac{1}{{\prod }_{i=1}^d{N}_i}\sum _{k=0}^{N-1}{x}_k\mathrm{Wal}(j,\frac{k}{N}).$$

$$\mathrm{Wal}(s,x)=\mathrm{Wal}(x,s)\text{ for all }s,x\in \mathbb{R},$$

$$\mathrm{Wal}(2^{k}s,x)=\mathrm{Wal}(s,2^{k}x)\text{ for all }s,x\in \mathbb{R},k\in \mathbb{N}$$

$$\mathrm{Wal}(s,x)\mathrm{Wal}(s,t)=\mathrm{Wal}(s,x\oplus t)\text{ for all }s,x\in \mathbb{R}.$$

$$\mathcal{W}\left\{af(t)+bg(t)\right\}=a\mathcal{W}\left\{f(t)\right\}+b\mathcal{W}\left\{g(t)\right\}\text{ for all }a,b\in \mathbb{R}\text{ and }f,g\in L^2([0,1{]}^d),$$

$$\mathcal{W}\left\{f(t\oplus x)\right\}(s)=\mathcal{W}\left\{f(t)\right\}(s)\mathrm{Wal}(x,s)\text{ for all }x\in {\mathbb{R}}^d\text{ and }f\in L^2([0,1{]}^d)$$

$$\mathcal{W}\left\{f(2^{m}t)\right\}(s)=\frac{1}{2^m}\mathcal{W}\left\{f(t)\right\}(\frac{s}{2^m})\text{ for all }m\in {\mathbb{N}}^d\text{ and }f\in L^2([0,1{]}^d).$$

$${\tilde{\varphi }}_n^{\text{left}}(x)=\sum _{l=0}^{2p-2}\left(\genfrac{}{}{0ex}{}{l}{n}\right)\varphi (x+l-p+1).$$

$${\tilde{\varphi }}_n^{\text{right}}(x)={\tilde{\varphi }}_{-1-n}^{\text{left}}(-x).$$

$$V_j^b=\mathrm{span}\left\{{\varphi }_{j,n}^b:n=0,\dots 2^j-1\right\},$$

Full text

1. Introduction

Reconstructing infinite-dimensional signals from a limited number of linear measurements is a key problem in sampling and approximation theory, and has received substantial attention over the last decades due to its many applications. The list of fields is comprehensive and includes Magnetic Resonance Imaging (MRI) [41, 30], electron tomography [38, 39], lensless cameras, fluorescence microscopy [51, 49], X-ray computed tomography [17, 48], surface scattering [35] among others. Efficient methods for such problems date back to Shannon’s sampling theorem [34, 50, 52], however, over the years, more modern approaches have been developed. Indeed, new methods include generalized sampling, which has been studied by Adcock, Hansen, Hrycak, Gröchenig, Kutyniok, Ma, Poon, Shadrin and others [1, 4, 2, 33, 32, 6, 42], its compressed versions investigated by Adcock, Hansen, Kutyniok, Lim, Poon and Roman [7, 37, 5, 47] as well as the predecessor; consistent sampling, analysed by Aldroubi, Eldar, Unser and others [8, 53, 23, 24, 25, 26]. Note that consistent sampling is very much related to the finite section method [13, 29, 31, 40]. More recently, new methods based on data assimilation have been successfully developed and analysed. A first approach for the same number of measurements $M$ and reconstructed coefficients $N$ was introduced under the name generalized empirical interpolation method by Maday, Patera, Penn and Yano in [44]. This was then further extended to Parametrized-Background Data-Weak (PBDW) approach with $M\geq N$ in [45, 43]. The PBDW approach was additionally analysed and shown to be optimal by Binev, Cohen, Dahmen, DeVore, Petrova, and Wojtaszczyk [10, 21, 11].

The problem is given as follows. An element $f\in\mathcal{H}$, where $\mathcal{H}$ is a separable Hilbert space, is to be reconstructed from measurements with linear functionals $(m_{i})_{i\in\mathbb{N}}:\mathcal{H}\rightarrow\mathbb{C}$ that can be represented by elements $s_{i}\in\mathcal{H}$ as $m_{i}(f)=\langle f,s_{i}\rangle$. The key issue is that the $m_{i}$ cannot be chosen freely, but are dictated by the modality of the sampling device, for example a Magnetic Resonance Imaging (MRI) scanner providing Fourier samples or a fluorescence microscope giving binary measurements. The goal is to reconstruct $f$ from the finite number of samples $\{m_{i}(f)\}_{i=1}^{M}$ for some $M\in\mathbb{N}$. The space of the functions $s_{i}$ is called the sampling space and is denoted by $\mathcal{S}=\overline{\operatorname{span}}\{s_{i}:i\in\mathbb{N}\}$, meaning the closure of the span. In practice, one can only acquire a finite number of samples. Therefore, we denote by $\mathcal{S}_{M}=\operatorname{span}\{s_{i}:i=1,\ldots,M\}$ the sampling space of the first $M$ elements. The reconstruction is typically done via a reconstruction space denoted by $\mathcal{R}$ and spanned by reconstruction functions $(r_{i})_{i\in\mathbb{N}}$, i.e. $\mathcal{R}=\overline{\operatorname{span}}\{r_{i}:i\in\mathbb{N}\}$. As in the case of the sampling space, it is impossible to acquire and save an infinite number of reconstruction coefficients. Hence, one has to restrict to a finite reconstruction space, which is denoted by $\mathcal{R}_{N}=\operatorname{span}\{r_{i}:i=1,\ldots,N\}$. The key is that the $r_{i}$ can be tailored to the desired signal type to be reconstructed. For example, spaces spanned by X-lets (wavelets, curvelets, contourlets, shearlets) [19, 20, 36, 14, 15, 16, 22, 46] may be preferable as reconstruction spaces in imaging applications, whereas polynomials may be useful when recovering very smooth functions.

The methods mentioned above can be described as follows: for $f\in\mathcal{H}$ and $N,M\in\mathbb{N}$, we define the reconstruction method of generalized sampling $G_{N,M}:\mathcal{H}\rightarrow\mathcal{R}_{N}$ by

[TABLE]

where $P_{\mathcal{S}_{M}}$ denotes the orthogonal projection on the subspace $\mathcal{S}_{M}$. Note that the stability and accuracy of this method depends on the subspace angle between the sampling and the reconstruction space, i.e.

[TABLE]

where we define the subspace angle between closed subspaces $U,V\in\mathcal{H}$

[TABLE]

$\omega(U,V)\in[0,\pi/2]$. Moreover, the condition number $\kappa$ of $G_{N,M}$ is also given by $\kappa(G_{N,M})=\mu(\mathcal{R}_{N},\mathcal{S}_{M})$. For the PBDW method one calculates

[TABLE]

and it can be shown that the accuracy then depends on subspace angle as follows

[TABLE]

Moreover, this is sharp because the constant $\mu(\mathcal{R}_{N},\mathcal{S}_{M})$ cannot be improved. It is clear that, in both approaches, the key to success lies in the ability to make sure that

[TABLE]

Thus, we need to balance the number of samples $M$ with the number of reconstruction vectors $N$, and this leads to the so-called stable sampling rate:

[TABLE]

The methods above can only be used efficiently when the stable sampling rate is known and reasonable. In particular, numerical calculations of the stable sampling rate are very time consuming. Moreover, if the stable sampling rate is worse than linear, the approximation quality of the reconstruction space must allow for rapid approximation to compensate for the ”slow” sampling rate. Fortunately, it is possible to obtain sharp results on describing $\Theta(N,\theta)$ for popular sampling and reconstruction spaces, and often, especially for the reconstruction with X-lets, one can establish linearity. In this paper we do so for sampling with Walsh functions and reconstructing with wavelets. Remark that this is not always the case. For example in the case of the reconstruction with Legendre polynomials, the stable sampling rate is known to be quadratic for Fourier measurements [33]. Additionally, for the PBDW-method there are approaches where the reconstruction space is fixed by the PDE and the sampling space is chosen adaptively. This is a very different setting from the one discussed here. Details and convergence rates can be found in [12].

1.1. Connection to previous work and novelty of the paper

The stable sampling rate is well understood when the samples $m_{i}(f)=\langle f,s_{i}\rangle$ are Fourier measurements. In other words, the $s_{i}$ are complex exponentials and $m_{i}(f)$ are the Fourier coefficients. In this case the stable sampling rate is linear for many X-lets including wavelets and shearlets. Fourier samples and X-lets are a natural starting point given the vast applications that are based on Fourier measurements (MRI, tomography problems with parallel beam, surface scattering, radio interferometry etc.). However, the next question regards binary measurements. By binary measurements we mean that the sampling functions $s_{i}$ can only take two values either $\{0,1\}$ or $\{-1,1\}$. Without loss of generality we can assume that the model uses $\{-1,1\}$, as one can, by adding one extra measurement with the constant function, convert from the $\{0,1\}$ setup to the $\{-1,1\}$ model.

Binary measurements are a mainstay in signal and image processing due to the ”on-off” nature of many physical sampling devices. Microscopy is an obvious application as well as the newly emerging techniques of lensless cameras. In the discrete setting binary measurements are often modelled with Hadamard matrices, and this is one of the reasons why Hadamard matrices are so important in signal processing. To model binary measurements we change the model from Fourier samples $\langle f,s_{i}\rangle$, where the $s_{i}$ are complex exponentials to letting the $s_{i}$ be Walsh functions. The Walsh functions are the binary counterpart to Fourier samples and complex exponentials. Thus, the key question is as follows: what is the stable sampling rate when sampling with Walsh functions and reconstructing with wavelets? The answer is that it is linear regardless of the dimension when we consider separable boundary wavelets. This means that sampling with binary measurements is as efficient (up to potentially a different constant) as sampling with Fourier samples when reconstructing with wavelets. We expect the techniques used in this paper to extend to other X-lets as well; however, the extension, as in the Fourier case, is non-trivial.

1.2. Main Theorem

We consider the sampling space $\mathcal{S}$ of Walsh functions, which will be described in more detail in Chapter 2, and let the reconstruction space $\mathcal{R}$ be the space of boundary-corrected Daubechies wavelets (see Chapter 3 for details). The main theorem states that the stable sampling rate is indeed linear in $N$.

Theorem 1.1**.**

Let $\mathcal{S}$ and $\mathcal{R}$ be the sampling and reconstruction space spanned by the d-dimensional Walsh functions and separable boundary wavelets respectively. Moreover, let $N=2^{dR}$ with $R\in\mathbb{N}$. Then for all $\theta\in(1,\infty)$ there exists $S_{\theta}$ such that for all $M\geq 2^{dR}S_{\theta}$ we have $\mu(\mathcal{R}_{N},\mathcal{S}_{M})\leq\theta$. In particular, one gets $\Theta\leq S_{\theta}N$. Hence, the relation $\Theta(N;\theta)=\mathcal{O}(N)$ holds for all $\theta\in(1,\infty)$.

2. Walsh functions - Defining the Sampling Space $\mathcal{S}_{M}$

Due to the fact that we are dealing with the $d$-dimensional case, we introduce multi-indices to make the notation more readable. Let $j=\left\{j_{1},\ldots,j_{d}\right\}\in\mathbb{N}^{d}$, $d\in\mathbb{N}$ be a multi-index. A natural number $n$ is in the context of a multi-index interpreted as a multi-index with the same entry, i.e. $n=\left\{n,\ldots,n\right\}$. Then we define the addition of two multi-indices for $j,r\in\mathbb{N}^{d}$ by the pointwise addition, i.e. $j+r=\left\{j_{1}+r_{1},\ldots,j_{d}+r_{d}\right\}$ and the sum

[TABLE]

where $k,r\in\mathbb{N}^{d}$. The multiplication of an multi-index with a real number is understood pointwise, as well as the division by a multi-index. The absolute value of a multi-index $j$ is given by $|j|=j_{1}+\ldots+j_{d}$. The $d$ dimensional functions that we use in this paper are constructed by the tensor product. For a function $f:\mathbb{R}\rightarrow\mathbb{R}$ and an input parameter $\left\{x_{i}\right\}_{i=1,\ldots,d}=x\in\mathbb{R}^{d}$ with $x_{i}\in\mathbb{R}$ we use the following notation to present the $d$-times tensor product of $f$, i.e.

[TABLE]

It should be clear from the input parameter, whether $f$ represents the function on $\mathbb{R}$ or $\mathbb{R}^{d}$.

2.1. Defining Walsh functions

The key property that makes Walsh functions attractive in many applications is that they take only the values $1$ and $-1$. However, as Walsh functions are defined in the dyadic analysis, some properties only hold for dyadic addition. Recalling the basics of dyadic addition, we represent elements $x\in\mathbb{R}_{+}$ with their dyadic representation as follows

[TABLE]

where $x_{i}\in\left\{0,1\right\}$ for all $i\in\mathbb{Z}$. The natural extension always ends in [math] for dyadic rational numbers and is infinite for dyadic irrational numbers. The representation is therefore unique. Elements of $\mathbb{R}_{-}$ are represented as in the decimal analysis with an additional $-$ in front of the representation. In the dyadic analysis the addition $\oplus:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$ is defined by

[TABLE]

where $x_{i}\oplus_{2}y_{i}$ is addition modulo two, i.e. $0\oplus_{2}0=0,0\oplus_{2}1=1,1\oplus_{2}0=1,1\oplus_{2}1=0$. For negative numbers one has $-x\oplus y=x\oplus-y=-(x\oplus y)$. Note that there is no closed form between the decimal and dyadic addition. In particular, for two numbers $x,h\in\mathbb{R}$ the expression of the decimal sum $x+h$ has a different expression in the dyadic addition for every pair of numbers. This leads to further investigation during the proof of the main theorem. One part on the way to control the subspace angle is not to deal with all timeshifts of the wavelet but instead transfer them to the Walsh function and deal with the Walsh polynomial. Unfortunately, all properties of Walsh functions rely on dyadic rather than decimal addition. Therefore, this difference of the additions needs special care. In chapter 4.1 we will see that under mild assumptions the additions can made be equal. In chapter 4.2 we tweak the wavelets to match these assumptions.

Now, we present the Walsh functions, which are used to represent the sampling space $\mathcal{S}_{M}$. Therefore, we give a definition of the classical Walsh functions that highlights the difference between the possible orderings.

Definition 2.1** ([27]).**

Let $s\in\mathbb{N}$ and $x\in[0,1)$. Then there exists a unique $n=n(s)\in\mathbb{N}$ such that $s=\sum_{i=0}^{n-1}s_{i}2^{i}$, in particular $s_{n-1}\neq 0$ and $s_{k}=0$ for all $k\geq n$. Let $s^{n}=\left\{s_{0},\ldots,s_{n-1}\right\}$ and for $x=\sum_{i=-\infty}^{-1}x_{i}2^{i}$ define $x^{n}=\left\{x_{-n},\ldots,x_{-1}\right\}$, and $\omega_{W}$: $\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ by

[TABLE]

The Walsh functions are then given by

[TABLE]

By changing the matrix $\omega_{W}$ one gets different orderings of the Walsh functions. For example, the identity matrix leads to the Walsh-Kronecker functions, which have the drawback that with a change of $n(s)$ all functions are altered, hence one has to fix the maximal $s$ in advance. The Walsh-Paley ordering is obtained by replacing $\omega_{W}$ with the reversal matrix, i.e.

[TABLE]

They overcome the previous problem, but the functions are not ordered such that the number of zero crossings increases with $s$. Both drawbacks are overcome with the Walsh-Kaczmarz ordering presented in the previous definition.

The classical Walsh functions can be extended to the generalized Walsh functions $\operatorname{Wal}:\mathbb{R}^{2}_{+}\rightarrow\left\{-1,1\right\}$ which are defined with the classical Walsh functions and the periodic continuation with period $1$ by

[TABLE]

where $s$ and $x$ have the dyadic representation $(s_{i})_{i\in\mathbb{Z}}$ and $(x_{i})_{i\in\mathbb{Z}}$ and $s_{0}$, $x_{0}$ are the corresponding elements of the sequence. This extension can also be defined by letting $\omega_{W}$ be infinite, i.e. be defined over $\mathbb{Z}$ instead of $\mathbb{N}$ and hence allow inputs with infinite dyadic representations over $\mathbb{Z}$. Moreover, the Walsh functions can also be extended to negative inputs. Therefore, we define the following equality as in [27]

[TABLE]

The Walsh functions in higher dimensions are obtained by the tensor product, i.e. for $s=\left\{s_{k}\right\}_{k=1,\ldots,d},x=\left\{x_{k}\right\}_{k=1,\ldots,d}\in\mathbb{R}^{d}$

[TABLE]

The Walsh functions can also be combined to Walsh polynomials similar to trigonometric polynomials.

Definition 2.2**.**

Let $A,B\in\mathbb{Z}^{d}$ such that $A_{i}\leq B_{i},i=1,\ldots,d$ and $\alpha_{j_{i}}\in\mathbb{R}$. Then for $z\in\mathbb{R}^{d}_{+}$ we define the Walsh polynomial of order $n=|B|$ by $\Phi(z)=\sum_{j=A}^{B}\alpha_{j}\operatorname{Wal}(j,z)$. The set of all Walsh polynomials up to degree $n$ is given by

[TABLE]

With the generalized Walsh functions one can define a continuous and discrete transform. To ensure that the following integral exists, let $f\in L^{2}([0,1]^{d})$, the generalized Walsh transform is given almost everywhere by

[TABLE]

This is suitable for our setting, because we consider only the Walsh transform of functions that are supported in $[0,1]^{d}$. In the discrete setting we have for $N=2^{n}$, $n\in\mathbb{N}$ and $x=\left\{x_{0},\ldots,x_{N-1}\right\}\in\mathbb{R}^{N}$ that the one dimensional discrete Walsh transform of $x$ is given by $X=\left\{X_{0},\ldots,X_{N-1}\right\}$ with

[TABLE]

This transform corresponds, as mentioned, to the multiplication with a Hadamard matrix. By the definition of $\operatorname{Wal}$ it corresponds to the Hadamard matrix in Walsh-Kaczmarz ordering. The addition here is again the dyadic addition. The discrete d-dimensional Walsh transformed of $x\in\mathbb{R}^{N_{1}\times\ldots\times N_{d}}$ where $x_{k_{i}}\in\mathbb{R}$, $k=\left\{k_{i}\right\}_{i=1,\ldots,d},k_{i}=0,\ldots,N_{i}-1$ is given by $X=\left\{X_{j}\right\}\in\mathbb{R}^{N_{1}\times\ldots\times N_{d}}$, where $X_{j_{i}}\in\mathbb{R}$, $j=\left\{j_{i}\right\}_{i=1,\ldots,d},j_{i}=0,\ldots,N_{i}-1$, with

[TABLE]

2.2. Properties of Walsh functions

The Walsh functions obey the following properties: They are symmetric,

[TABLE]

and they obey the scaling property as well as the multiplicative identity, i.e

[TABLE]

and

[TABLE]

These properties can be directly transferred to the $d$-dimensional Walsh functions and the continuous Walsh transform, i.e. it holds, that the continuous Walsh transform is linear

[TABLE]

and obeys the following shift and scaling property, i.e.

[TABLE]

and

[TABLE]

3. Wavelets - Defining the Reconstruction Space $\mathcal{R}_{N}$

3.1. Boundary Wavelets

3.1.1. Boundary Wavelet space in one dimension

Daubechies boundary wavelets are deduced from general Daubechies wavelets. They have the advantage that they keep desirable properties, such as smoothness and vanishing moments, from their mother wavelet. In contrast, other approaches to find an orthonormal wavelet basis for $L^{2}([0,1])$, such as extension with zero, periodising or folding, loose smoothness. For the construction of the Daubechies boundary wavelets, as presented in [18], one starts with the Daubechies scaling functions. First, we deal with the scaling functions on the positive line $[0,\infty)$. Remember that a Daubechies scaling function $\phi$ of order $p$ has the support $[-p+1,p]$. Then the functions $\phi_{n}(x)=\phi(x-n)$ have their support completely in $[0,\infty)$ for $n\geq p-1$. However, these functions do not even generate the polynomials on $[0,\infty)$, so they do not represent smooth functions well. Therefore, the following functions are added to circumvent this issue:

[TABLE]

It is shown in [18] that these functions together with the inner ones, i.e. the translates of the scaling function, whose support is completely contained in the positive real line, span all polynomials with degree smaller or equal to $p-1$ on $[0,\infty)$. The same construction can be done for the negative line $(-\infty,0]$ and then be shifted by $1$ to get to the desired interval. This means in detail that the scaling function on the right hand side can be deduced from those on the left side, i.e. the construction for the right hand side results from a shift of $1$ in the functions that do intersect with the right end of the interval and a reflection. We have that

[TABLE]

In the next step we bring both systems together on $[0,1]$. To make sure that each shift of the scaling function is either an inner a left or a right scaling function, we consider scaling functions at a level $j\geq J_{0}$, where $2^{J_{0}}\geq 2p-1$. This way the support size of the scaling function at that scale is smaller than $1$. Therefore, the scaling function can intersect only with [math] or $1$ and hence the correction is well defined. The functions are now all corrected on the boundaries and they span the desired space $L^{2}([0,1])$. To form an orthonormal basis, we simply apply a Gram-Schmidt procedure. The new functions, after the orthonormalisation, are denoted by $\phi_{n}^{\text{left}},\phi_{n}^{\text{right}}$. The functions have staggered support, i.e. $\operatorname{supp}\phi_{n}^{\text{left}}=[0,p+n]$. Therefore, all $\phi$ have support length at most $2p-1$. Hence, the change to the boundary wavelet preserves the favourable property of a small support size. The dilated boundary scaling functions can be deduced from this construction, like the scaling functions for the real line. With this construction we obtain $2^{j}+2$ scaling functions at every scale $j$, but in many applications one prefers to have $2^{j}$ scaling functions. Therefore, we remove the two outermost interior scaling functions, i.e. those with the support closest to [math] and $1$ but not intersecting with them. This results in the subspaces

[TABLE]

where

[TABLE]

In [18] it is proven that we can define the wavelet space at every scale $j$ similar to the case on the real line by

[TABLE]

The original wavelet functions $\psi_{j,k}$ from the real line are in $W_{j}^{b}$ for $k=p,\ldots,2^{j}-p-1$. Because of the dimension of the scaling space one knows that $\dim W_{j}^{b}=2^{j}$. Hence, one has to add $2p$ wavelets. As they do not play an important role for the investigation of the main theorem we point the interested reader to [18] for detailed information.

Now that the necessary information about wavelets is introduced, we discuss the reconstruction space. The data is usually sparsely represented in the wavelet scheme, i.e. they only have large coefficients up to a certain scale. Therefore, the reconstruction space contains only the wavelets up to some scale $R$. Moreover, the low frequency part can be represented by the scaling space at some level $J\geq J_{0}$. In theory the choice of the lowest level $J$ is free. Nevertheless, it is common to use $J=J_{0}$. This results in the following reconstruction space. For $R\in\mathbb{N}$, the space of wavelets up to a scaling of $R$ is given by

[TABLE]

and has $N=2^{R}$ elements. Due to the construction the ”left” scaling functions are translates of the mother scaling functions and the ”right” scaling functions are reflected translated scaling functions, denoted by $\phi^{\#}$. Therefore,

[TABLE]

and every $\varphi\in\mathcal{R}_{N}$ with $||\varphi||=1$ has the representation

[TABLE]

3.1.2. Boundary wavelets in higher dimensions

In this paper, we also consider the $d$-dimensional case. For the reconstruction in $d$-dimensions we focus on separable boundary wavelets. Therefore, the $d$-dimensional wavelets can be derived from the one dimensional case by tensoring the scaling space and then studying the according wavelet space.

From (3.4) we have the one dimensional boundary scaling function. With the tensor product we get the $d$ dimensional one, i.e. $\phi_{j,n}^{d}=\phi_{j,n_{1}}^{b}\otimes\ldots\otimes\phi_{j,n_{d}}^{b}$ for $n=\left\{n_{1},\ldots,n_{d}\right\}\in\mathbb{N}^{d}$, $j\geq J_{0}$. To make this easier to read we set $\phi_{j,n}:=\phi_{j,n}^{d}$, as the dimension is defined by the context. Then the $d$-dimensional scaling space is given by

[TABLE]

For the purpose of constructing higher dimensional boundary wavelets we exploit the MRA structure. We have that

[TABLE]

Therefore, we can divide in higher dimensions the scaling space at one level in the scaling space and the wavelet space in the lower level

[TABLE]

This way we have defined the $d$ dimensional boundary-corrected wavelet space $W_{j-1}^{b,d}$ by

[TABLE]

Due to (3.6) we only have to focus on the scaling space, as the sum over the wavelet spaces can be represented by the scaling space at highest scale. Therefore, we do not explain details about the wavelets here. We have with (3.7)

[TABLE]

According to the size of the one dimensional scaling space, we know that the $d$-dimensional scaling space has size $2^{dj}$. The reconstruction space for $N=2^{dR}$ is then

[TABLE]

In order to get the scaling space for the boundary wavelets in one dimension we had to reflect the scaling function for translates $k=2^{j}-p,\ldots,2^{j}-1$. This means that $V_{j}^{b,d}$ is spanned by the translates of $2^{d}$ functions, which are constructed by tensoring of the original scaling function $\phi^{0}:=\phi$ and the translated version $\phi^{1}:=\phi^{\#}$. Define $K_{0}=\left\{0,\ldots,2^{j}-p-1\right\}$ and $K_{1}=\left\{2^{j}-p,\ldots,2^{j}-1\right\}$, then the mapping $m:\left\{0,\ldots,2^{j}-1\right\}\mapsto\left\{0,1\right\}$ is given by

[TABLE]

This allows us to represent $\varphi\in V_{j}^{b,d}$ with $||\varphi||=1$ by

[TABLE]

where

[TABLE]

So at this point we have the wavelets and the scaling functions which span the space $L^{2}([0,1]^{d})$ and therefore the reconstruction from Walsh functions in the wavelet space is guaranteed.

4. The main theorems and its proof

Taken together, we can now prove the main result.

4.1. Useful lemmas about Walsh functions

For the proof of theorem 1.1 we have to combine the properties of Walsh functions in the dyadic analysis with the properties of the wavelets in the decimal analysis. Therefore, we consider under which conditions the decimal and dyadic additions are equal. This is important to combine the multiplicative identity of the Walsh functions with the translates of the wavelets.

Lemma 4.1**.**

Let $x\in[0,1)$ and $m\in\mathbb{N}$ then $x+m=x\oplus m.$

Proof.

The dyadic representation of $x$ is $\left\{\ldots,0,x_{-1},x_{-2},\ldots\right\}$ and the dyadic representation of

[TABLE]

Because the representations do not have non-zero elements at the same position, one achieves the following

[TABLE]

∎

Next, we look at the inverse element for the dyadic addition. This is also discussed in [28] and will be used in Corollary 4.3.

Lemma 4.2**.**

The dyadic sum of two numbers $x,y\in\mathbb{R}_{+}$ is [math] if and only if $x=y$.

Proof.

Let $x,y\in\mathbb{R}_{+}$ with the dyadic representation $\left\{x_{i}\right\}_{i\in\mathbb{Z}}$ and $\left\{y_{i}\right\}_{i\in\mathbb{Z}}$. Then

[TABLE]

if and only if $x_{i}\oplus_{2}y_{i}=0$ for all $i\in\mathbb{Z}$. This is the case if and only if $x_{i}=y_{i}$ for all $i\in\mathbb{Z}$, i.e. $x=y$. ∎

With this the relation between the decimal addition and the multiplicative identity of the Walsh functions can be found.

Corollary 4.3**.**

Let $t\in\mathbb{N}$ and $x\in[0,1)$, then the following holds:

[TABLE]

Proof.

With Lemma 4.1 we have that $x\oplus t=x+t$. This allows

[TABLE]

∎

Next, we analyse the sum of Walsh functions with equally distributed inputs. This will be used in Lemma 4.5.

Lemma 4.4**.**

Let $N=2^{n}$, $n\in\mathbb{N}$, then for all $s\in\mathbb{N}_{0}$ the following property holds:

[TABLE]

Proof.

The first case for $s=0$ follows directly by the definition of the Walsh function as $\operatorname{wal}(0;x)\equiv 1$ for all $x\in[0,1)$. For the second part we use the equal distribution of Walsh functions in intervals where it takes the values $-1$ and $1$, i.e. for $s\leq 2^{m}$, $m\in\mathbb{N}_{0}$ the Walsh function $\operatorname{wal}(s;x)$ takes the value $1$ on $2^{m-1}$ intervals of length $\frac{1}{2^{m}}$ and $-1$ on the same number of intervals of that length [27]. As the sequence $\left\{i/N\right\}_{i=0,\ldots N-1}$ is equally distributed on this interval, the sum equals [math]. ∎

With this information in hand we can now prove the following lemma, which shows a relation between the values of the discrete Walsh transform and the signal itself. This will then be used in Lemma 4.6.

Lemma 4.5**.**

Let $N=\left\{N_{i}\right\}_{i=1,\ldots,d}$, where $N_{i}=2^{n_{i}},n_{i}\in\mathbb{N}$ and $i=1,\ldots,d$. Let $x\in\mathbb{R}^{N_{1}\times\ldots\times N_{d}}$, where $x=\left\{x_{k}\right\}$ and $x_{k_{i}}\in\mathbb{R}$, $k=\left\{k_{i}\right\}_{i=1,\ldots,d},k_{i}=0,\ldots,N_{i}-1$ be a discrete d-dimensional signal. Given the discrete d-dimensional Walsh transformed by $X=\left\{X_{j}\right\}\in\mathbb{R}^{N_{1}\times\ldots\times N_{d}}$, where $X_{j_{i}}\in\mathbb{R}$, $j=\left\{j_{i}\right\}_{i=1,\ldots,d},j_{i}=0,\ldots,N_{i}-1$, with

[TABLE]

it follows that

[TABLE]

Proof.

First, one observes that by definition and the fact that $N_{i}=2^{n_{i}},i=1,\ldots,d$ the following holds:

[TABLE]

Next, recalling (4.5) and Lemma 4.2 we directly get the desired property by

[TABLE]

∎

With this we can achieve a useful relation between the Walsh polynomial, i.e. $\Phi(z)=\sum_{j=A}^{B}\alpha_{j}\operatorname{Wal}(j,z)$ with $A,B\in\mathbb{Z}^{d}$ and $\alpha_{j_{i}}\in\mathbb{R}$ for all $j_{i}=A_{i},\ldots,B_{i}$, $i=1,\ldots,d$, and its coefficients $\alpha$ similar to the trigonometric polynomial in [3].

Lemma 4.6**.**

Let $A,B\in\mathbb{Z}^{d}$ such that $A_{i}\leq B_{i},i=1,\ldots,d$ and consider the Walsh polynomial $\Phi(z)=\sum_{j=A}^{B}\alpha_{j}\operatorname{Wal}(j,z)$ for $z\in\mathbb{R}^{d}_{+}$. If $L=\left\{L_{1},\ldots,L_{d}\right\}$ with $L_{i}=2^{n_{i}}$, $n_{i}\in\mathbb{N},i=1,\ldots,d$ such that $2L_{i}\geq B_{i}-A_{i}+1$, then

[TABLE]

Proof.

For the proof Lemma 4.5 is used. Therefore, let $x=\left\{x_{k}\right\}$, where $x_{k_{i}}\in\mathbb{R}$, $k=\left\{k_{i}\right\}_{i=1,\ldots,d},k_{i}=0,\ldots,2L_{i}-1$ and $x\in\mathbb{R}^{2L_{1}\times\ldots\times 2L_{d}}$ with the discrete Walsh transformed $X=\left\{X_{j}\right\}\in\mathbb{R}^{2L_{1}\times\ldots\times 2L_{d}}$, where $X_{j_{i}}\in\mathbb{R}$, $j=\left\{j_{i}\right\}_{i=1,\ldots,d},j_{i}=0,\ldots,2L_{i}-1$ and $X\in\mathbb{R}^{2L_{1}\times\ldots\times 2L_{d}}$. Consider the sums $k_{i}+A_{i}$, $i=1,\ldots,d$, there exists a number $\tilde{A}_{i}(k)$, such that $k_{i}+A_{i}=k_{i}\oplus\tilde{A}_{i}(k_{i})$ for all $i=1,\ldots,d$. As before, we denote by $\tilde{A}(k)$ the multi-index containing all $\tilde{A}_{i}(k_{i})$. Define the coefficients

[TABLE]

and the sequence $x$ as follows:

[TABLE]

Then one gets with the scaling property (2.16) and the multiplicative identity (2.17)

[TABLE]

With that one can conclude

[TABLE]

∎

4.2. Changes of Wavelets

As mentioned in chapter 2 the decimal and dyadic addition do not correspond directly to each other. In particular, the representation of the decimal addition with a number $h$ to a number $x$ depends on both parts of the sum. However, in Corollary 4.3 we have seen that for $x\in[0,1)$ and $n\in\mathbb{N}$ the dyadic and decimal addition coincide. In the proof of the main theorem we want to transfer the time shifts of the wavelet to the Walsh function, i.e.

[TABLE]

Therefore, to enable us to use Corollary 4.3 and make the last equation an equality, the domain of the wavelets needs to be restricted to $[0,1]^{d}$. This is not a contradiction to the construction of the previous chapter, because the functions $\phi_{R,n}$ are indeed supported in $[0,1]^{d}$. However, the scaling function at level [math] is not only supported there and that is the function that we are dealing with after the change of variables in the integral. To solve this problem we represent the scaling function as a sum of functions that are supported in $[0,1]$, i.e.

[TABLE]

and

[TABLE]

This can also be done accordingly for the reflected function $\phi^{\#}$. In the higher dimensional case we have

[TABLE]

and $\phi_{i_{k}}$ defined as above. This way the multiplicative identity holds also for the decimal time shift of the wavelets.

4.3. Useful lemma about wavelets

For the proof of the main theorem we have to bound the decay rate of the Wavelets under the Walsh transform. We analyse the Haar wavelets and the boundary-corrected Daubechies wavelet separately.

For the Haar wavelet we have that the scaling function $\phi=\mathcal{X}_{[0,1]}$. Hence,

[TABLE]

Therefore we have

[TABLE]

Next, we use the result from [9] regarding the decay rate of Hölder continuous functions. We have that the $\phi_{i}$ and $\phi_{i}^{\#}$ have the same smoothness properties as the original function $\phi$. Hence, we have for the Daubechies wavelet of order $2$ that it is Hölder continuous with coefficient $\alpha=0.55$. The higher order wavelets are all Hölder continuous with coefficient $\alpha=1$. The results in [9] hold only for inputs in $\mathbb{N}$. We will use the same proof strategy and notation for inputs in the form of $\frac{k}{2^{R}}$ for some $k,R\in\mathbb{N}$.

Lemma 4.7**.**

Let $\phi$ be a Daubechies wavelet of order $p>1$. And let $m\geq 1$,$R\in\mathbb{N}$, $L=2^{R}$ and $j=0,\ldots,L-1$. Then we have that

[TABLE]

and

[TABLE]

for some constants $C$ and $C^{\#}$.

Proof.

First, observe that

[TABLE]

We want to use the same technique as in [9] to divide the integral into parts where the Walsh function takes the values $+1$ and $-1$ on intervals of the same length. Let $q\in\mathbb{N}$ such that $2^{q}\leq m<2^{q+1}$, then we have that $2^{q+R}\leq j+Lm<2^{q+R+1}$. Now, define the intervals $\Delta_{k}^{q}=\left[2^{-q}k,2^{-q}(k+1)\right)$. It was shown that the function $\operatorname{Wal}(j+Lm,s)$ is constant on the interval $\Delta_{2l}^{q+R+1}$ and $\Delta_{2l+1}^{q+R+1}$ and takes the values $+1$ on one of them and $-1$ on the other. Hence, $\operatorname{Wal}(j+Lm,\frac{x}{L})$ takes these values on the intervals $\Delta_{2k}^{q+1}$ and $\Delta_{2k+1}^{q+1}$. Due to the Hölder continuity we have that there exists a constant $C$ such that $\phi_{i}(x)\leq\phi_{i}(s)+C|x-s|^{\alpha}$ for all $s\in[0,1]$. With this one gets

[TABLE]

Hence,

[TABLE]

Thus, we get for the complete integral on $[0,1]$

[TABLE]

The proof holds analogously for $\phi_{i}^{\#}$. ∎

4.4. Proof of the main theorem

With the tools established above, we can now prove the main result. To make the exposition easier to read we first prove the theorem in one dimension and then make the generalisation to several dimensions in a separate proof. Given the setup with the multi-indices framework, this can be done reasonably smoothly.

Proof of Theorem 1.1 in one dimension.

The aim of this proof is to find for every $\theta\in(1,\infty)$ an integer $S_{\theta}\in\mathbb{N}$, such that for all $M\geq S_{\theta}N$ the subspace angle is bounded, i.e. $\mu(\mathcal{R}_{N},\mathcal{S}_{M})\leq\theta$. Let $R\in\mathbb{N}$ be the number of reconstructed levels, i.e. $N=2^{R}$. We start with a suitable representation of $\cos(\omega(\mathcal{R}_{N},\mathcal{S}_{M}))$. There exists $\varphi\in\mathcal{R}_{N}$ with $||\varphi||=1$ such that $\inf_{f\in\mathcal{R}_{N},||f||=1}||P_{\mathcal{S}_{M}}f||=||P_{\mathcal{S}_{M}}\varphi||$, because the closed unit ball in $\mathcal{R}_{N}$ is compact and $P_{\mathcal{S}_{M}}$ is continuous. By (3.8) we can represent $\varphi$ as

[TABLE]

and

[TABLE]

The first equation (4.38) allows us to deal only with the scaling function instead of both the wavelets at different scales and the scaling function. The second one (4.39) enables us to bound $P_{\mathcal{S}_{M}}^{\perp}\varphi$ from above in lieu of $P_{\mathcal{S}_{M}}\varphi$ from below.

Instead of dealing with all different shifts of the scaling function, we aim for a closed form that only depends on the functions $\phi$ and $\phi^{\#}$. An essential part in the construction of this is the use of the scaling property in Corollary 4.3. Therefore, it is necessary, that $m\in\mathbb{N}$ and $x\in[0,1)$. For this sake, the functions $\phi_{i}$ were defined in (4.22) and we define

[TABLE]

with $z\mapsto p_{R}(z)$ and $p_{R}(z)$ being the smallest integer such that $p_{R}(z)2^{R}+z>0$. This yields

[TABLE]

where we used in the last line the fact that the Walsh functions are $1$-periodic, if the other input data is an integer. Then, we have that $x\in[0,1]$ and $n+i-1+2^{R}p_{R}(n+i-1)\in\mathbb{N}$. Hence, Corollary 4.3 can be used in the third line to get

[TABLE]

With

[TABLE]

it results in

[TABLE]

Analogously this can be done for the reflected function $\phi^{\#}$. Thus, by using

[TABLE]

similar to the above, we get

[TABLE]

This representation is very useful. Indeed, one only has to pay attention to the decay rate of the Walsh transform of the pieces of the mother scaling function, its reflection and the Walsh polynomial. Moreover, the pieces fulfil $f(t)=0$ for $t<0$, such that the Walsh transform with one kernel can be used for the analysis. Consider

[TABLE]

Next, one uses the linearity of the orthogonal projection to change the order of the summands, such that the sum over the scaling function pieces can be dealt with in the end. First, we take out the sum over the parts of the wavelet $\phi_{i,R,n}$ to handle every cut-out of the wavelet separately. In particular, by (4.56)

[TABLE]

Second, it follows, by using (4.51), (4.54) and the Cauchy-Schwarz inequality, that

[TABLE]

We will only deal with the first summand and the other follow analogously. In the following step the $1$- periodicity of the Walsh function is used. Let $S\in\mathbb{R}_{+}$ such that $M=S2^{R}$, then by replacing $k\geq M$ with $k=mL+j$ with $L=2^{R}$, $m\geq S$ and $j=0,\ldots,L-1$ we have

[TABLE]

The last sum can be estimated via Lemma 4.7 and (4.27) by

[TABLE]

with $\alpha=1$ for the Haar wavelet and $\alpha$ being the Hölder coefficient for the other Daubechies wavelets of order $p>1$. For the first sum $\sum\limits_{j=0}^{L-1}\frac{1}{L}\left|\Phi_{i}(\frac{j}{L})\right|^{2}$ we have with $\alpha_{n}=0$ for $n=-p+1,\ldots,-1$ that

[TABLE]

and $L=2^{R}$, such that we can use Lemma 4.6 as $2^{R}-p+i-1+2^{R}p_{R}(i-1)-(-p+1+i-1+2^{R}p_{R}(i-1))+1=2^{R}-p-(-p+1)+1=2^{R}\geq L$ and obtain

[TABLE]

Altogether this gives with Lemma 4.7

[TABLE]

and similarly

[TABLE]

Using (4.69), (4.70) and (4.71) yields the following estimation

[TABLE]

Thus, $||P_{\mathcal{S}_{M}}^{\perp}\varphi||\leq\gamma$ whenever

[TABLE]

where $C=4\max\left\{A^{2},{A^{\#}}^{2}\right\}$. It follows from (4.39) that $\cos(\omega(\mathcal{R}_{N},\mathcal{S}_{M}))\geq 1-\gamma\geq\frac{1}{\theta}$ ,i.e. $\mu(\mathcal{R}_{N},\mathcal{S}_{M})\leq\theta$, whenever the constant $S$, which is dependent on $\theta$ and therefore denoted by $S_{\theta}$, fulfils (4.74) with $\gamma=1-\frac{1}{\theta}$, i.e

[TABLE]

and $M=S_{\theta}L=S_{\theta}N$. ∎

Proof of Theorem 1.1 in $d$-dimensions.

In higher dimensions we represent $\varphi\in\mathcal{R}_{N}$ in terms of the sum over $2^{d}$ different tensor products. Then, we need to investigate the inner products of these summands with the Walsh functions as in the one dimensional case. At this point the results from the one dimensional case come into play. Next, one investigates the parts of the set $I_{M}^{\perp}$, where $I_{M}=\left\{l=\left\{l_{1},\ldots,l_{d}\right\},l_{k}=0,\ldots,M_{k}-1\right\}$ and $M=\left\{M_{1},\ldots,M_{d}\right\}\in\mathbb{N}^{d}$ and $I_{M}^{\perp}=\mathbb{N}^{d}\setminus I_{M}$, which correspond to the largest estimates of the inner products of the summands of the wavelet and the Walsh function. Finally, these can be bounded with estimates from the one dimensional case and additional care for the finite sums.

Now, we present the described steps in more detail. Let $\varphi\in\mathcal{R}_{N}$ with $||\varphi||=1$. Then we can represent $\varphi$ as in (3.16) in the following sum

[TABLE]

In the one dimensional case we derived the representation of $\sum_{n=0}^{2^{R}-p-1}\alpha_{n}\langle\phi_{i,R,n},\operatorname{Wal}(k,\cdot)\rangle$ in terms of the Walsh transform of the wavelet and the Walsh polynomial, i.e. 2^{-R/2}\stackon[2pt]{\phi_{i}}{\stretchto{\scalerel*[\widthof{\phi_{i}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptscriptstyle W}}}}(\frac{k}{2^{R}})\Phi_{i}(\frac{k}{2^{R}}). This equality from (4.51) is used to represent the inner product in higher dimensions. For this we need to define $p_{R}$ from (4.41) for higher dimensions. In particular, let $p_{R}:\mathbb{Z}^{d}\rightarrow\mathbb{N}^{d}$ with $\left\{z_{i}\right\}_{i=1,\ldots,d}=z\mapsto p_{R}(z)=\left\{p_{R}(z)_{i}\right\}_{i=1,\ldots,d}$ and $p_{R}(z)_{i}$ being the smallest integer such that $p_{R}(z)_{i}2^{R}-z_{i}>0$ for all $i=1,\ldots,d$. Further, let $l=(l_{1},\ldots,l_{d})\in\mathbb{Z}^{d}$. This yields

[TABLE]

Here, the problem is reduced to the one dimensional case and we can apply (4.42) and (4.46) to get

[TABLE]

This is defined now as in (4.50)

[TABLE]

Note that the different definitions from the one dimensional case for $\Phi$ and $\Phi^{\#}$ are combined in the notation with the $K_{s}$. We get with this the presentation of the inner products as desired:

[TABLE]

For the representation of indices that correspond to the sampling functions, let $I_{M}=\{l=\{l_{1},\ldots,l_{d}\},l_{k}=0,\ldots,M_{k}-1\}$, where $M=\{M_{1},\ldots,M_{k}\}\in\mathbb{N}^{d}$ is the number of samples. Then $l\notin I_{M}$ corresponds to $l>m$ in the one dimensional case. We now want to analyse the orthogonal projection on the orthogonal complement of the sampling space

[TABLE]

This way we ensured the use of Corollary 4.3. Next we change the order again to deal with the different cut-out functions separately. This was seen already in (4.58). We get

[TABLE]

With the Cauchy Schwarz inequality and a careful reordering we get as in (4.61)

[TABLE]

Now, let $S\in\mathbb{N}$ be given, such that the number of samples $M=\left\{M_{1},\ldots,M_{k}\right\}\in\mathbb{N}^{d}$ is $M_{k}=S2^{R}$. Then if $l=\left\{l_{1},\ldots,l_{d}\right\}\notin I_{\mathcal{M}}$ at least one $l_{k}>M_{k}$. The sum is the largest if only one $l_{k}$ fulfils this estimate. Hence, without loss of generality let $l_{1}>M_{1}$ and $l_{k}\leq M_{k}$ for $k=2,\ldots,d$. Now let $l_{k}=j_{k}+u_{k}2^{R}$ and $\alpha=0.55$ or $\alpha=1$ depending on the chosen wavelet. Then we get similar to (4.65):

[TABLE]

The last sum can be estimated with the help of Lemma 4.6. In (4.66) and (4.69) this was derived in the one dimensional case that can be directly used here, such that

[TABLE]

From the fact that the $\phi_{R,n}^{s}$ form an orthonormal bases and $||\varphi||=1$ we have

[TABLE]

This together with the fact that $|\left\{0,1\right\}^{d}|=2^{d}$ gives

[TABLE]

By replacing $C^{d-1}$ with $C$ we have

[TABLE]

Thus, $||P_{\mathcal{S}_{M}}^{\perp}\varphi||\leq\gamma$ whenever

[TABLE]

It follows from (4.39) that

[TABLE]

i.e. $\mu(\mathcal{R}_{N},\mathcal{S}_{M})\leq\theta$, whenever $S_{\theta}$ fulfils (4.93) with $\gamma=1-\frac{1}{\theta}$, i.e.

[TABLE]

and $|M|=M_{1}\cdot\ldots\cdot M_{d}=S_{\theta}^{d}2^{dR}=S_{\theta}^{d}N$. ∎

5. Numerical Experiments

In this chapter we underline the theoretical results with numerical experiments. We first calculate the stable sampling rate for different stabilities $\theta$ and Daubechies wavelets. Then, we see that the reconstruction with generalized sampling leads to much better results then the direct inversion with the Walsh transform. Moreover, we point out that it is important to consider the stable sampling rate, as otherwise the reconstruction gets very unstable with meaningless results.

First, we see in Figure 1 the stable sampling rate for the different wavelets and stabilities in the one-dimensional case. One can see that it is indeed linear with jumps according to the levels of the wavelets. Moreover, it is easy to detect that the constant $S_{\theta}$ is considerably low such that the number of samples needed is only marginally larger than the number of coefficients that we reconstruct. It is not surprising that the stable sampling rate gets larger for smaller $\theta$. In the theory of the reconstruction from Fourier measurements we have a direct relation between the smoothness of the wavelets and the size of the stable sampling rate. Similar relations are not known for the Walsh wavelet case.

In Figure 2, we demonstrate the reconstruction with generalized sampling. For this sake we consider two different functions on $[0,1]$. First, we look at the cosine function in Figure 2a, taking $77$ Walsh samples. In Figure 2c the direct inversion is shown. It is clear that the reconstruction has a lot of block artefacts, whereas the reconstruction with generalized sampling of $64$ Daubechies $8$ wavelets has nearly no visible artefacts. The same artefacts can be seen in 2f. In this case $192$ Walsh samples were taken and $128$ wavelet coefficients were reconstructed. The artefacts with the direct Walsh inverse are much stronger than the common Gibbs phenomena for the Fourier case. Because of this, reconstructions with Walsh functions are not feasible in practice. They are also the reason why one does not use Haar wavelets as they obey the same block artefacts. This underlines the need of a reconstruction technique that refers the data from the sampling space to a much more appropriate reconstruction space, where the data is represented sparsely. In this case we get fewer artefacts. For completely continuous functions as in 2b the reconstruction with generalized sampling has nearly no artefacts. In case of discontinuities as in the second function 2d one gets some artefacts at the discontinuities as can be seen in 2e. Even so, the overall reconstruction quality is much better and the reconstruction still obeys the regularity properties of each part of the function.

Nevertheless, it is important to take the stable sampling rate in mind. If one tries to reconstruct with fewer samples then needed, the reconstruction gets very unstable and one gets meaningless results. This can be seen in Figure 3, where a function on $[0,1]$ is reconstructed from $512$, which is much more than the $77$ and $192$ for the other functions.

6. Conclusion

We were able to investigate a very important part of the error estimate for different reconstruction methods. Moreover, we showed that binary measurements modelled by Walsh functions are well suited to reconstruct images with wavelets. This, together with the results in [1, 3], gives a broad knowledge about the accuracy and stability for two major applications of sampling theory, i.e. systems with Fourier samples and those with binary measurements.

Acknowledgement

LT acknowledges support by the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L016516/1 for the University of Cambridge Centre for Doctoral Training, the Cambridge Centre for Analysis. AH acknowledges support from a Royal Society University Research Fellowship.