Aleksandar Ignjatović
email: ignjat@cse.unsw.edu.au

Chromatic Derivatives, Chromatic Expansions and Associated Function Spaces

We present the basic properties of the chromatic derivatives and the chromatic expansions as well as a motivation for introducing these notions. The chromatic derivatives are special, numerically robust linear differential operators; the chromatic expansions are the associated local expansions, which possess the best features of both the Taylor and the Nyquist expansions. This makes them potentially useful in fields involving sampled data, such as signal and image processing.

Motivation

The Nyquist-(Whittaker-Kotelnikov-Shannon) expansion f (t) = ∞ n=-∞ f (n) sin π(t -n)/π(t -n) of a π-band limited signal of finite energy f (t) ∈ BL(π) is of global nature, because it requires samples of the signal at integers of arbitrarily large absolute value. On the other hand, since signals from BL(π) are analytic functions, they can also be represented by the Taylor expansion, f (t) = ∞ n=0 f (n) (0) t n /n!. Such expansion is of local nature, because the values of the derivatives f (n) (0) are determined by the values of the signal in an arbitrarily small neighborhood of zero. While the Nyquist expansion has a central role in digital signal processing, the Taylor expansion is of very limited use there, for several reasons.

(1) Numerical evaluation of higher order derivatives of a signal from its samples is very noise sensitive; in general, one is cautioned against numerical differentiation of signals given by empirical samples. [START_REF] Ignjatovic | Chromatic derivatives and associated function spaces[END_REF] The Taylor expansion of a signal f ∈ BL(π) converges non-uniformly; its truncations are unbounded and have rapid error accumulation.

(3) The Nyquist expansion of a signal f ∈ BL(π) converges to f in BL(π) and thus the action of a filter A on any f ∈ BL(π) can be expressed using the samples of f and the impulse response A[sinc ] of A, i.e.,

A[f ](t) = ∞ n=-∞ f (n) A [sinc ] (t -n).
(

) 1
In contrast, the polynomials obtained by truncating the Taylor series do not belong to BL(π) and nothing similar to (1) holds for the Taylor expansion.

The chromatic derivatives and the chromatic expansions and approximations were introduced to obtain local signal representations which do not suffer from these problems.

Chromatic Derivatives

To explain our notions, we first consider normalized and rescaled Legendre polynomials P L n (ω) which satisfy

1 2π π -π P L n (ω) P L m (ω)dω = δ(m -n),
and then define operator polynomials

K n t = 1 i n P L n i d dt . ( 2
)
It is easy to verify that for f ∈ BL(π) and its Fourier transform f (ω) we have

K n [f ](t) = 1 2π π -π i n P L n (ω) f (ω) e i ωt dω.
Figure 1 compares the plots of P L n (ω) and ω n /π n for n = 15 to n = 18, which are the transfer functions (save a factor of i n ) of the operators K n and of the (normalized) derivatives 1/π n d n /dt n , respectively. While the transfer functions of the normalized "standard derivatives" 1/π n d n /dt n obliterate the spectrum of the signal, leaving only its edges which in practice contain mostly noise, the transfer functions of operators K n form a family of well separated, interleaved and increasingly refined comb filters. Due to their spectrum preserving property, we call the operators K n the chromatic derivatives associated with the Legendre polynomials. Both analytic estimates and empirical tests have shown that the chromatic derivatives can be accurately and robustly evaluated from samples of the signal taken at a small multiple (2 to 4) of the usual Nyquist rate, thus solving problem (1) associated with numerical evaluation of the standard derivatives, mentioned above. Chromatic expansions, on the other hand, were introduced to solve problems (2) and (3).

Chromatic Approximations

Proposition 1 Let K n be the chromatic derivatives associated with the Legendre polynomials, and let f (t) be any analytic function; then for all t,

f (t) = ∞ n=0 (-1) n K n [f ](u) K n [sinc ](t -u). (3) If f ∈ BL(π)
k ≤ n, f (k) (u) = d k /dt k CA[f, n, u](t) t=u .
Figure 2 compares the behavior of the chromatic approximation (black) of a signal f ∈ BL(π) (gray) with the behavior of the Taylor approximation of f (t) (dashed). Both approximations are of order sixteen. The plot reveals that, when approximating a signal f ∈ BL(π), a chromatic approximation has a much gentler error accumulation when moving away from the point of expansion than the Taylor approximation of the same order. Functions K n [sinc ](t) appearing in the chromatic expansion associated with the Legendre polynomials are given by K n [sinc ](t) = (-1) n √ 2n + 1 j n (πt), where j n is the spherical Bessel function of the first kind of order Thus, unlike the monomials that appear in the Taylor formula, functions K n [sinc ](t) belong to BL(π) and satisfy |K n [sinc ](t)| ≤ 1 for all t ∈ R. Consequently, the chromatic approximations are bounded on R and belong to BL(π). Also, as Proposition 1 asserts, the chromatic approximation of a signal f ∈ BL(π) converges in BL(π). Thus, if A is a filter, then A commutes with the differential operators K n and for every f ∈ BL(π), we have the following analogue of (1):

A[f ](t) = ∞ n=0 (-1) n K n [f ](0) K n [A[ sinc ]](t).
Thus, while local, the chromatic expansion possesses the features that make the Nyquist expansion useful in signal processing. This, together with numerical robustness of the chromatic derivatives, makes chromatic approximations applicable in fields involving empirically sampled data, such as digital signal and image processing. The next proposition demonstrates another remarkable feature of the chromatic derivatives which is relevant to signal processing.

Proposition 2 Let K n be the chromatic derivatives associated with the (re-scaled and normalized) Legendre polynomials, and f, g ∈ BL(π). Then

∞ n=0 K n [f ](t) 2 = ∞ -∞ f (x) 2 dx; ∞ n=0 K n [f ](t)K n [g](t) = ∞ -∞ f (x)g(x)dx; ∞ n=0 K n [f ](t)K n t [g(u -t)] = ∞ -∞ f (x)g(u -x)dx.
Thus, the sums on the left hand side of the above equations do not depend on the choice of the instant t.

Note that the above equations provide local representations of the usual norm, the scalar product and the convolution, respectively, which are defined in L 2 globally, as improper integrals.

Given the above properties of the Legendre polynomials, it is natural to ask if other families of orthonormal polynomials have similar properties. This question was answered in [START_REF] Ignjatovic | Local approximations based on orthogonal differential operators[END_REF].

General Chromatic Derivatives

Let M : P ω → R be a linear functional on the vector space P ω of real polynomials in the variable ω. Such M is called a moment functional and

µ n = M(ω n ) is the moment of M of order n. Definition 1 A moment functionals M is chromatic if it satisfies the following conditions (condition (iii) is not es- sential, but simplifies the technicalities): (i) M is positive definite; (ii) lim sup n→∞ µ 1/n n /n < ∞; (iii) M is symmetric, i.e., µ 2n+1 = 0 for all n.
For functionals M which satisfy conditions (i) and (iii) there exists a family of real polynomials {P M n (ω)} n∈N , such that P M n (ω) contains only powers of ω of the same parity as n and which are orthonormal with respect to M; i.e., for all m, n,

M(P M m (ω) P M n (ω)) = δ(m -n).
The family {P M n (ω)} n∈N is a family of orthonormal polynomials which corresponds to a symmetric positive definite moment functional M just in case there exists a sequence of positive reals {γ n } n∈N such that

P M n+1 (ω) = 1 γn ω P M n (ω) -γn-1 γn P M n-1 (ω). (4)
For every positive definite moment functional there exists a non-decreasing bounded function a(ω), called an mdistribution function, such that for the associated Stieltjes integral we have

∞ -∞ ω n da(ω) = µ n , (5) ∞ -∞ P M n (ω) P M m (ω) da(ω) = δ(m -n). ( 6
) If M is chromatic, then condition (3) implies that {P M n (ω)} n∈N is a complete system in L 2 a(ω) . Let ϕ ∈ L 2 a(ω)
; we can define a corresponding function

f ϕ : R → C by f ϕ (t) = ∞ -∞ ϕ(ω)e iωt da(ω), (7)
and one can show that (7) can be differentiated under the integral sign any number of times. Setting

K n = 1 i n P M n (ω) i d dt
we get that for all t

K n [f ϕ ](t) = ∞ -∞ i n P M n (ω) ϕ(ω) e iωt da(ω), (8) i.e., ϕ(ω)e iωt , P M n (ω) a(ω) = (-i) n K n [f ϕ ](t). Thus, ϕ(ω)e iωt = (-i) n K n [f ϕ ](t)P M n (ω)
, and by Parseval's Theorem, for every t ∈ R,

∞ n=0 |K n [f ϕ ](t)| 2 = ϕ(ω)e iωt 2 a(ω) = ϕ(ω) 2 a(ω) . Thus, if ϕ ∈ L 2 a(ω) , then the sum ∞ n=0 |K n [f ϕ ](t)| 2 con- verges to a constant function on R. If we let m(t) = ∞ -∞ e iωt da(ω), (9)
then (5) implies m (k) (0) = i k µ k . It can be shown that condition (iii) of Definition 1 implies that m(t) is analytic at every t ∈ R (moreover, it is analytic on a strip in C; see [START_REF] Ignjatovic | Chromatic derivatives and associated function spaces[END_REF]). For the chromatic approximation associated with M,

CA M [f, n, u](t) = n k=0 (-1) k K k [f ](u)K k [m](t -u), one can show that |f ϕ (t) -CA M [f ϕ , n, u](t)| < ∞ k=n+1 K k [f ϕ ](u) 2 .
Thus,

f ϕ (t) = ∞ k=0 (-1) k K k [f ϕ ](u) K k [m](t -u)
, and the convergence is uniform on R.

Definition 2 L M

2 denotes the space of functions analytic on R which satisfy

∞ k=0 K k [f ](0) 2 < ∞. Let f (t) ∈ L M 2 ; then ϕ f (ω) = ∞ k=0 (-i) k K k [f ](0)P M k (ω) belongs to L 2 a(ω)
and for all t, f (t) = ∞ -∞ ϕ f (ω) e iωt da(ω). On the space L M 2 one can now introduce locally defined norm, inner product and convolution using equations from Proposition 2, and for every fixed u, the chromatic expansion of an f ∈ L M 2 is just the Fourier series of f in the orthonormal and complete base {K n u [m(t -u)]} n∈N .

