On the Unification of Schemes and Software for Wavelets on the Interval

We revisit the construction of wavelets on the interval with various degrees of polynomial exactness, and explain how existing schemes for orthogonal- and Spline wavelets can be extended to compactly supported delay-normalized wavelets. The contribution differs substantially from previous ones in how results are stated and deduced: linear algebra notation is exploited more heavily, and the use of sums and complicated index notation is reduced. This extended use of linear algebra eases translation to software, and a general open source implementation, which uses the deductions in this paper as a reference, has been developed. Key features of this implementation is its flexibility w.r.t. the length of the input, as well as its generality regarding the wavelet transform.


Introduction
Wavelets on the interval are well studied, with several existing constructions addressing various degrees of polynomial exactness at the primal and dual sides. Developments in this respect can be traced back to [7,17] (see also [1,15]). The most common construction of orthogonal wavelets on the interval is possibly [9], while the first constructions of biorthogonal spline wavelets on the interval stem from [11]. More recent constructions of such wavelets (see for instance [5,6,19]) aim at improving the condition of the bases.
Software utilizing the mentioned constructions for wavelets on the interval is limited, however. Most software involving wavelet transforms typically abandon polynomial exactness in favour of simpler extension strategies at the boundaries, such as periodic or sym-metric extensions [24]. 1 To the authors knowledge, there does not exist openly available software embracing the constructions mentioned above. 2 This paper is an attempt to establish the fundaments for software supporting recent constructions of wavelets on the interval, and is closely tied to an open source implementation. 3 In the process some of the constructions will be revisited, and their proofs will be rewritten (in terms of linear algebra and notation, see Sect. 1.1) so that they differ substantially from that found in existing literature, and so that translation to software is straightforward. When rewriting the proofs it will be seen that they also apply for more general compactly supported wavelets, not only to the orthogonal-or Spline cases mostly found in the literature: For most cases of delay-normalized wavelets (see Sect. 2), they also apply, as well as the well-known method of stable completion [4]. It will also be explained how the more recent construction in [19] can be put into this context of delay-normalized wavelets.
The reader will note that the presented deductions follow the same line as those in [11] (with substantial changes to the notation). As the bases therein suffer in terms of conditioning, this choice deserves some comments. The improved conditioning in recent constructions forces a primal multiresolution analysis (including the boundary functions, and their number) to be fixed. However, this puts constraints on the dimensions of the input to the corresponding Discrete Wavelet Transform: As an example, an m-level DWT on the interval as defined in [9] is possible only if the input length is divisible with 2 m . [11] offers flexibility in this respect, however, so that it is adaptable to the input length and the number of levels.
In essence this flexibility lies in absorbing some of the internal scaling functions in the constructed boundary functions, which also partially explains the resulting bad conditions. In summary, one either would like to transform a vector over a given number of levels 1. with no constraints on the length, for which the strategy from [11] can be applied (although bad conditioning may result), 2. with a highly constrained length, that meets the requirement in more recent contributions such as [19].
If one accepts the constraint mentioned in 2., the software implementation can easily support the construction in [19] as well.
A comment is also in order as to why one also should address wavelets on the interval outside the biorthogonal spline case. By far spline wavelets are the most common in the literature due to their favourable properties, and the large machinery available for them. In practice, however, there may be specific requirements on the wavelet transform in question, thereby excluding spline wavelets. The deductions given here then can be useful in such cases.
As a final note, the reader should be aware that the purpose of this paper is to serve as reference deductions for an implementation. As such the paper does not construct new wavelets, it does not address recent successful applications of wavelets to the numerical solutions of PDE's, and multiwavelets are not addressed. The software implementation hopefully can serve as a playground for researchers experimenting with wavelets, and for others to extend.