Examples

Example 1. (Legendre polynomials/Spherical Bessel functions) Let L n (ω) be the Legendre polynomials; if we set

P L n (ω) = √ 2n + 1 L n (ω/π), then π -π P L n (ω)P L m (ω) dω 2π = δ(m -n).
The corresponding recursion coefficients in equation ( 4) are given by the formula γ n = π(n+1)/ 4(n + 1) 2 -1.

In this case m(t) = sinc t, and

K n [m](t) = ( -
T n (ω) = √ 2 T n (ω/π)
for n > 0. In this case

π -π P T n (ω)P T m (ω) dω π 2 1-( ω π ) 2 = δ(n -m).
The corresponding function ( 9) is m(t) = J 0 (πt) and

K n [m](t) = (-1) n √ 2 J n (πt) for n > 0, where J n (t)
is the Bessel function of the first kind of order n. In the recurrence relation (4) the coefficients are given by γ 0 = π/ √ 2 and γ n = π/2 for n > 0. The corresponding space L M 2 consists of analytic functions whose Fourier transform f (ω) is supported in (-π, π) and satisfies

π -π 1 -(ω/π) 2 | f (ω)| 2 dω < ∞. The chromatic expansion of a function f (t) is the Neumann series f (t) = f (0)J 0 (πt) + √ 2 ∞ n=1 K n [f ](0)J n (πt).
Thus, the chromatic expansions corresponding to various families of orthogonal polynomials can be seen as generalizations of the Neumann series, while the families of corresponding functions {K n [m](t)} n∈N can be seen as generalizations (and a uniform representation) of some familiar families of special functions.

Example 3. (Hermite polynomials/Gaussian monomial functions) Let H n (ω) be the Hermite polynomials; then the polynomials given by

P H n (ω) = (2 n n!) -1/2 H n (ω) satisfy ∞ -∞ P H n (ω)P H m (ω) e -ω 2 √ π dω = δ(n -m).
The corresponding function defined by (9) is m(t) = e -t 2 /4 and K n [m](t) = (-1) n t n e -t 2 /4 / √ 2 n n!. The corresponding recursion coefficients are given by γ n = (n + 1)/2. The corresponding space L M 2 consists of analytic functions whose Fourier transform f (ω) satisfies

∞ -∞ | f (ω)| 2 e ω 2 dω < ∞.
The chromatic expansion of f (t) is just the Taylor expansion of f (t) e t 2 /4 , multiplied by e -t 2 /4 .

Weakly Bounded Moment Functionals

To study local (i.e., non-uniform) convergence of chromatic expansions, we somewhat restrict the class of moment functionals we consider.

Definition 3 Let M be a symmetric positive definite moment functional and let γ n > 0 be such that (4) holds.

(i) M is weakly bounded if there exist some M ≥ 1, some 0 ≤ p < 1 and some integer r, such that for all n ≥ 0,

1/M ≤ γ n ≤ M (n + r) p and γ n /γ n+1 ≤ M 2 .
(ii) M is bounded if there exists some M ≥ 1 such that 1/M ≤ γ n ≤ M for all n ≥ 0.

Thus, every bounded moment functional is also weakly bounded with p = 0. Functionals in our Example 1 and Example 2 are bounded. For bounded moment functionals the corresponding m-distribution a(ω) has a finite support and consequently m(t) is a band-limited signal. However, m(t) can be of infinite energy (i.e., not in L 2 ) as is the case in our Example 2. Moment functional in Example 3 is weakly bounded but not bounded (p = 1/2). We note that all important examples of classical orthogonal polynomials which correspond to weakly bounded moment functionals in fact satisfy a stronger condition 0 < lim n→∞ γ n /n p < ∞ for some 0 ≤ p < 1.

Lemma 3 If M is a weakly bounded moment functional, then lim k→∞ (µ k /k!) 1/k = 0. Thus, M is chromatic; moreover, m(z) = ∞ n=0 i n µ n z n /n! is an entire function on C.
Lemma 4 Let M be weakly bounded and p < 1 as in Definition 3(i); then for every integer k ≥ 1/(1 -p) there exists K > 0 and a polynomial P (x) such that for every n ∈ N and every z ∈ C,

|K n [m](z)| < |Kz| n P (|z|)e |Kz| k /n! 1-p .
This Lemma is used to prove the following Proposition.

Proposition 5 Let M be as in Lemma 4, f (z) an entire function and u ∈ C. If lim n→∞ |f (n) (u)/n! 1-p | 1/n = 0, then the chromatic expansion of f (z) centered at u converges everywhere to f (z), and the convergence is uniform on every disc of finite radius.

Thus, if M is bounded (p = 0) and f is an entire function, then the chromatic expansion CE[f, u](t) converges to f (t) for all t. Many well known equalities for the Bessel functions J n (t) are just the special cases of chromatic expansions. For example, the chromatic expansions of f

(t) = e iωt , f (t) = 1 and f (t) = m(t + u) yield e iωt = ∞ n=0 i n P M n (ω) K n [m](t); m(t) + ∞ n=1 n k=1 γ 2k-2 γ 2k-1 K 2n [m](t) = 1, m(t + u) = ∞ n=0 (-1) n K n [m](u)K n [m](t),
which generalize the following well known equalities:

e i ωt = J 0 (t) + 2 ∞ n=1 i n T n (ω)J n (t); J 0 (t) + 2 ∞ n=1 J 2n (t) = 1; J 0 (t + u) = J 0 (u)J 0 (t) + 2 ∞ n=1 (-1) n J n (u)J n (t).

Non-Separable Inner Product Spaces

If M is weakly bounded, the periodic functions do not belong to L M 2 ; for example, ∞ n=0 K n [sin ωt] 2 diverges. We now consider some inner product spaces in which pure harmonic oscillations have finite positive norms ( [START_REF] Ignjatovic | Chromatic derivatives and local approximations[END_REF][START_REF] Ignjatovic | Chromatic derivatives and associated function spaces[END_REF]).

Definition 4 Assume again that M is weakly bounded and let p be as in Definition 3. We denote by C M the vector space of analytic functions such that the sequence

ν f n (t) = 1/(n + 1) 1-p n k=0 K k [f ](t) 2 converges uniformly on every finite interval. Proposition 6 Let f, g ∈ C M and σ f g n (t) = 1/(n + 1) 1-p n k=0 K k [f ](t)K k [g](t)
; then the sequence {σ f g n (t)} n∈N converges to a constant function. In particular, ν f n (t) is constant. Corollary 7 Let C M 0 be the vector space consisting of analytic functions f (t) such that lim n→∞ ν f n (t) = 0; then in the quotient space C M 2 = C M /C M 0 the limit lim n→∞ σ f g n (t) is independent of t and defines a scalar product on C M 2 . Proposition 8 Let M correspond to Chebyshev polynomials as in our Example 2; then functions f ω (t) = √ 2 sin ωt and g ω (t) = √ 2 cos ωt for all 0 < ω < π form an uncountable orthonormal system of vectors in C M 2 . Proposition 9 Let M correspond to Hermite polynomials as in our Example 3; then for all ω > 0 functions f ω (t) = sin ωt and g ω (t) = cos ωt form an uncountable orthogonal system of vectors in C M 2 , and f ω M = g ω M = e ω 2 /2 / 4 √ 2π.

Conjecture 1 Assume that for some 0 ≤ p < 1 the recursion coefficients γ n in (4) are such that γ n /n p converges to a finite positive limit. Then, for the corresponding family of orthogonal polynomials we have 0 < lim n→∞ 1/(n + 1) 1-p n k=0 P M k (ω) 2 < ∞ for all ω in the support sp(a) of the corresponding mdistribution function a(ω). Thus, in the corresponding space C M 2 all pure harmonic oscillations with positive frequencies ω ∈ sp(a) have finite positive norm and are mutually orthogonal.

Detailed presentation of the theory of chromatic derivatives can be found in our references; preprints of some unpublished manuscripts are available at http://www.cse.unsw.edu.au/˜ignjat/diff.

Figure 1 :

 1 Figure 1: Graphs of P L n (ω) (left) and ω n /π n (right), for n = 15 -18.

Figure 2 :

 2 Figure 2: Chromatic approximation (black) and Taylor's approximation (dashed) of a signal from BL(π) (gray).

  1) n (2n + 1) j n (πt), where j n (x) is the spherical Bessel function of the first kind of order n. The corresponding space L M 2 consists of all analytic functions which belong to L 2 and have a Fourier Transform supported in [-π, π]. Example 2. (Chebyshev polynomials of the first kind/Bessel functions) Let P T n (ω) be the family of orthonormal polynomials obtained by normalizing and rescaling the Chebyshev polynomials of the first kind, T n (ω), by setting P T 0 (ω) = 1 and P

  the series converges uniformly and in L 2 .

	The series in (3) is denoted by CE[f, u](t) and is called
	the chromatic expansion of f (t) associated with the Leg-
	endre polynomials; a truncation of this series up to first
	n + 1 terms is denoted by CA[f, n, u](t) and is called a
	chromatic approximation of f (t). Just like a Taylor ap-
	proximation, a chromatic approximation is also a local ap-
	proximation: its coefficients are the values of differential
	operators K m [f ](u) at a single instant u, and for all


Vesa Lehtinen
email: vesa.lehtinen@tut.fi

Markku Renfors
email: markku.renfors@tut.fi

Structures for Interpolation, Decimation, and Nonuniform Sampling Based on Newton's Interpolation Formula

The variable fractional-delay (FD) filter structure by Tassart and Depalle performs Lagrange interpolation in an efficient way. We point out that this structure directly corresponds to Newton's interpolation (backward difference) formula, hence we prefer to refer to it as the Newton FD filter. This structure does not function correctly when the fractional delay is made time-variant, e.g., in sample rate conversion. We present a simple modification that enables time-variant usage such as fractional sample rate conversion and nonuniform resampling. We refer to the new structure as the Newton (interpolator) structure. Almost all advantages of the Newton FD structure are preserved. Furthermore, we suggest that by transposing the Newton interpolator we obtain the transposed Newton structure which can be used in decimation as well as reconstruction of nonuniformly sampled signals, analogously to the transposed Farrow structure. The presented structures are a competitive alternative for the Farrow structure family when low complexity and flexibility are required.

Introduction

In [START_REF] Tassart | Fractional delays using Lagrange interpolators[END_REF][2] [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF], Tassart and Depalle as well as Candan derive an efficient implementation structure for FD filters, depicted in Fig. 1, from Lagrange's interpolation formula. It turns out that the obtained filter structure directly corresponds to Newton's (backward difference) interpolation formula [START_REF] Weisstein | Newton's Backward Difference Formula[END_REF] (with some subexpression sharing) which indeed is equivalent with Lagrange interpolation [START_REF] Meijering | A chronology of interpolation: From ancient astronomy to modern signal and image processing[END_REF]. Newton's backward difference formula is ,

where [START_REF] Tassart | Analytical approximations of fractional delays: Lagrange interpolators and allpass filters[END_REF] is the rising factorial, and is the backward difference operator such that and , resulting in [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF] Newton's backward difference formula provides an efficient means to realise piecewise-polynomial interpolation for DSP. Its complexity is only O(M) (where is the interpolator order)-cf. equivalent Lagrange implementations based on the Farrow structure [START_REF] Farrow | A continuously variable digital delay element[END_REF] having O(M 2 ) complexity [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF]. The subfilters are multiplier-free and extremely simple. The structure is modular, as highlighted by the grey shading in Fig. 1, and the interpolator order can be changed in real time [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF].

Unfortunately, the structure presented in Fig. 1 does not function correctly in sample rate conversion (SRC). Because the multiplications are performed between the subfilters, making them time-variant will result in incorrect output. This is because each output sample should only depend on the current value of the delay parameter D; in Fig. 1, past values of D contribute to the output through the delayed paths through the subfilters. Therefore, the structure in Fig. 1 is only useful in single-rate, time-invariant or slowly-varying fractional-delay filtering.

We propose a slightly modified structure that allows arbitrary resampling, including increasing the sample rate by arbitrary, also fractional, factors (fractional interpolation). We also point out that the structure can be transposed to obtain a decimator structure that possesses all the advantages of the Newton interpolation structure.

This work was supported by the Graduate School in Electronics, Telecommunications and Automation (GETA).

f t τ + ( ) τ m ( ) ∆ m f t ( ) m! ---------------------------- m 0 = ∞ ∑ = τ m ( ) τ k + ( ) k 0 = m 1 - ∏ = ∆ ∆ m f t ( ) ∆ m 1 -f t ( ) ∆ m 1 -f t 1 - ( ) - = ∆ 0 f t ( ) f t ( ) = ∆ m f t ( ) 1 - ( ) k m k     f t k - ( ). k 0 = m ∑ = M

The Newton structure for interpolation

In order to allow fractional SRC and arbitrary resampling, the Newton structure must work correctly with a time-variant fractional delay. This is achieved through two simple steps: (i) We invert the summation order at the output part of the structure from that presented in [START_REF] Tassart | Fractional delays using Lagrange interpolators[END_REF][3] (this was already done in [START_REF] Tassart | Analytical approximations of fractional delays: Lagrange interpolators and allpass filters[END_REF]). (ii) The time-varying multiplications can now be implemented in the high-rate part between the adders. The improved structure is shown in Fig. 2. We refer to it as the Newton interpolator structure or the Newton structure for short. Also the improved structure is modular, permitting changing the interpolator order in real time. In single-rate FD filtering, the improved structure is equivalent to [START_REF] Tassart | Fractional delays using Lagrange interpolators[END_REF][2] [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF]. In Fig. 2, the H&S blocks stand for hold & sample, i.e., each output sample obtains the value of the previously arrived input sample.

In fractional interpolation, i.e., increasing the sample rate by a fractional factor, we use the common notation illustrated in Fig. 3. The time interval between the previous input sample and the next output sample to be generated is expressed using the fractional interval variable which is normalised with respect to the input sample interval so that . Interpolation of uniformly spaced input samples can be modelled as convolution [START_REF] Meijering | A chronology of interpolation: From ancient astronomy to modern signal and image processing[END_REF], leading to the generic model depicted in Fig. 4 [START_REF] Crochiere | Multirate Digital Signal Processing[END_REF]. The continuous-time (CT) linear time-invariant (LTI) model filter is piecewise polynomial, with pieces, each with duration equal to the input sample interval . Hence the impulse response length is .

The composite transfer function of cascaded subfilters is [START_REF] Weisstein | Newton's Backward Difference Formula[END_REF] cf. [START_REF] Candan | An efficient filtering structure for Lagrange interpolation[END_REF]. The output of the interpolator is µ H CT f ( )

1-z

1 - D+M -1 - M ------------------------ . . . . . . 1-z 1 - -D+1 2 -------------- 1 z 1 - - D - 1-z 1 - -D+2 3 --------------
1-z 1 - H&S D t ( )-M +1 1 M ⁄ . . . . . . 1-z 1 - H&S D t ( )-2 1 3 ⁄ 1-z 1 - H&S D t ( ) 1 1-z 1 - H&S D t ( )-1 1 2 ⁄ - - - - µ µ [0 1) , ∈ M 1 + T in M 1 + ( )T in
lT out l+1 ( )T out k T in k +1 ( )T in µ l+1 T in k -1 ( )T in l-1 ( )T out µ l-1 T in Input samples Output samples T in
DT CT x n [ ] x CT t ( )
y CT t ( ) y n [ ] y CT nT out ( ) = x CT t ( ) 1 F in -------- x n [ ]δ t n F in -------- -     n ∑ = @ F in @ F out m 1 z 1 - - ( ) m 1 - ( ) n m n     z n -, n 0 = m ∑ = (2.1)
where [START_REF] Meijering | A chronology of interpolation: From ancient astronomy to modern signal and image processing[END_REF] for , and

is the falling factorial. The delay of the interpolator is . The parameter can be chosen quite freely, but the best amplitude response and linear phase response are obtained with [START_REF] Tassart | Fractional delays using Lagrange interpolators[END_REF]. The continuous-time model impulse response of the interpolator is then (cf. the expression of the filter input in Fig. 4) [START_REF] Crochiere | Multirate Digital Signal Processing[END_REF] The reversed summation order in the high-rate part comes with a price: the structure is more costly to pipeline than those in [START_REF] Tassart | Fractional delays using Lagrange interpolators[END_REF][3] because the signal paths cannot share pipeline registers.

The transposed Newton structure

There exists a duality 1 between decimation and interpolation that allows transforming a decimator into an interpolator and vice versa through network transposition [START_REF] Crochiere | Multirate Digital Signal Processing[END_REF]. By transposing the Newton interpolator, we obtain the structure depicted in Fig. 5. We refer to this as the transposed Newton structure. The transpose is obtained by inverting the flow direction of all signals and replacing each block with its dual. For instance, the H&S block is replaced with the accumulate & dump (A&D) block, which sums up all its input samples since the previous output sample. This is also the most straightforward way to obtain the transposed Farrow structure from the Farrow structure 2 [START_REF] Babic | Implementation of the transposed Farrow structure[END_REF].

The output samples of the transposed Newton structure are uniformly spaced, but the input samples may arrive at arbitrary time instants. The generic SRC model (Fig. 4) is valid also for the transposed Newton structure. The model impulse response is again piecewise-polynomial, now with the piece duration equal to the output sample interval. The model impulse response is obtained by replacing with in [START_REF] Crochiere | Multirate Digital Signal Processing[END_REF] and redefining according to Fig. 6 (reflecting the duality between decimation and interpolation). For an input sample arriving at time instant , the fractional interval is [START_REF] Hentschel | Continuous-time digital filters for sample-rate conversion in reconfigurable radio terminals[END_REF] For fractional decimation, the fractional interval for the l th input sample is [START_REF] Babic | Implementation of the transposed Farrow structure[END_REF] The impulse response in the generic model is now [START_REF] Marvasti | Recovery of signals from nonuniform samples using iterative methods[END_REF] 1. There exist a number of definitions for duality, including the adjoint. Here we use the generalised duality/transpose as defined in [START_REF] Crochiere | Multirate Digital Signal Processing[END_REF]. ----------------------m -----------------------. A&D

y k µ + ( )T in ( ) h n µ + ( )T in ( ) x k n - [ ] n 0 = M ∑ = x k n - [ ] 1 - ( ) n 1 - ( ) m m n     D 0 µ - ( ) m m! -
n = M ∑ n 0 = M ∑ = m n     0 n 0 < n m > ∨ , = m 0 ≥ x ( ) m x k - ( ) k 0 = m 1 - ∏ = D 0 T in D 0 D 0 M 1 + ( ) 2 ⁄ = h n µ + ( )T in ( ) 1 T in ------- 1 - ( ) n m + m n     D 0 µ - ( ) m m! -
1-z 1 - A&D D t ( )-M +1 1 M ⁄ . . . . . . 1-z 1 - A&D D t ( )-1 1 2 ⁄ 1-z 1 - A&D D t ( )-2 1 3 ⁄ 1-z 1 - A&D D t ( ) 1 - - - - 2.
The structure in [START_REF] Hentschel | Continuous-time digital filters for sample-rate conversion in reconfigurable radio terminals[END_REF] (transposed structure I in [START_REF] Babic | Implementation of the transposed Farrow structure[END_REF]) is not the true transpose of the Farrow structure even though the duality of responses holds. µ

T in T out µ

lT in l-1 ( )T in k T out k -1 ( )T out µ l-1 T out k +1 ( )T out l+1 ( )T in µ l+1 T out
Input samples Output samples

T out t µ t ( ) t T out --------- t T out --------- - [0 1). , ∈ = µ l lT in T out --------- lT in T out ---------. - = h n µ + ( )T out ( ) 1 T out --------- 1 - ( ) n m + m n     D 0 µ - ( ) m m! ------------------------ m n = M ∑ = with integer Again,
for the best response. In the frequency response, the model filter has zeros at each (nonzero) integer multiple of the output sample rate, hence realising antialiasing regardless of the decimation factor.