Notation
The paper follows notation in [21,22], which introduce the reader to signal processing and wavelets in a linear algebra friendly way, and in a style very different from that common for wavelets. The books also use a similar software implementation, and actually build it from scratch. The interval notation [a, b] = {a, a + 1, . . . , b} will be used to denote all integers between the two integers a and b. If b < a, [ These can also be combined, i.e., for k 1 , k 2 ∈ Z, one has k 1 + k 2 [a, b] = {k 1 + k 2 a, k 1 + k 2 (a + 1), . . . k 1 + k 2 b}.
This notation will be used to refer to segments of matrices and vectors. It should be clear from the context whether a range of integers, or an actual interval on the real line, is meant. This notation will eliminate much of the extensive indexing in wavelet literature. In particular it will simplify referring to segments of the DWT/IDWT matrices, as will often be needed. In the literature a DWT/IDWT is often expressed in terms of the filter coefficients, since these represent all entries in those matrices. This brings one away from simple matrixvector expressions, and our deductions will therefore avoid this. Wavelet bases for L 2 (R) contain an infinite number of basis functions at each resolution, whereas wavelet bases on the interval have finitely many. It will therefore be convenient to mix notation for finite and infinite matrices, and allowing finite matrices and vectors to have any given legal range of row-and column indices. In particular, in an expression on the form ⎛ it will be assumed that the column vector on the left hand side has row indices [0, N − 1], and that the column vector on the right hand side has row indices [−R + 1, N − 1]. The matrix C can be any infinite matrix, but when written as above it will be assumed that the range of column-and row indices in C are [0, N − 1] × [−R + 1, N − 1], i.e., that the indices match. Since any range of row-and column-indices may be legal, entries with index 0 or (0, 0) will occasionally be underlined (as in filter notation in signal processing), to make positions clear. The MATLAB notation that a simple colon denotes all elements along an axis, will also be followed.

Organization of the Paper
The paper is organized as follows. In Sect. 2 the general setup for wavelets is introduced, and in Sect. 3 the setup is specialized to the interval. In Sects. 4 and 5 the scaling functions and the corresponding mother wavelets are constructed. While those sections were adapted to the left end of the interval, Sect. 6 explains how delay-normalizedness ensures that the construction at the right end can be obtained from a simple mirroring operation of the left end. In Sect. 7 the result in [19] are put into the context of this contribution, and some notes on the software implementation can be found in Sect. 8. A more detailed explanation of this implementation can be found in the technical report [2].

Setup for Wavelets on the Entire Real Line
Let φ and ψ be the scaling function and the mother wavelet of a compactly supported wavelet. Assume also that φ is exact of order N (meaning that all polynomials of degree less than N can be written as linear combinations of the translates {φ(t − n)} n ). Similarly letφ,ψ , andÑ be the corresponding quantities for the dual wavelet. The resolution space V 0 is the space spanned by the translates φ 0,n (t) = φ(t − n), while the detail space W 0 is the space spanned by ψ 0,n (t) = ψ(t − n). For m > 0, the resolution-and detail spaces V m and W m are the spaces spanned by the dilated functions respectively. Similar definitions apply forφ andψ . One also writes so that V m = span(φ m ) and W m = span(ψ m ). When φ gives rise to a multiresolution analysis (i.e., where the dilated scaling functions and mother wavelets are listed in alternating order) is also a basis for V m . This alternating order of the basis functions is non-standard in wavelet literature, where all φ m−1,n -functions usually precede the ψ m−1,n . This reordering has the advantage that the index n into the basis C m represents time, and that change of coordinate matrices involving those bases will be banded.
On the dual side one similarly definesφ m ,ψ m ,Ṽ m ,W m , andC m . The Gramm matrix of two bases B = {b i } i and C = {c j } j , denoted ( B, C ), is the matrix with entries b i , c j . If ( B, C ) = I one also says that B and C are biorthogonal. A wavelet is called biorthogonal if the corresponding bases are biorthogonal, i.e., ( φ m ,φ m ) = ( C m ,C m ) = I . Some of the most used biorthogonal wavelets were established in [8]. Some of the most used orthonormal wavelets, for which φ =φ, ψ =ψ , and ( φ m , φ m ) = ( C m , C m ) = I (i.e., both φ m and C m are orthonormal bases for V m ) were established in [12]. Denoting by supp(f ) the support interval of the function f , a convention therein is that The change of coordinates from φ m to C m is called the (one-level) Discrete Wavelet Transform, or DWT, and denoted H (i.e., H = P Cm←φ m ). Its inverse is the IDWT, denoted by G (i.e., G = P φ m ←Cm ), and can be written as Since the bases here are doubly infinite, the component with index zero is emphasized by underlining it, i.e., the coordinate vector of f (t) = c −1 φ 0,−1 + c 0 φ 0,0 + c 1 φ 0,1 in φ 0 will be written as [f ] φ 0 = (c −1 , c 0 , c 1 ). Coefficients which are zero were not listed here, as is common in signal processing filter notation. H and G can be expressed in terms of filters as follows [21,22,Chap. 3]: 1. The even-indexed rows of H coincide with those of a (low-pass) filter matrix, denoted H 0 .
2. The odd-indexed rows of H coincide with those of a (high-pass) filter matrix, denotes H 1 .
3. The even-indexed columns of G coincide with those of a (low-pass) filter matrix, denoted G 0 . 4. The odd-indexed columns of G coincide with those of a (high-pass) filter matrix, denoted G 1 .
Thus, H can be alternatively defined as the unique matrix compatible with filters H 0 and H 1 , and G as the unique matrix compatible with filters G 0 and G 1 . It is known (Exercise 5.10 in [21,22]) that if the filters of a wavelet are finite impulse response (FIR), then there exist an integer d and α ∈ R so that Since the alternating sign corresponds to a shift in frequency by π , this says that, up to multiplication with a scalar, 1. H 1 is the high-pass filter corresponding to the low-pass filter G 0 , 2. G 1 is the high-pass filter corresponding to the low-pass filter H 0 .
When d = 0 in (5), (φ, ψ) is said to be delay-normalized [24]. Clearly there is no loss in generality in assuming this, as changing d simply reorders the mother wavelet basis functions with a shift. Delay-normalized wavelets will be assumed in the following, as this will simplify some proofs. Wavelets with symmetric filters are clearly delay-normalized.
These formulas tell us which scaling functions at scale 1 contribute in ψ andψ , a fact which will be useful. It is straightforward to find the supports of the mother wavelets from the supports of the filters (see for instance Exercise 5.16 in [21,22]). In particular, a delaynormalized wavelet can be recognized in terms of the supports by the requirement Clearly (φ 0,n , ψ), as well as (φ 0,n ,ψ), are also delay-normalized for any n, as translating φ andφ with the same n gives scaling functions for a new biorthogonal wavelet. For an orthonormal wavelet the filters and the dual filters equal, and H is orthogonal. From the deductions above one sees that supp(G 0 ) = supp(G T 1 ), in order for an orthonormal wavelet to be delay-normalized. It is straightforward to check that assumption (3) implies that (6) holds, so that this support assumption guarantees a delay-normalizedness.

Setup for Wavelets on the Interval
When restricting to an interval of the form [0, M], the wavelet bases φ m , ψ m , and the functions φ m,n , ψ m,n , will be replaced with new bases φ b m , ψ b m , and modified functions φ b m,n , ψ b m,n . Here b is short for boundary, and only those functions supported near the boundaries are modified (called left-and right edge functions). Initial candidates for the left edge scaling functions will first be defined. It will then be seen how changes of coordinates can be applied to make those functions orthonormal/biorthogonal. The right edge functions will be obtained by repeating the left edge analysis, following a mirroring operation. The following definition is a generalization of that in [9].
and the sets Some comments are in order.
-The first part of these functions are replacements of the {φ 0,k } K−1 k=K−N . Moreover, span {φ b 0,k } k≥0 is independent of the choice of {c k } N−1 k=0 , and will contain all polynomials of degree < N on (0, ∞). This follows since ∞ n=−∞ c k (n)φ 0,n is a polynomial on (−∞, ∞), and its restriction to (0, ∞) can be written as -The second part of these functions are translates of the scaling function, all supported on [0, ∞) since K ≥ −L,K ≥ −L. They are called internal functions. The index shift K − N =K −Ñ is present for technical reasons: In order to compute an m-level DWT on the interval, we will see that K must be chosen accordingly.
0 . This is reflected in the IDWT matrix in that rows are shifted K − N entries downwards. However, K − N basis functions will be added later, so that the net effect is that there is no shift. Let C be the matrix with entries C n,k = c k (n) for (n, k) C clearly has linearly independent columns, and thus full rank N . Any N rows of C give a nonsingular matrix, since any polynomial of degree N − 1 is uniquely identified from N distinct points.
Letting α be the column vector with entries α i , and using (9) one gets , it is clear that Cα = 0. But since C has linearly independent columns it follows that α = 0, so that the are linearly dependent as well. 1. now follows from the obvious fact that {φ b 0,k } N−1 k=0 and {φ 0,k } k≥K = {φ b 0,k } k≥N are linearly independent on [0, ∞). 2. and 3. follow also easily, since It is known that the modified edge functions satisfy a refinement relation (This fact will be reproved in our setting in the following), so that the new spaces V b m also give rise to a multiresolution analysis. One can thus define change of coordinate matrices H b and G b as before, replacing the counterparts on the entire real line. If the first N scaling functions/mother wavelets need modification, (4) translates to with all but the last two listed functions being modified versions. Since the unmodified functions inherit known refinability relations, the two columns listed last above are known.
For the first columns one will write with X representing the contribution of the modified functions, Z that of the internal functions, i.e.
Denoting the even-and odd-indexed columns in X and Z, by X e , Z e and X o , Z o (throughout the paper the letters e and o will indicate even and odd indices), respectively, (10) is equivalent to Since the left edge functions span the same space, regardless of the choice of polynomials, it makes sense to consider changes of coordinates between different candidates for edge functions. One can apply several such coordinate changes, in order to obtain functions φ b 0,n with desired properties. The following result concerns how X e , Z e , and C are updated by such coordinate changes (note that Lemma 1 guarantees the uniqueness of X e and Z e in a factorization of the form (11)).

Lemma 2 Assume that a change of coordinates is applied to the left edge functions. Let
0,k } the change of coordinate matrix from the second to the first basis. If (14) transforms X e , Z e , and C according to Proof One has that One also has that ⎛ (13) follows.