The transposed Newton structure is able to receive input samples at arbitrary time instants, which makes it a potential building block for reconstruction of signals from nonuniformly spaced samples (e.g., in algorithms like [START_REF] Marvasti | Recovery of signals from nonuniform samples using iterative methods[END_REF] [START_REF] Marvasti | Reconstruction of speech signals with lost samples[END_REF]), as earlier suggested for the transposed Farrow structure in [START_REF] Babic | Reconstruction of non-uniformly sampled signal using transposed Farrow structure[END_REF].

The transposed Newton structure shares the advantages and disadvantages of the Newton interpolator, such as modularity, complexity and the inefficient zero locations.

Computational complexity

In interpolation by factor the Newton structure will perform additions and multiplications per input sample on average. In decimation by the transposed Newton structure will perform additions and multiplications per output sample. The first term in the addition count comes from the A&D block. Multiplication by a constant inverse of a small integer requires only few additions/subtractions.

Unambiguous complexity comparison between the proposed structures and alternatives, mainly the Farrow family, would require specifying the implementation technology and the SRC factor. However, the following points can be made: (i) The basis multipliers are more complex in the Newton structures (integer part present in the time-variant coefficients) than in Farrow structures (no integer part). Hence, large SRC factors are unfavourable to the Newton family. (ii) If the Lagrange response suffices, the ultimate simplicity of the subfilters makes the Newton family superior to the Farrow structure when the SRC factor is small. (iii) The response of the Newton structures can be improved only by increasing the order (i.e., number of stages). In designs with a low oversampling factor and/or strict performance requirements, this may lead to a very high filter order. In such cases, an optimised Farrow design with a non-Lagrange response will have a lower complexity and smaller delay.

Conclusions

The proposed structures allow efficient piecewise Newton interpolation for SRC and arbitrary resampling as well as its dual for decimation and reconstruction of nonuniformly sampled signals. The advantages of the proposed structures include low, O(M) complexity (high orders are feasible at the cost of a long delay), very simple subfilters and run-time adjustability of the filter order. As a drawback, the basis multipliers running at the high-rate end of the filter have longer wordlengths than in the Farrow counterparts.

Due to their simplicity, the Newton structures may be useful as building blocks of more complicated algorithms for interpolation, decimation, and reconstruction of nonuniformly sampled signals.

Figure 1 .

 1 Figure 1. The fractional-delay filter structure proposed in [1][3], based on Newton's interpolation formula.

Figure 2 .

 2 Figure 2. The Newton interpolator structure suitable for sample rate conversion. The hold & sample (H&S) blocks perform the sampling at the output sample instants.

H&S

Figure 3 .

 3 Figure 3. Definition of the fractional interval for interpolation.

Figure 4 .

 4 Figure 4. The generic model for SRC by arbitrary factors.

Figure 5 .

 5 Figure 5. The transposed Newton structure for decimation and reconstruction of signals from nonuniformly spaced samples.

Figure 6 .

 6 Figure 6. Definition of the fractional interval for the transposed structure (dual of interpolation).


Wang-Q Lim
email: wlim@mathematik.uni-osnabrueck.de

Discrete Shearlet Transform : New Multiscale Directional Image Representation

It is now widely acknowledged that analyzing the intrinsic geometrical features of an underlying image is essentially needed in image processing. In order to achieve this, several directional image representation schemes have been proposed. In this report, we develop the discrete shearlet transform (DST) which provides efficient multiscale directional representation. We also show that the implementation of the transform is built in the discrete framework based on a multiresolution analysis. We further assess the performance of the DST in image denoising and approximation applications. In image approximation, our adaptive approximation scheme using the DST significantly outperforms the wavelet transform (up to 3.0dB) and other competing transforms. Also, in image denoising, the DST compares favorably with other existing methods in the literature.

Introduction

Sharp image transitions or singularities such as edges are expensive to represent and intergrating the geometric regularity in the image representation is a key challenge to improve state of the art applications to image compression and denoising. To exploit the anisotropic regularity of a surface along edges, the basis must include elongated functions that are nearly parallel to the edges. Several image representations have been proposed to capture geometric image regularity. They include curvelets [START_REF] Candes | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF], contourlets [START_REF] Do | The contourlet transform: An efficient directional multiresolution image representation[END_REF] and bandelets [START_REF] Peyre | Discrete Bandelets with Geometric Orthogonal Filters[END_REF]. In particular, the construction of curvelets is not built directly in the discrete domain and they do not provide a multiresolution representation of the geometry. In consequence, the implementation and the mathematical analysis are more involved and less efficient. Contourlets are bases constructed with elongated basis functions using a combination of a multiscale and a directional filter bank. However, contourlets have less clear directional features than curvelets, which leads to artifacts in denoising and compression. Bandelets are bases adapted to the function that is represented. Asymptotically, the resulting bandelets are regular functions with compact support, which is not the case for contourlets. However, in order to find bases adapted to an image, the bandelet transform searches for the optimal geometry. For an image of N pixels, the complexity of this best bandelet basis algorithm is O(N 3/2 ) which requires extensive computation [START_REF] Peyre | Discrete Bandelets with Geometric Orthogonal Filters[END_REF]. Recently, a new representation scheme has been introduced [START_REF] Labate | Sparse Multidimensional Representation using Shearlets[END_REF]. These so called shearlets are frame elements which yield (nearly) optimally sparse representations [START_REF] Guo | Optimally Sparse Multidimensional Representation using Shearlets[END_REF]. This new representation system is based on a simple and rigorous mathematical framework which not only provides a more flexible theoretical tool for the geometric representation of multidimensional data, but is also more natural for implementations. As a result, the shearlet approach can be associated to a multiresolution analysis [START_REF] Labate | Sparse Multidimensional Representation using Shearlets[END_REF]. However constructions proposed in [START_REF] Labate | Sparse Multidimensional Representation using Shearlets[END_REF] do not provide compactly supported shearlets and this property is essentially needed especially in image processing applications. In fact, in order to capture local singularities in images efficiently, basis functions need to be well localized in the spatial domain. In this report, we construct compactly supported shearlets and show that there is a multiresolution analysis associated with this construction. Based on this, we develop the fast discrete shearlet transform (DST) which provides efficient directional representations.

Shearlets

A family of vectors {ϕ n } n∈Γ constitutes a frame for a Hilbert space H if there exist two positive constants A, B such that for each f ∈ H we have

A f 2 ≤ n∈Γ | f, ϕ n | 2 ≤ B f 2 .
In the event that A = B, the frame is said to be tight. Let us next introduce some notations that we will use throughout this paper. For f ∈ L 2 (R d ), the Fourier transform of f is defined by

f (ω) = R d f (x)e -2πix•ω dx.
Also, for t ∈ R d and A ∈ GL d (R), we define the following unitary operators:

T t (f )(x) = f (x -t) and D A (f )(x) = |A| -1 2 f (A -1 x). Finally, for q ∈ ( 1 2
, 1] and a > 1, we define

A 0 = a q 0 0 a 1 2
and

B 0 = 1 1 0 1 (1)
and

A 1 = a 1 2 0 0 a q and B 1 = 1 0 1 1 . ( 2
)
We are now ready to define a shearlet frame as follows.

For c ∈ R + , ψ 1 0 , . . . , ψ L 0 , ψ 1 1 , . . . , ψ L 1 ∈ L 2 (R 2 ) and φ ∈ L 2 (R 2 ), we define

Ψ 0 c = {ψ i,0 jkm : j, k ∈ Z, m ∈ Z 2 , i = 1, . . . , L}, Ψ 1 c = {ψ i,1 jkm : j, k ∈ Z, m ∈ Z 2 , i = 1, . . .

, L}, and

Ψ 2 c = {T cm φ : m ∈ Z 2 } ∪{ψ i,0 jkm : j ≥ 0, -2 j ≤ k ≤ 2 j , m ∈ Z 2 , i = 1, . . . , L} ∪{ψ i,1 jkm : j ≥ 0, -2 j ≤ k ≤ 2 j , m ∈ Z 2 , i = 1, . . . , L} where ψ i, jkm = D A -j B -k T cm ψ i (3) for = 0, 1, m ∈ Z 2 , i = 1, . . . , L and j, k ∈ Z. If Ψ p c is a frame for L 2 (R 2
), then we call the functions ψ i, jkm in the system Ψ p c shearlets. Observe that each element ψ i, jkm in Ψ p c is obtained by applying an anisotropic scaling matrix A and a shear matrix B to fixed generating functions ψ i . This implies that the system Ψ p c can provide window functions which can be elongated along arbitrary directions. Therefore, the geometrical structures of singularities in images can be efficiently represented and analyzed using those window functions. In fact, it was shown that 2-dimensional piecewise smooth functions with C 2 -singularities can be approximated with nearly optimal approximation rate using shearlets. We refer to [START_REF] Guo | Optimally Sparse Multidimensional Representation using Shearlets[END_REF] for details. Furthermore, one can show that shearlets can completely analyze the singular structures of piecewise smooth images [START_REF] Guo | Edge Analysis and identification using the Continuous Shearlet Transform[END_REF]. In fact, this property of shearlets is useful especially in signal and image processing, since singularities and irregular structures carry essential information in a signal. For example, discontinuities in the intensity of an image indicate the presence of edges. Figure 1 displays examples of shearlets which can be elongated along arbitrary direction in the spatial domain.

Construction of Shearlets

In this section, we will introduce some useful sufficient conditions to construct compactly supported shearlets. Using these conditions, we will show that the system Ψ p c can be generated by simple separable functions associated with a multiresolution analysis. Furthermore, this leads to the fast DST, and we will discuss this in the next section. We first discuss sufficient conditions for the existence of compactly supported shearlets. For this, let α > max (1, (1 -p)γ) and γ > max α+1 p , 1 1-p be fixed positive numbers for 0 < p < 1. We choose α , γ > 0 such that α ≥ α + γ and γ ≥ α -α + γ. Then we obtain the following results [START_REF] Lim | Compactly Supported Shearlet Frames and Their Applications[END_REF].

ψ i 0 (x 1 , x 2 ) = γ i (x 1 )θ(x 2 ) such that |γ i (ω 1 )| ≤ K 1 |ω 1 | α (1 + |ω 1 | 2 ) γ /2 and | θ(ω 1 )| ≤ K 2 (1 + |ω 1 | 2 ) -γ /2 . If ess inf |ω 1 |≤1/2 | θ(ω 1 )| 2 ≥ K 3 > 0 (4)
and