Finding the Left Edge Scaling Functions
First the refinement relations satisfied by the modified edge functions is established.

Lemma 3 For each choice of polynomial basis {c
with where and where C † is the generalized inverse of C. X e is nonsingular.
Here it is assumed that C † has row indices equal to the column indices of C, and vice versa.
Proof The first part of this proof corresponds to Lemma 3.1 in [11]. Since V b 0 contains all polynomials of degree less than N , C can be chosen so that Inserting 2t for t and multiplying with √ 2 one also has Comparing these and using matrix notation one sees that ⎛ where (14) for this choice of φ b 0,k , with X e and Z e having the stated indices. Since clearly X e = D, it must be nonsingular. Since (14) holds and X e is nonsingular for this initial basis, Lemma 2 says that this will be the case for any other polynomial basis as well. Now, (9) can be written ⎛ where the matrix product with G T was split into two parts on the third line. Noticing that (14) can also be rewritten as ⎛ Comparing with (17) and using the linear independence of {{φ 1,n } n≥−R+1 } on [0, ∞), one sees that Multiplying with C † to the left in the first equation gives the first equation in (15).
Some remarks on the initial choice of polynomials can be found in the technical report [2].

Change of Coordinates for Staggered Supports
We will now try to make a change of coordinates so that the new bases satisfy One says that the supports are staggered. 0,k have staggered supports, φ b 0,k will lie in this k-dimensional subspace, so that standard orthogonalization procedures give us a unique (up to signs) orthonormal basis of functions with staggered supports. Staggered supports can thus be used to single out unique boundary functions (as partially noted in [9]).
More generally we will say that This more general definition also comprises the setting in [19], and also possible supports for the mother wavelets. A coordinate change transforming the f i to functions with staggered support can clearly be interpreted in terms of bringing a matrix to echelon form.
To obtain bases with staggered supports, Lemma 2 says that one needs to find a change of coordinates P so that the lower N × N -block of CP is upper triangular. Clearly this can be achieved by means of a QR-factorization, or an LU factorization. The technical report [2] contains further details.

Change of Coordinates for Orthogonalization
In the following we will assume that N =Ñ (the deductions are a bit more complicated when N =Ñ , see Sect. 13 in the technical report [2] for further details). Since  4 One of the two sets may here contain internal scaling functions. We will see that a coordinate change can be made so that it does not affect these.
the Gramm matrix of the initial bases, one sees that, when N =Ñ (which will be assumed for simplicity in the following), Solving AV B T = F is equivalent to solving the linear system (A ⊗ B)vec(V ) = vec(F ) [16], where ⊗ is the (left) Kronecker product, and where vec(X) is the vector where the rows of X have been stacked horizontally and then transposed to a column vector. Equation (19) can therefore be written as (see also Theorem 3.2 in [18]. This paper also elaborates on the general computation of Gramm matrices), where I is the N 2 × N 2 identity matrix. 5 Denote bases by and let P = P B←C andP = PB ←C be the corresponding coordinate changes (i.e., from the old to the new bases). It is straightforward to show that Since one wants bases C andC so that ( C,C ) = I , and since upper triangular coordinate changes preserve staggered supports, one seeks upper triangular matrices P andP so that P T ( B,B )P = I , i.e., so that Now, if Y = LU is an L 1 U-factorization 6 of Y , 7 one can choose our upper triangular coordinate changes as P = (L −1 ) T andP = U −1 . In the orthogonal case where B =B and C =C, Y is positive semidefinite, and thus has a unique Cholesky factorization Y = LL T , so that one can choose P =P = (L −1 ) T as our coordinate change. 8 In the following it will always be assumed that φ b 0 andφ b 0 are biorthogonal, both with staggered supports. 5 X e has eigenvalues {2 −k−1/2 } N −1 k=0 . It follows that ρ(X e ) < 1, so that I − (X e ) T ⊗ (X e ) T is nonsingular, so that Y is unique. 6 L 1 U means that L is lower triangular with ones on the diagonal, U is upper triangular. Similarly LU 1 means that U is upper triangular with ones on the diagonal. 7 It is a major issue whether Y is nonsingular in the general case. Some special cases are handled in [11,18,20]. For Y to have a unique L 1 U-factorization one also needs the principal leading submatrices to be nonsingular. The software implementation handles these issues only numerically. Other factorizations such as a LU 1 or LDU could also be chosen. 8 It has been noted in the literature that Y can be badly conditioned. [10] proposes to use a Singular Value Decomposition of Y to address this problem. [10] does not assume staggered supports.

Stable Completion and the Left Edge Mother Wavelets
The method of stable completion finds additional mother wavelets at the boundaries, thereby modifying the bases ψ m to bases ψ b m . Only those mother wavelets supported near the boundaries are changed. Defining C b m similarly as C m in Equation (2) (i.e., replacing ψ m with ψ b m ), it is again required that both φ b m and C b m are bases for V b m (just as φ m and C m are for V m ). As initial candidates for the functions in ψ 0 andψ 0 we define ψ b 0,k = ψ 0,k+K−N and ψ b 0,k =ψ 0,k+K−Ñ . This aligns φ b 0,n and ψ b 0,n in the same way as φ 0,n and ψ 0,n . First the ψ bandψ b -functions that satisfy old refinement relations will be characterized.

Lemma 4
The following hold.
In the following it will be assumed that N 0 andÑ 0 satisfy the properties in Lemma 4, and we will set N 0 = max(N 0 ,Ñ 0 ). The set {ψ b 0,k } k≥N 0 = {ψ 0,k } k≥N 0 +K−N accounts for all but the first N 0 + K − N mother wavelets {ψ 0,k } k≥0 in V b 1 , and similarly on the dual side. Define and write A = span(A),Ã = span(Ã). The previous lemma states that A ⊂ V b 1 , andÃ ⊂ V b 1 and A andÃ, are clearly biorthogonal by construction. Denote by A ⊥ the orthogonal complement of A inṼ b 1 . If g ∈ A,g will denote the vector inÃ with the same index, and vice versa. We define, for f ∈ A andf ∈Ã,

Note that
-P equals the identity on A, and equals zero onÃ ⊥ .
-I − P equals the identity onÃ ⊥ , and equals zero on A.
-The spaces A andÃ ⊥ are linearly independent: If v ∈ V b 1 lies in both these spaces it must be on the form v = g∈A α g g, and must for allh ∈Ã satisfy g∈A α i g,h = α h = 0, where biorthogonality was used.
is the unique decomposition of v ∈ V b 1 in A ⊕Ã ⊥ . P is thus a generalization of orthogonal projection, for which A =Ã.

Lemma 5
Assume that N 0 ≥ R andÑ 0 ≥R, and that they satisfy the properties in Lemma 4. Assume also that the supports of φ b 1 andφ b 1 are staggered. Then the following hold.
,Ã is a linearly independent set.

Remark 1
The requirements N 0 ≥ R,Ñ 0 ≥R) are a bit difficult to grasp. If these are not fulfilled, the set [0, K − N + R) will not be contained in S (the way this set is defined), and this will imply that one can't find enough functions (using the strategy of the proof) to find a new basis for V b . Since one would like N 0 andÑ 0 to be as small as possible (to inherit as many of the old refinement relations as possible), the implementation computes these as Proof Only the first statement is considered, since the second statement follows from the same line of arguments. It is easily checked that |S| = N 0 + K − N , and that the last and largest entry of S is Due to staggered supports and Lemma 3, the highest index among boundary basis functions at resolutions 1 which contribute in φ b 0,0 is Also, the highest index among boundary basis functions at resolutions 1 which contribute in φ b 0,k is K − N + R + 2k. It follows that the coordinate matrix of {φ b 0,k } k≥0 and {φ b 1,k } k∈S relative to φ b 1 has different highest contributing indices. In particular, any finite set of these columns must be linearly independent, a fact which will be used in the final part of the proof.