ess inf a -q ≤|ω 1 |≤1 L i=1 |γ i (ω 1 )| 2 ≥ K 4 > 0, (5)
then there exists c 0 > 0 such that Ψ 0 c is a frame for L 2 (R 2 ) for all c ≤ c 0 .

Observe that the functions ψ 1 0 , . . . , ψ L 0 are separable functions, and the one-dimensional scaling function θ and wavelets γ i can be chosen with sufficient vanishing moments in this case. We now show some concrete examples of compactly supported shearlets using Theorem 3.1. Assume that a = 4 and q = 1 in (1) and ( 2). Let us consider a box spline [START_REF] Candes | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] of order m defined as follows.

θm (ω 1 ) = sin πω 1 πω 1 m+1 e -i ω1 ,
where = 1 if m is even, and = 0 if m is odd. Obviously, we have the following two scaling equation:

θm (2ω 1 ) = m 0 (ω 1 ) θm (ω 1 )
and m 0 (ω 1 ) = (cos πω 1 ) m+1 e -i πω 1 .

Let α and γ be positive real numbers as in Theorem 3.1. We now define

ψ1 0 (ω) = (i) √ 2 sin πω 1 θm (ω 1 ) θm (ω 2 )
and

ψ2 0 (ω) = (i) sin πω 1 2 θm ( ω 1 2 ) θm (ω 2 ),
where ≥ α and m + 1 ≥ γ . Then, by Theorem 3.1, ψ 1 0 and ψ 2 0 generate a frame Ψ 0 c for c ≤ c 0 with some c 0 > 0. There are infinitely many possible choices for and m. For example, one can choose = 9 and m = 11.

Define φ(x 1 , x 2 ) = θ m (x 1 )θ m (x 2 ), ψ1 1 (ω) = (i) √ 2 sin πω 2 θm (ω 2 ) θm (ω 1 )
and

ψ2 1 (ω) = (i) sin πω 2 2 θm ( ω 2 2
) θm (ω 1 ).

Then similar arguments show that ψ 1 1 and ψ 2 1 generate a frame Ψ 1 c for c ≤ c 0 with some c 0 > 0. Furthermore, the functions φ, ψ i for = 0, 1 and i = 1, 2 generate a frame Ψ 2 c with c ≤ c 0 for some c 0 > 0.

Discrete Shearlet Transform

In the previous section, we constructed compactly supported shearlets generated by separable functions associated with a multiresolution analysis. In this section, we will show that this multiresolution analysis leads to the fast DST which computes f, ψ i, jkm . To be more specific, we let a = 4 and q = 1 in (1) and ( 2). For notational convenience, we let n

= (n 1 , n 2 ), m = (m 1 , m 2 ), d = (d 1 , d 2 ) ∈ Z 2 and
I 2 be a 2 by 2 identity matrix. Let θ ∈ L 2 (R) be a compactly supported function such that {θ(• -n 1 ) : n 1 ∈ Z} is an orthonormal sequence and

θ(x 1 ) = n 1 ∈Z h(n 1 ) √ 2θ(2x 1 -n 1 ). (6)
Define

γ(x 1 ) = n 1 ∈Z g(n 1 ) √ 2θ(2x 1 -n 1 ) (7)
such that γ has sufficient vanishing moments and the pair of the filters h and g is a pair of conjugate mirror filters. We assume that γ and θ satisfy decay conditions ( 4) and ( 5) in Theorem 3.1. We also define

φ(x 1 , x 2 ) = θ(x 1 )θ(x 2 ), ψ 1 (x 1 , x 2 ) = γ(x +1 )θ(x 2-) (8)
and

ψ 2 (x 1 , x 2 ) = 2 -1 2 γ( x +1 2 )θ(x 2-) (9)
for = 0, 1. Then Theorem 3.1 can be easily generalized to show that the functions ψ 1 0 , ψ 2 0 , ψ 1 1 , ψ 2 1 and φ generate a shearlet frame Ψ 2 c with c < c 0 for some c 0 > 0. Let J be a positive odd integer. Based on a multiresolution analysis associated with the two-scale equation ( 6), we can now easily derive a fast algorithm for computing shearlet coefficients f, ψ i, jkm for = 0, 1,j = 1, . . . , J-1 2 , and -2 j ≤ k ≤ 2 j as follows. First, assume that where

f = n∈Z 2 f J (n)D 2 -J I2 T n φ (10)
f J (n) = f, D 2 -J I 2 T n φ . For h = 0, 1, let us define maps D k,j h : 2 (Z 2 ) → 2 (Z 2 ) by (D k,j h x)(d) = m∈Z 2 d k,j h (d, m)x(m)
where

d k,j h (d, m) = D B k/2 j h T m φ, T d φ and x ∈ (Z 2 ).
Also we define

H(ω 1 ) = n1 h(n 1 )e -2iπω 1 and G(ω 1 ) = n 1 g(n 1 )e -2iπω1 .
Finally, we let h j , g 0 j and g 1 j be the Fourier coefficients of

       H j (ω 2 ) = J-j-1 k=0 H 2 k ω 2 for J -j > 0, G 0 j (ω 1 ) = J-2j-2 k=0 H(2 k ω 1 )G(2 J-2j-1 ω 1 ), G 1 j (ω 1 ) = J-2j-1 k=0 H(2 k ω 1 )G(2 J-2j ω 1 ), (11) respectively
. Then we obtain          f, ψ 1,0 jkm = (((D k,j 0 f J ) * r h j ) ↓2 J-j * c g 0 j ) ↓2 J-2j (m), f, ψ 2,0 jkm = (((D k,j 0 f J ) * r h j ) ↓2 J-j * c g 1 j ) ↓2 J-2j+1 (m), f, ψ 1,1 jkm = (((D k,j 1 f J ) * c h j ) ↓2 J-j * r g 0 j ) ↓2 J-2j (m), f, ψ 2,1 jkm = (((D k,j 1 f J ) * c h j ) ↓2 J-j * r g 1 j ) ↓2 J-2j+1 (m), (12)
where * c and * r are convolutions along the vertical and horizontal axes respectively, ↓ 2 j is the downsampling by 2 j and h(n) = h(-n) for given filter coefficients h(n). From (12), we observe that the shearlet transform f, ψ i, jkm is the application of the shear transformation D B k/2 j to f ∈ L 2 (R 2 ) followed by the wavelet transform associated with anisotropic scaling matrix A . In this case, applying D k,j to f J ∈ 2 (Z 2 ) corresponds to applying the shear transform D B k/2 j in the discrete domain. Thus we simply replace the operator D k,j by the discrete shear transform P k,j for f J ∈ 2 (Z 2 ), where we define the discrete shear transforms P 0 k,j and P 1 k,j as follows:

(P 0 k,j f J )(n) = f J n 1 + (k/2 j )n 2 , n 2 , (P 1 k,j f J )(n) = f J n 1 , n 2 + (k/2 j )n 1 . ( 13
)
Let M be a fixed positive integer. Since P 0 k,j and P 1 k,j are unitary operators on (Z 2 ), we can extend the shearlet transform defined in (12) to a linear transform S consisting of finitely many orthogonal transforms S M k and SM k where

S M k (f J ) = WP 0 k,M (f J ) and SM k (f J ) = WP 1 k,M (f J )
and W and W are the wavelet transform associated with an anisotropic sampling matrices A 0 and A 1 , respectively. For the precise definitions of W and W, we refer to [START_REF] Lim | Compactly Supported Shearlet Frames and Their Applications[END_REF].

In this case, the linear transform S, which we call DST, is defined by

S = (S M -2 M , . . . , S M 2 M , SM -2 M , . . . , SM 2

Image Approximation Using DST

In this section, we present some results of the DST in image compression applications. In this case, we use adaptive image representation using the DST. The main idea of this is similar to the matching pursuit introduced by Mallat and Zhong [START_REF] Mallat | Matching Pursuits With Time-Frequency Dictionaries[END_REF]. The matching pursuit selects vectors one by one from a given basis dictionary at each iteration step. On the other hand, our approximation scheme searches the optimal directional index k 0 at each iteration step so that corresponding the orthogonal transform S M k0 or SM k0 provides an optimal nonlinear approximation with P nonzero terms among all possible 2 M +2 + 2 orthogonal transforms in S. For a detailed description of this algorithm, we refer to [START_REF] Lim | Compactly Supported Shearlet Frames and Their Applications[END_REF]. For numerical tests, we compare the performance of the DST to other transforms such as the discrete biorthogonal CDF 9/7 wavelet transform (DWT) [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] and contourlet transform (CT) [START_REF] Do | The contourlet transform: An efficient directional multiresolution image representation[END_REF] in image compression (see Figure 3). We used only 2 directions (horizontal and vertical) and 4 level decomposition for our DST. In this case, our numerical tests indicate that only a few iterations (1-5) can give significant improvement over other transforms and computing time is comparable to the wavelet transform. For more results, we refer to [START_REF] Mallat | Matching Pursuits With Time-Frequency Dictionaries[END_REF].

Conclusion

We have constructed compactly supported shearlet systems which can provide efficient directional image representations. We also have developed the fast discrete implementation of shearlets called the DST. This algorithm consists of applying the shear transforms in the discrete domain followed by the anisotropic wavelet transforms. Applications of our proposed transform in image approximation and denoising were studied. In image approximation, the results obtained with our adaptive image representation using the DST are significantly superior to those of other transforms such as the DWT and CT both visually and with respect to PSNR. In denoising, we studied the performance of the DST coupled with a (partially) translation invariant hard tresholding estimator. Our results indicate that the DST consistently outperforms other competing transforms. For detailed numerical results, we refer to [START_REF] Lim | Compactly Supported Shearlet Frames and Their Applications[END_REF].

Figure 1 :

 1 Figure 1: Examples of shearlets in the spatial domain. The top row illustrates shearlet functions ψ i,0 jk0 associated with matrices A 0 and B 0 in (1). The bottom row shows shearlet functions ψ i,1 jk0 associated with matrices A 1 and B 1 in (2)..

Figure 2 :

 2 Figure 2: Examples of anisotropic discrete wavelet decomposition: (a) Anisotropic discrete wavelet decomposition by W, (b) Anisotropic discrete wavelet decomposition by W.

  M ) for a given M ∈ Z + . Notice that redundancy of the DST is K = 2 M +2 + 2 and the DST merely requires O(KN ) operations for an image of N pixels. It is obvious that the inverse DST is simply the adjoint of S with normalization.