Since clearly {ψ
where α is a vector indexed over S, γ a vector indexed over [N 0 , N 0 ). But then This means that for some vector β, and where γ was expanded to a vector with indices from [N 0 , ∞). Viewing φ b 0,n as functions on [L − R + 1, ∞) (see Equation (7)), and similarly for φ b 1,k , we can ensure that (23) holds also on (−∞, ∞), by adding a finite linear combination of functions φ 1,k , k ≤ −R (i.e., scaling functions supported on (−∞, 0) on the left hand side: Taking inner products withφ b 0,n one sees that β n = 0 for n sufficiently large. Taking inner products over (−∞, ∞) with {ψ 0,n } n≥N 0 +K−N on both sides (these may not be supported on Since also φ b 0,n ,ψ b 0,k (−∞,∞) = 0 for all k, it follows that γ n = 0 for n ≥ N 0 . That also α k = β n = 0 for k ∈ S, and for smaller n, follows since, as noted above, from the fact that any finite set of columns in the coordinate matrix of In the next section this analysis will be repeated at the right edge, and it will follow from a simple dimension count that the two mentioned linearly independent sets are in fact bases for V b 1 andṼ b 1 . If S = {k 1 , . . . , k |S] }, one therefore defines for N − K ≤ n < N 0 , We point out that when N − K < 0 this gives negative index sets. In particular the matrix G b takes the form

Lemma 6
The coordinate matrix of the {{(I − P )φ b 1,k } k∈S relative to φ b 1 is where T is the largest integer so that 2T ≤ 2N 0 − R −L − 1, i.e., Proof In the proof above it was shown that φ b 1,k ,ψ b 0,n = 0 for k ∈ S and n ≥ N 0 , so that Now, theφ b 0,n can be expressed in terms of {φ b 1,r } r≥2n+K−Ñ+L , and since the largest index in S is 2N 0 + K − N − R − 1, one can have contribution in the sum above only when n satisfies Since N − K =Ñ −K, this occurs when 2n ≤ 2N 0 − R −L − 1. This gives the expression for T . One obtains It follows that This gives the individual columns in the coordinate matrix of the {{(I − P )φ b 1,k } k∈S relative to φ b 1 , so that this matrix is I : Remark 2 AT also needs to be computed for the dual wavelet (i.e.,T =Ñ 0 + −R +L+1 2 ).

Remark 3
The previous lemma does not provide any row limits. To deduce such limits, note first that φ b 0,T can be expressed in terms of . Note also that, since by definition 2T ≥ 2N 0 − R −L − 2, it follows that SinceL + 1 ≤ R, this is ≤ 2T + R + K − N . Since the largest entry in S is 2N 0 + K − N − R − 1, and after dropping rows that are zero, this proves that the coordinate matrix of the {{(I − P )φ b 1,k } k∈S can also be written relative to {φ b 1,0 , . . . , φ b 1,2T +K−N+R }, and as with the 1-level DWT the change of coordinates from φ b i.e., ⎛ ⎜ ⎝ The Proof One gets as above (29) follows.
(29) is applied twice. First the supports of {ψ b 0,k } N−1 k=0 is made staggered by finding a coordinate change which brings X o Z o to echelon form. Secondly the mother wavelets at the left edge are bi-orthogonalized, preserving their staggeredness. One now needs the Similarly to Equation (19) it follows that Y = (X o ) TX o + (Z o ) TZ o (there is no Y on the right hand side here, however, since there are no ψ -functions on the right hand side). With Y found, one proceeds as in the end of Sect. 4 to find the required changes of coordinates.

The Right Edge
Up to now K andK have denoted the number of scaling functions at resolution 1 needed to synthesize {φ b 0,k } N−1 k=0 and {φ b 0,k }Ñ −1 k=0 . In analyzing the right edge, the flexibility in these numbers needs to be exploited, in order to obtain a Discrete Wavelet Transform on the interval for the input length in question. K andK will now be assigned different values at the left and right edge, and K L andK L will be written for the left edge values, K R andK R for the right edge values. Since the value N 0 depended on K as well, its notation is changed to N 0,L and N 0,R .
All functions will be assumed to be defined on [0, M]. The operation m(f ) = f (M − t) "mirrors" functions on [0, M] so that the left and right edges are swapped. To reuse the left edge analysis, right edge functions will be found so that their mirrors are on the form that has been considered, i.e., delay-normalized, possibly with the same supports. If this is the case, the right edge analysis simply boils down to repeating the left edge analysis with reversed filter coefficients. The following result addresses this. Some comments are in order. 1. The case L + R =L +R = 0 includes wavelets with symmetric filters, while the case L + R =L +R = 1 includes orthogonal wavelets with the support assumption (3).
In both cases the new pairs of delay-normalized functions are adjacent basis functions, the only difference being that their internal order differs (something an implementation must take into account). In summary, when L + R =L +R one can always assume equal supports at the left and right edge, and that the mirroring process preserves the ordering of the basis functions.
3. It is straightforward to prove that m(φ m,2 m M−(L+R) ) has the same support as φ m,0 , in the same way one proves the third statement. If L + R =L +R and (φ, ψ) is delay-normalized, it also follows that m(ψ m,2 m M−1 ) has the same support as ψ m,0 . In other words, when constructing a multiresolution analysis it is desirable to consider the sets {φ m,n } , since a mirroring operation on each of these simply produce new functions with the same supports. 4. When L+R =L+R = 0 the right edge analysis is still possible, but the supports at the right edge will be different from those at the left edge. Also, the relative ordering at the right edge needs to be handled differently for the wavelet and the dual wavelet, making things more complicated. The software implementation therefore handles only the case L + R = L +R, and this will be assumed in the following.
Note that, assuming that L + R = 0 or 1, from the two filters G 0 and G 1 one can compute everything (L, R,L, andR, as well as H 0 and H 1 , N andÑ ). In particular there is no need to specify the positions of the filter coefficients, since these can also be inferred. The software implementation takes advantage of this fact. In order for m(ψ 0,n ) to ensure delay-normalized mirrors, Equation (6)  Both these are linearly independent (simply repeat the arguments from Sect. 5 by also including modified functions at the right edge), and one can define

Proof of
m ) (and their duals, and one still has biorthogonality) as before. A simple count gives that It is easily checked that dim m and φ b m are alternative bases for V b m (redefining C b m in the obvious way), so that one can define the DWT/IDWT as before. The following result provides a requirement on K L + K R in order for an m-level DWT to be computable.

Remark 4
The value of M is not needed in the computations, and is eliminated in favour of dim(x). The proof of this theorem makes it clear that M must be chosen as in order for x to be the coordinates in a wavelet basis on [0, M]. In addition to Equation (31), we need also take into account that K L and K R must satisfy the requirements of Definition 1.

Proof of
which we also will have use for. Since K L − N =K L −Ñ , and L + R =L +R, the same requirement is needed for a dual m-level DWT as well.
The condition from Theorem 2 is not sufficient for an m-level DWT to be computable, however: It may be that the smallest resolution spaces do not have room for all the boundary functions needed in the construction. Details on this can be found in [2].
In summary, the software first finds (for a given m and dim(x)) K L + K R from Equation (31), and then chooses K L and K R so that they are as equal as possible.
Of particular interest are the cases where no shift in the basis functions is needed, i.e., when we can set K L = K R = N . Important cases are 1. dim(x) = M2 m + 1 when L + R =L +R = 0 (for instance wavelets with symmetric filters), 2. dim(x) = M2 m when L + R =L +R = 1 (for instance orthonormal wavelets). [19] gives a more refined construction of Biorthogonal Spline wavelets on the interval where the primal boundary functions are fixed, and directly defined from the Schoenberg-Spline basis with equidistant knots on the interval [23]. All internal primal scaling functions are included, thereby changing the multiresolution minimally. This differs from the previous part of the paper, where the leftmost internal functions may be absorbed in the constructed boundary functions, leading to more changes to the multiresolution at the boundaries (in terms of boundary functions with wilder oscillations). This is particularly the case when K andK are increased to adapt the multiresolution to a given input length. The strategy in [19], however, puts requirements on the input length (since the number of primal boundary functions is fixed). To get around this, a combination with that of the previous part of the paper is needed.

Biorthogonal Spline Wavelets on the Interval as Defined in [19]
Let us explain how one can extend the results from [19] to the general delay-normalized case, including how to make stable completions. One clearly needs s = min(R − L − 1, N) boundary adapted functions at the primal side, ands = min(R −L− 1,Ñ) boundary adapted functions at the dual side. With d := L + s −L −s ≥ 0, this modifies equations (7) and (8) in Definition 1 to so that we need to construct d functions {φ b } d−1 k=0 on the dual side (it is easily checked that d = N − 2 in the case of Spline wavelets, agreeing with the construction in [19]). These are defined recursively as This is simply a restatement of Equation (4.3) in [19] with all values of k there combined: G b :,2[0,d) represents the refinement relations for the functions {φ b } d−1 k=0 we seek to construct (these refinement relations were previously denoted byX e andZ e ). We will assume that the nonzero indices in column n ofG b e come at indices [n + 1, x + 2n] (i.e., we assume staggered supports, and that the φ b 0,n can be defined in terms of φ b 0,n+1 , . . . , φ b 0,d−1 , i.e., iteratively in terms of previously constructed functions). This is again all compatible with [19], with the additional complication of finding a number x so that unique refinement relations can be found. The next result is the delay-normalized generalization of Theorem 4.13 in [19], and explains how to find the columns inG b iteratively, and the value for x.
, and assume that z is a solution to the linear system Define column n ofG b e so that (G b e ) [n+1,x+2n],n = z, and zero elsewhere, and definẽ Then φ b 0,k ,φ b 0,n = δ k,n for all k. Moreover, the linear system is square if x = s + 1 (in which case t = n + s + 1).
The value of x here is again compatible with the statements for the Spline case in [19]. We do not prove that the coefficient matrix of the above system is invertible. As commented in [19], numerical experiments show this to be true for the most common values of N anÑ in the Spline case, but a general proof for this was not given.
0,n is one of the functions we seek to construct), only the In particular, Using the above, for 2k > 2n + x + s (i.e., k ≥ n + x+s+1 which is zero by assumption (since 2k − s > 2n + x ≥ x). Note that t = x + n if and only if which is seen to hold if x = s + 1 (this agrees with x = N in the spline case). We have that  .
Since this is upper triangular, one can proceed as in [19] to biorthogonalize while preserving staggered supports. It is not too hard to adapt the arguments from Sect. 5 to obtain stable completions also here. Simply note that the quantity K − N needs to be replaced by −L − s (compare Definition 1 and Equation (35)). This modifies the definition of the sets S andS as well. The set S is "punctured" with d values due to the addition of the {φ b 0,k } d−1 k=0 . Currently the software implementation does not support a delay-normalized generalization of the results in [19]. There are two reasons for this. First of all, keeping the primal boundary functions fixed is not compatible with adaptability to the input length, one of the main features of the software implementation (and there is no consensus on how to combine the strategy from [19] with absorption of inner functions in the boundary functions). Secondly, one would need some canonical boundary functions at the primal side in the general delay-normalized case, similarly to those one has from the Schoenberg Spline basis in the Spline case. Generally it is not clear what such canonical boundary functions should be. 9

Notes on the Implementation
This contribution spares the reader for many tedious calculations needed for the software implementation. Necessary details can be found in the technical report [2], and can be summarized as follows.
-What are the smallest input sizes so that a DWT/IDWT is possible (i.e., so that the left and right boundary functions do not interact)? -Adoption of a lifting-based approach to the interval. -Preconditioning (this was addressed in [9] for the orthogonal case, but [2] addresses this more generally). -One has freedom in how to scale the modified boundary functions. [2] computes Gramm matrices as described in [18], in order to find the norms of the modified functions, and uses these to scale them accordingly. -Computation of the Gramm matrices when N =Ñ .
[2] also explains the interface, and provides some code examples on how to use it. The code repository itself can be found at https://github.com/oyvindry/wl.

Conclusion
A unified scheme for wavelets on the interval was established, which extends known cases (such as orthonormal-and Spline-) of wavelets, with various degrees of polynomial exactness. We also considered the method of stable completion within the new scheme, and offered some new perspectives on how to establish such completions. We made a big point of using simplified notation in the construction, and building heavily on linear algebra concepts. As a consequence, the way towards a software implementation was shortened, and interested readers are encouraged to experiment further with the software implementation accompanying the paper.

Funding Note Open access funding provided by University of Oslo (incl Oslo University Hospital).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.