Figure 3 :

 3 Figure 3: Compression results of 'Barbara' image of size 512 × 512: The image is reconstructed from 5024 most significant coefficients. Top left: Zoomed original image, Top right: Zoomed image reconstructed by the DWT (PSNR = 25.11), Bottom left: Zoomed image reconstructed by the CT (PSNR = 25.88), Bottom right: Zoomed image reconstructed by the DST with only 1 iteration step (PSNR = 26.73).


S Dahlke
email: dahlke@mathematik.uni-marburg.de

G Steidl
email: steidl@math.uni-mannheim.de

G Teschke
email: teschke@hs-nb.de

The Continuous Shearlet Transform in Arbitrary Space Dimensions, Frame Construction, and Analysis of Singularities

This note is concerned with the generalization of the continuous shearlet transform to higher dimensions. Similar to the two-dimensional case, our approach is based on translations, anisotropic dilations and specific shear matrices. We show that the associated integral transform again originates from a square-integrable representation of a specific group, the full n-variate shearlet group. Moreover, we verify that by applying the coorbit theory, canonical scales of smoothness spaces and associated Banach frames can be derived. We also indicate how our transform can be used to characterize singularities in signals.

Introduction

Modern technology allows for easy creation, transmission and storage of huge amounts of data. Confronted with a flood of data, such as internet traffic, or audio and video applications, nowadays the key problem is to extract the relevant information from these sets. To this end, usually the first step is to decompose the signal with respect to suitable building blocks which are well-suited for the specific application and allow a fast and efficient extraction. In this context, one particular problem which is currently in the center of interest is the analysis of directional information. Due to the bias to the coordinate axes, classical approaches such as, e.g., wavelet or Gabor transforms are clearly not the best choices, and hence new building blocks have to be developed. In recent studies, several approaches have been suggested such as ridgelets [START_REF] Candès | Ridgelets: a key to higher-dimensional intermittency?[END_REF], curvelets [START_REF] Candès | Curvelets -A surprisingly effective nonadaptive representation for objects with edges[END_REF], contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF], shearlets [START_REF] Guo | Sparse multidimensional representations using anisotropic dilation und shear operators[END_REF] and many others. For a general approach see also [START_REF] Guo | Wavelets with composite dilations and their MRA properties[END_REF]. Among all these approaches, the shearlet transform stands out because it is related to group theory, i.e., this transform can be derived from a square-integrable representation π : S → U(L 2 (R 2 )) of a certain group S, the socalled shearlet group, see [START_REF] Dahlke | The uncertainty principle associated with the continuous shearlet transform[END_REF]. Therefore, in the context of the shearlet transform, all the powerful tools of group representation theory can be exploited. So far, the shearlet transform is well developed for problems in R 2 . However, for analyzing higher-dimensional data sets, there is clearly an urgent need for further generalizations and applications. This is exactly the concern of this paper. One particular field of application is the geometrical structure analysis of multi-dimensional data, e.g. multimodal spectral measurements in meteorology. To our best knowledge, it seems that there exist only few results in this direction: some important progress has been achieved for the curvelet case in [START_REF] Borup | Frame decomposition of decomposition spaces[END_REF] and for surfacelets in [START_REF] Lu | Multidimensional directional filterbanks and surfacelets[END_REF]. However, for the shearlet approach the question has been completely open.

Multivariate Continuous Shearlet Transform

In this section, we introduce the shearlet transform on L 2 (R n ). This requires the generalization of the twodimensional parabolic dilation matrix and of the shear matrix, respectively. Let I n denote the (n, n)-identity matrix and 0 n , resp. 1 n the vectors with n entries 0, resp. 1. For a ∈ R * := R \ {0} and s ∈ R n-1 , we set

A a := a 0 T n-1 0 n-1 sgn(a)|a| 1 n I n-1
and In the following, we use only the left Haar measure and use the abbreviation

S s := 1 s T 0 n-1 I n-1 . Lemma 1 The set R * × R n-1 × R n endowed
dµ = dµ l . For f ∈ L 2 (R n ) we define π(a, s, t)f (x) = f a,s,t (x) := |a| 1 2n -1 f (A -1 a S -1 s (x -t)). (1) It is easy to check that π : S → U(L 2 (R n )) is a mapping from S into the group U(L 2 (R n )) of unitary operators on L 2 (R n ).
Recall that a unitary representation of a locally compact group G with the left Haar measure µ on a Hilbert space H is a homomorphism π from G into the group of unitary operators U(H) on H which is continuous with respect to the strong operator topology.

Lemma 2

The mapping π defined by (1) is a unitary representation of S.

A nontrivial function ψ ∈ L 2 (R n ) is called admissible, if S | ψ, π(a, s, t)ψ | 2 dµ(a, s, t) < ∞.
If π is irreducible and there exits at least one admissible

function ψ ∈ L 2 (R n ), then π is called square integrable.
The following result shows that the unitary representation π defined in ( 1) is square integrable.

Theorem 3 A function ψ ∈ L 2 (R n ) is

admissible if and only if it fulfills the admissibility condition

C ψ := R n | ψ(ω)| 2 |ω 1 | n dω < ∞. (2)
Then, for any f ∈ L 2 (R n ), the following equality holds true:

S | f, ψ a,s,t | 2 dµ(a, s, t) = C ψ f 2 L2(R n ) . (3)
In particular, the unitary representation π is irreducible and hence square integrable.

An example of a continuous shearlet can be constructed as follows: Let

ψ 1 be a continuous wavelet with ψ1 ∈ C ∞ (R) and supp ψ1 ⊆ [-2, -1 2 ] ∪ [ 1 2 , 2], and let ψ 2 be such that ψ2 ∈ C ∞ (R n-1 ) and supp ψ2 ⊆ [-1, 1] n-1 . Then the function ψ ∈ L 2 (R n ) defined by ψ(ω) = ψ(ω 1 , ω) = ψ1 (ω 1 ) ψ2 1 ω 1 ω
is a continuous shearlet. The support of ψ is depiced for ω 1 ≥ 0 in Fig. 1.

Multivariate Shearlet Coorbit Theory

In this section we want to establish a coorbit theory based on the square integrable representation (1) of the shearlet group. We mainly follow the lines of [START_REF] Dahlke | Shearlet Coorbit Spaces and Associated Banach Frames[END_REF]. For further information on coorbit space theory, the reader is referred to [START_REF] Feichtinger | A unified approach to atomic decompositions via integrable group representations[END_REF][START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decomposition I[END_REF][START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decomposition II[END_REF][START_REF] Feichtinger | Nonorthogonal wavelet and Gabor expansions and group representations[END_REF][START_REF] Gröchenig | Describing functions: Atomic decompositions versus frames[END_REF]].

E T © f f f f f f f 1 2 2 ω 1 ω 3 ω 2
Figure 1: Support of the shearlet ψ for ω 1 ≥ 0.

Shearlet Coorbit Spaces

We consider weight functions w(a, s, t) = w(a, s) that are locally integrable with respect to a and s, i.e., w ∈ L loc 1 (R n ) and fulfill w ((a, s, t) • (a , s , t )) ≤ w(a, s, t)w(a , s , t ) and w(a, s, t) ≥ 1 for all (a, s, t), (a , s , t ) ∈ S. For 1 ≤ p < ∞, let L p,w (S) := {F measurable : F Lp,w(S) := S |F (g)| p w(a, s, t) p dµ(a, s, t)

1 p < ∞},
and let L ∞,w be defined with the usual modifications. In order to construct the coorbit spaces related to the shearlet group we have to ensure that there exists a function

ψ ∈ L 2 (R n ) such that SH ψ (ψ) = ψ, π(a, s, t)ψ ∈ L 1,w (S). (4)
Fortunately, it turns out that (4) can be satisfied in our setting.

Theorem 4 Let ψ be a Schwartz function such that

supp ψ ⊆ ([-a 1 , -a 0 ] ∪ [a 0 , a 1 ]) × Q b ,where Q b := [-b 1 , b 1 ] × • • • × [-b n-1 , b n-1 ]. Then we have that SH ψ (ψ) ∈ L 1,w (S), i.e., ψ, π(•)ψ L1,w(S) = S |SH ψ (ψ)(a, s, t)| w(a, s, t) dµ(a, s, t) < ∞.
For ψ satisfying (4) we can consider the space

H 1,w := {f ∈ L 2 (R n ) : SH ψ (f ) ∈ L 1,w (S)}, (5)
with norm f H1,w := SH ψ f L1,w(S) and its antidual H ∼ 1,w , the space of all continuous conjugate-linear functionals on H 1,w . The spaces H 1,w and H ∼ 1,w are π-invariant Banach spaces with continuous embeddings H 1,w → H → H ∼ 1,w , and their definition is independent of the shearlet ψ. Then the inner product on L 2 (R n ) × L 2 (R n ) extends to a sesquilinear form on H ∼ 1,w × H 1,w , therefore for ψ ∈ H 1,w and f ∈ H ∼ 1,w the extended representation coefficients SH ψ (f )(a, s, t) := f, π(a, s, t)ψ H ∼ 1,w ×H1,w are well-defined. Now, for 1 ≤ p ≤ ∞, we define the shearlet coorbit spaces

SC p,w := {f ∈ H ∼ 1,w : SH ψ (f ) ∈ L p,w (S)} (6)
with norms f SCp,w := SH ψ f Lp,w(S) . It holds that SC 1,w = H 1,w and SC 1,1 = L 2 (R n ).

Shearlet Banach Frames

The Feichtinger-Gröchenig theory provides us with a machinery to construct atomic decompositions and Banach frames for our shearlet coorbit spaces SC p,w . In a first step, we have to determine, for a compact neighborhood U of e ∈ S with non-void interior, so-called U -dense sets. A (countable) family X = ((a, s, t) λ ) λ∈Λ in S is said to be U -dense if ∪ λ∈Λ (a, s, t) λ U = S, and separated if for some compact neighborhood Q of e we have (a i , s i , t i )Q ∩ (a j , s j , t j )Q = ∅, i = j, and relatively separated if X is a finite union of separated sets.

Lemma 5 Let U be a neighborhood of the identity in S, and let α > 1 and β, γ > 0 be defined such that

[α 1 n -1 , α 1 n ) × [-β 2 , β 2 ) n-1 × [-γ 2 , γ 2 ) n ⊆ U. (7) Then the sequence {( α j , βα j(1-1 n ) k, S βα j(1-1 n ) k A α j γm) : j ∈ Z, k ∈ Z n-1 , m ∈ Z n , ∈ {-1, 1}} (8)
is U -dense and relatively separated.

Next we define the U -oscillation as

osc U (a, s, t) := sup u∈U |SH ψ (ψ)(u • (a, s, t)) -SH ψ (ψ)(a, s, t)|. (9)
Then, the following decomposition theorem, which was proved in a general setting in [START_REF] Feichtinger | A unified approach to atomic decompositions via integrable group representations[END_REF][START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decomposition I[END_REF][START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decomposition II[END_REF][START_REF] Feichtinger | Nonorthogonal wavelet and Gabor expansions and group representations[END_REF][START_REF] Gröchenig | Describing functions: Atomic decompositions versus frames[END_REF], says that discretizing the representation by means of an U -dense set produces an atomic decomposition for SC p,w .

Theorem 6 Assume that the irreducible, unitary representation π is w-integrable and let an appropriately nor-

malized ψ ∈ L 2 (R n ) which fulfills M ψ, π(a, s, t) := sup u∈(a,s,t)U | ψ, π(u)ψ | ∈ L 1,w (S) (10
) be given. Choose a neighborhood U of e so small that osc U L1,w(S) < 1.

Then for any U -dense and relatively separated set X = ((a, s, t) λ ) λ∈Λ the space SC p,w has the following atomic decomposition:

If f ∈ SC p,w , then f = λ∈Λ c λ (f )π((a, s, t) λ )ψ (12)
where the sequence of coefficients depends linearly on f and satisfies

(c λ (f )) λ∈Λ p,w ≤ C f SCp,w (13)
with a constant C depending only on ψ and with p,w being defined by p,w := {c = (c λ ) λ∈Λ : c p,w := cw p < ∞},

where w = (w((a, s, t) λ )) λ∈Λ . Con- versely, if (c λ (f )) λ∈Λ ∈ p,w , then f = λ∈Λ c λ π((a, s, t) λ )ψ is in SC p,w and f SCp,w ≤ C (c λ (f )) λ∈Λ p,w . (14)
Given such an atomic decomposition, the problem arises under which conditions a function f is completely determined by its moments f, π((a, s, t) λ )ψ and how f can be reconstructed from these moments. This is answered by the following theorem which establishes the existence of Banach frames.

Theorem 7 Impose the same assumptions as in Theorem 6. Choose a neighborhood U of e such that osc U L1,w(S) < 1/ SH ψ (ψ) L1,w(S) .

Then, for every U -dense and relatively separated family X = ((a, s, t) λ ) λ∈Λ in G the set {π((a, s, t) λ )ψ : λ ∈ Λ} is a Banach frame for SH p,w . This means that

i) f ∈ SC p,w if and only if ( f, π((a, s, t) λ )ψ H ∼ 1,w ×H1,w ) λ∈Λ ∈ p,w ; ii) there exist two constants 0 < D ≤ D < ∞ such that D f SCp,w ≤ ( f, π((a, s, t) λ )ψ H ∼ 1,w ×H1,w ) λ∈Λ p,w ≤ D f SCp,w ; (16)
iii) there exists a bounded, linear reconstruction operator R from p,w to SC p,w such that

R ( f, ψ((a, s, t) λ )ψ H ∼ 1,w ×H1,w ) λ∈Λ = f.
It can be checked that the conditions ( 10), ( 11) and ( 15) can be satisfied, see [START_REF] Dahlke | The continuous shearlet transform in arbitrary space dimensions[END_REF] for details.

Analysis of Singularities

In this section, we deal with the decay of the shearlet transform at hyperplane singularities. The following analysis generalizes techniques and results presented in [START_REF] Kutyniok | Resolution of the wavefront set using continuous shearlets[END_REF] for two dimensions. An (n -m)-dimensional hyperplane in R n , 1 ≤ m ≤ n -1, not containing the x 1 -axis can be written w.l.o.g. as

   x 1 . . . x m    x A + P    x m+1 . . . x n    x E =    0 . . . 0    , P :=    p T 1 . . . p T m    ∈ R m,n-m .
Then we obtain for

ν m := δ(x A + P x E )
with the Delta distribution δ that νm (ω) = R n δ(x A + P x E )e -2πi( x A ,ω A + x E ,ω E ) dx = R n-m e -2πi(-P x E ,ω A + x E ,ω E ) dx E = δ(ω E -P T ω A ).

The following theorem describes the decay of the shearlet transform at hyperplane singularities. We use the notation SH ψ f (a, s, t) ∼ |a| r as a → 0, if there exist constants 0 < c ≤ C < ∞ such that c|a| r ≤ SH ψ f (a, s, t) ≤ C|a| r as a → 0. as a → 0.

(18)

Otherwise, the shearlet transform SH ψ ν m decays rapidly as a → 0.

Similar results can be derived for point singularities, see again [START_REF] Dahlke | The continuous shearlet transform in arbitrary space dimensions[END_REF] for details.

Theorem 8

 8 Let ψ ∈ L 2 (R n ) be a shearlet satisfying ψ ∈ C ∞ (R n ). Assume further that ψ(ω) = ψ1 (ω 1 ) ψ2 (ω/ω 1 ), where supp ψ1 ∈ [-a 1 , -a 0 ] ∪ [a 0 , a 1 ] for some a 1 > a 0 ≥ α > 0 and supp ψ2 ∈ Q b . If (s m , . . . , s n-1 ) = (-1, s 1 , . . . ,s m-1 ) P and (t 1 , . . . , t m ) = -(t m+1 , . . . , t n ) P T , then SH ψ ν m (a, s, t) ∼ |a| 1-2m 2n


Masaru Kamada
email: kamada@mx.ibaraki.ac.jp

Semih Özlem

Hiromasa Habuchi
email: habuchi@mx.ibaraki.ac.jp

Orthogonal Exponential Spline Pulses with Application to Impulse Radio

With application to the impulse radio communications in mind, a locally supported and zero-mean pulse which is orthogonal to its shifts by integers is sought among the exponential splines having the knot interval 1

2 . An example pulse is obtained that complies with the regulation imposed by the US Federal Communications Commission and will potentially enable an impulse radio communications system as fast as 6G pulses per second.

Introduction

The M-shaped linear spline

M (t) =            √ 3t, 0 ≤ t ≤ 1 2 √ 3(2 -3t), 1 2 ≤ t ≤ 1 √ 3(3t -4), 1 ≤ t ≤ 3 2 √ 3(2 -t), 3 2 ≤ t ≤ 2 0, elsewhere (1)
plotted in Fig. 1 is not a wavelet in the sense of muntiresolutional analysis because M (t) is not orthogonal to its contracted version M (2t). But it has three remarkable properties that (i) it is locally supported, (ii) its integration over the domain is zero, and (iii) its shifts by integers are orthogonal to one another [START_REF] Jerri | Wavelets -Detailed Treatment with Applications[END_REF]. Those properties are exactly what is required of pulses for the impulse radio communications [START_REF] Win | Impulse radio: how it works[END_REF]. The three properties are required (i) for the sake of real-time communications, (ii) for the pulse to be feasible as a radio waveform, and (iii) for pulse detection to be robust against noise in the sense of least-square estimation, respectively.

We shall look for this kind of pulse functions in the broader family of exponential splines [START_REF] Unser | Cardinal exponential splines: Part I-Theory and filtering algorithms[END_REF][START_REF] Unser | Cardinal exponential splines: Part II-Think analog, act digital[END_REF] which have the

1 -1 0 M(t) 1 2
Figure 1: M-shaped linear spline.

advantage that they can be shaped through linear dynamical systems [START_REF] Unser | Cardinal exponential splines: Part II-Think analog, act digital[END_REF] . The pulse functions, if they are found, will work as practical pulses which carry information in the impulse radio communications.

The problem is simple: we are to find a locally supported and zero-mean exponential spline q(t) with the knot inter-

val 1 2 that satisfies ∫ ∞ -∞ q(t)q(t -k)dt = { 1, k = 0 0, k ̸ = 0 (2)
for any ingeter k. This paper presents a procedure to find such a pulse function and its application to the impulse radio.

Construction of orthogonal pulses

Any exponential spline can be represented by a linear combination of the exponential B-spline and its shifts [START_REF] Unser | Cardinal exponential splines: Part I-Theory and filtering algorithms[END_REF][START_REF] Unser | Cardinal exponential splines: Part II-Think analog, act digital[END_REF]. An exponential B-spline with the knot interval 1 2 is the output

β(t) = S(b)(t) (3)
of a linear dynamical system S having the transfer function

G(s) = µ n-1 s n-1 + • • • + µ 1 s + µ 0 (s -λ 0 )(s -λ 1 ) • • • (s -λ n-1 ) (4)
for the input being a series of delta functions

b(t) = n ∑ l=0 b l δ(t -l/2) ( 5
)
such that

B(z) = n ∑ l=0 b l z -l 2 = (1-z -1 2 e λ 0 2 )(1-z -1 2 e λ 1 2 ) • • • (1-z -1 2 e λ n-1 2
). [START_REF] Win | Impulse radio: how it works[END_REF] This exponential B-spline is locally supported as

β(t) = 0, t / ∈ ( 0, n 2 ) . ( 7
)
In order to keep the splines zero-mean, instead of the original exponential B-spline β(t), we shall use

α(t) = β(t) -β ( t - 1 2
) (8) which has the zero mean ∫ ∞ -∞ α(t)dt = 0 (9)
and is locally supported as

α(t) = 0, t / ∈ ( 0, n + 1 2 ) . ( 10
)
Another representation of this α(t) is the output

α(t) = S(a)(t) (11)
of S for the input

a(t) = n ∑ l=0 a l δ(t -l/2), ( 12
)
where

A(z) = n+1 ∑ l=0 a l z -l 2 = (1-z -1 2 e λ 0 2 ) • • • (1-z -1 2 e λ n-1 2 )(1 -z -1 2 ).( 13
)
Let the desired pulse function be represented in the form

q(t) = n-1 ∑ l=0 c l α(t -l/2). ( 14
)
Then it is automatic that q(t) is locally supported as

q(t) = 0, t / ∈ (0, n) (15)
and has the zero mean

∫ ∞ -∞ q(t)dt = 0. ( 16
)
The remaining request is that its autocorrelation

r(x) = ∫ ∞ -∞ q(t)q(t -x)dt (17)
should satisfy the orthogonality conditions

r(k) = { 1, k = 0 0, k = ±1, ±2, • • • (18)
with respect to shift by integers. Here the number n of {α(t -l/2)} n-1 l=0 employed for composing q(t) in ( 14) is chosen so that the number n of the unknown coefficients {c l } n-1 l=0 be the same as that of the essential conditions

r(k) = { 1, k = 0 0, k = 1, 2, • • • , n -1 (19)
reduced from (18) by ( 15) and the equality r(x) = r(-x).

Now we have only to find the coefficients {c l } n-1 l=0 that make (19) hold good. Define

c(t) = n-1 ∑ l=0 c l δ(t -l/2) and C(z) = n-1 ∑ l=0 c l z -l 2 (20) by {c l } n-1 l=0
, and prepare time-reversed functions

ã(t) = a(-t), c(t) = c(-t), q(t) = q(-t) (21)
and the "mirror" system S having the transfer function G(-s). Then we can express the correlation by

r(k) = (q * q)(k) = (S • S)(a * ã * c * c)(k), (22)
where * denotes the convolution integral, and we can write

D(z) = C(z)C(z -1 ) in the form C(z)C(z -1 ) = d 0 + n-1 ∑ j=1 d j (z -j 2 + z j 2 ) ( 23
)
which implies

(c * c)(t) = d 0 δ(t)+ n-1 ∑ j=1 d j (δ(t-j/2) + δ(t+j/2)) .(24)
In the meantime, a locally supported exponential spline

ϕ(x) = (S • S)(a * ã)(x) (25)
associated with the composite system S • S satisfies

ϕ(x) = ϕ(-x). ( 26
)
By ( 22), ( 24), ( 25) and ( 26), we can reduce the orthogonality conditions (19) to the linear equations

d 0 ϕ(k)+ n-1 ∑ j=1 d j (ϕ(k -j/2) + ϕ(k + j/2)) = { 1, k = 0 0, k = 1, 2, • • • , n -1. (27)
Solvability of (27) for {d j } n-1 j=0 can be checked by numerical computation in practice. A simpler condition in terms of dynamical parameters is yet to be established.

We assume that ( 27) is solvable since we cannot proceed unless this is the case. Then, C(z)C(z -1 ) determined by (23) from {d j } n-1 j=0 can be factorized in the form

C(z)C(z -1 ) = γ 0 (z -1 2 -γ 1 )(z 1 2 -γ 1 )(z -1 2 -γ 2 )(z 1 2 -γ 2 ) • • • (z -1 2 -γ n-1 )(z 1 2 -γ n-1 ). ( 28
)
Taking half the factors, we can find

C(z) = ± √ γ 0 (z -1 2 -γ 1 )(z -1 2 -γ 2 ) • • • (z -1 2 -γ n-1 ) (29)
that gives the sought coefficients {c l } n-1 l=0 by (20). Exciting the system S with the input series of delta functions

v(t) = n-1 ∑ l=0 c l a(t -l/2), ( 30
)
we obtain the desired pulse function

q(t) = S(v)(t) = n-1 ∑ l=0 c l α(t -l/2). ( 31
)
In the case G(s) = 1 s , the problem is trivial and the resulting pulse is the Haar function

H(t) =    1, 0 < t ≤ 1 2 -1, 1 2 < t ≤< 1 0, elsewhere. (32)
The case G(s) = 1 s 2 yields M (t) of (1) as expected. Because it happens that M (t) = √ 3(H * H)(t), we might speculate that the pulse associated with G(s) = 1 s 3 could be proportional to (H * H * H)(t). But that is not true since

(H * H * H)(t) is not orthogonal to (H * H * H)(t -2).
It is interesting as well as disappointing that we obtain a complex-valued pulse in the case G(s) = 1 s 3 . A nice example pulse will appear in the next section in the context of its application to the impulse radio communications.

Application to Impulse Radio

While the series of delta functions a(t) does not exist in the real world, its integration

∫ t -∞ a(τ )dτ =    0, t < 0 ∑ l k=0 a k , l 2 < t < l+1 2 , l = 0, 1, • • • , n ∑ n+1 k=0 a k = A(1) = 0, n+1 2 < t (33)
is a locally supported piecewise constant function that can be easily generated by electric current switches.

The system S excited by the piecewise constant function

u(t) = ∫ t -∞ v(τ )dτ = n ∑ l=0 c l ∫ t-l/2 -∞ a(τ )dτ (34)
shapes the pulse

p(t) = S(u)(t) (35)
which is locally supported as

p(t) = 0, t / ∈ (0, n) (36)
and has the relationship

p(t) = n ∑ l=0 c l ∫ t-l/2 -∞ α(τ )dτ = ∫ t -∞ q(τ )dτ. (37)
Besides the simple and practical system (35) to shape p(t) from the piecewise constant seed u(t), the pulse p(t) has the remarkable property

∫ ∞ -∞ d 2 dt 2 p(t) p(t -k)dt = - ∫ ∞ -∞ q(t)q(t -k)dt = { -1, k = 0 0, k = ±1, ±2, • • • (38)
which follows from (17), ( 18), (36), (37) and the partial integration formula. This property gives the foundation to transmission and detection of the pulse p(t) in the impulse radio communications.

Given data bits {w l }, we transmit the waveform

w(t) = S ( ∞ ∑ l=-∞ w l u(t -l) ) = ∞ ∑ l=-∞ w l p(t -l) (39)
as illustrated in Fig. 2. Since a good broadband antenna is well approximated [START_REF] Win | Impulse radio: how it works[END_REF] by d dt , the transmitted signal w(t) is differentiated once by the transmitter antenna to be the radio signal

d dt w(t) = ∞ ∑ l=-∞ w l d dt p(t -l) (40)
and again by the receiver antenna to arrive at the receiver as

d 2 dt 2 w(t) = ∞ ∑ l=-∞ w l d 2 dt 2 p(t -l). ( 41
)
Correlating the received signal d 2 dt 2 w(t) with the template pulse p(t-k), which is the same as the transmission pulse, for its duration (k, k + n), we have the bit w k recovered by

∫ k+n k d 2 dt 2 w(t) p(t-k)dt = ∫ ∞ -∞ d 2 dt 2 w(t) p(t-k)dt = ∞ ∑ l=-∞ w l ∫ ∞ -∞ d 2 dt 2 p(t-l) p(t-k)dt = -w k (42)
because of the property (38).

It should be noted that, because of (38), the detection formula (42) virtually performs the least-squares approximation of the radio waveform d dt w(t) by d dt p(t-k) = q(t-k) to detect w k . Additive noises superimposed on d dt w(t) will then be most suppressed in the sense of least-squares estimation.

An example pulse associated with the transfer function

G(s) = 1 (s+18)(s+11.1i+10 -13 )(s-11.1i+10 -13 ) (43)
and its derivatives are plotted in Fig. 3. The correlation in Fig. 4 becomes 1 and 0 at the origin and at the other integers, respectively, to verify (38). The power spectral density of the radio pulse d dt p(t) = q(t) is plotted in Fig. 5 along with the spectral mask (plotted by the boxy line) for the indoor ultra-wideband communications systems [START_REF]Revision of Part 15 the Commission's rule regarding ultra-wideband transmission systems[END_REF] imposed by the US Federal Communications Commission as the upper bound which no practical pulses are allowed to exceed. The frequency axis of the mask is scaled down by 6 GHz for the purpose of comparison, or equivalently, the pulse repetition rate is assumed to be 6 G pulses per second, which is much faster than the 1.32G pulses per second of the high speed direct sequence ultra-wideband protocol discussed in the IEEE 802.15.3a standard.

The fast transmission is possible because the pulses are orthogonal even though they are densely overlapping. But dense pulses are prone to interfere with one another in the situation that several reflected pulses arrive with various delays. Multipath compensation by digital filtering is crucial in order to effectively exploit the dense pulses we obtained. Transmitting a sounder pulse and digitizing the observed correlations, we have the end-to-end impulse response of the multipath channel. Digital filtering by an FIR approximation of the inverse impulse response will work as a kind of rake receiver. This compensation requires an analog-to-digital converter and a digital filter that work at the pulse rate and thus costs more hardware. But this cost should be justified since all the pulse-based systems cannot be faster without having denser pulses in the first place. A detailed analysis of the multipath effects, channel modeling error, and pulse synchronization is available in [START_REF] Özlem | Construction of orthogonal overlapping pulses for impulse radio communications[END_REF].

We may ignore the multipath effects and channel modeling error in the extreme situation that antennas are inductively coupled at a very short distance less than one inch. TransferJet technology has been working in the same situation at the maximum transmission rate of 560Mbps since 2008. A faster system will hopefully be the first application of the dense pulses obtained in this paper.

Conclusions

Inspired by the M-shaped orthogonal pulse, we derived a procedure to construct an exponential spline pulse with the knot interval 1 2 that is locally supported, has its mean zero, and is orthogonal to its shifts by integers. An exam- ple pulse was obtained that will potentially enable an impulse radio communications system as fast as 6G pulses per second under the FCC regulation for the indoor ultracommunications.

Acknowledgment

This work was partially supported by JSPS grant-in-aid No. 17560357.

Figure 2 :

 2 Figure 2: Schematic diagram of the transceiver.

Figure 3 :

 3 Figure 3: Pulses for impulse radio.

Figure 4 :

 4 Figure 4: Correlation of the pulse.

Figure 5 :

 5 Figure 5: Power spectral density of the pulse and the FCC spectral mask.