Non-uniform sampling problem S.H. Kulkarni, R. Radha and S. Sivananthan Department of Mathematics, Indian Institute of Technology Madras, Chennai - 600 036, India. e-mail:
[email protected],
[email protected], ssiva
[email protected] Abstract If ϕ is a continuous function satisfying certain decay properties, then it is shown that any f belongs to the shift invariant space V (ϕ) can be reconstructed uniquely and stably from its sample points which are perturbed version of integers. The reconstruction procedure is discussed based on the convergence of the projection method of Laurent operators and the computation of the solutions of finite dimensional truncated matrix equations recursively. The implementation is carried out using MatLab. Some illustrations are given.
1
Introduction
A sampling problem involves the reconstruction of a function f from its samples f (xj ), j ∈ Z. Unless apriori some conditions are assumed on f , the problem is not meaningful. In fact, one has to obtain some conditions on a sampling set X = {xj ∈ R/j ∈ J}, J a countable set, to recover a certain class of functions f from the sample set X in a stable and unique manner. Since sampling problem deals with real world applications, it is also important to find a stable numerical algorithm to obtain a reconstruction procedure. The classical Shannon sampling theorem states that if f ∈ L2 (R) ∩ B[− 1 , 1 ] , then f 2 2 can be reconstructed from its samples {f (k)/k ∈ Z} by the formula X f (x) = f (k) sinc(x − k), k∈Z
where sinc(y) = sinπyπy and B[− 1 , 1 ] denotes the class of functions in L1 (R) whose Fourier 2 2 transforms have support in [− 12 , 21 ]. Kadec in [10] proved that if X = {xk ∈ R : |xk − k| ≤ L < 14 , k ∈ Z} then any f ∈ L2 (R) ∩ B[− 1 , 1 ] can be completely recovered from the sample 2 2 points in X. Another important result in this direction is stated in terms of Beruling density. Let #[X ∩ (y + [0, r])] D(X) = lim inf . r→∞ y∈R r If D(X) > 1, then any f ∈ L2 ∩B[− 1 , 1 ] can be reconstructed from its samples {f (xj )/xj ∈ 2 2 X} stably and uniquely. Conversely, if every f ∈ L2 ∩ B[− 1 , 1 ] can be uniquely and stably 2 2
1
reconstructed from its samples then D(X) ≥ 1 (see [11]). Notice that in all these cases f is assumed to be a band limited function, that is, its Fourier transform has a compact support. The problem of nonuniform sampling in general shift invariant spaces was studied by Aldroubi and Gr¨ochenig in [1]. In this paper they proved that if f belongs to a shift invariant space V (ϕ), where ϕ is an nth order B-spline function, then f can be reconstructed stably and uniquely from a non-uniform sampled data. They have also obtained equivalent conditions for the stable set of sampling for V (ϕ) for certain ϕ (see proposition 1 in [1]). A more detailed study of shift invariant spaces can be found in [2]. In this paper, we prove that if ϕ belongs to a more general class namely ϕ ∈ F1 given in definition 3.3, and if X is a non-uniform sample consisting of perturbed version of integers, then f ∈ V (ϕ) can be reconstructed uniquely and stably. Before stating other main results of this paper, we recall the following definitions. Definition 1.1. The Wiener amalgam space W (Lp ) (1 ≤ p < ∞) is defined to be the class of measurable functions f satisfying X ess sup{|f (x + k)|p /x ∈ [0, 1]} < ∞. || f ||pW (Lp ) := k∈Z
For the importance of these spaces and further details, we refer to [6] and [7]. Let Wc (Lp ) denote the class of functions in W (Lp ) which are continuous. Definition 1.2. A sampling set X = {xj ∈ R/j ∈ J} is said to be a stable set of sampling for a class of functions in V ⊆ L2 (R) if there exists A, B > 0 such that X A||f ||22 ≤ |f (xj )|2 ≤ B||f ||22 ∀ f ∈ V. (1.1) xj ∈X
Definition 1.3. Shift-invariant space V (ϕ): Let ϕ be a continuous real valued function on R. Assume that ∃ A0 , B 0 > 0 satisfying X 0 0 A ||c||`2 (Z) ≤ || ck ϕ(· − k)||2 ≤ B ||c||`2 (Z) , ∀ c = (ck ) ∈ `2 (Z) (1.2) k∈Z
Let τk ϕ be defined by τk ϕ(x) = ϕ(x − k), where k ∈ Z and X V (ϕ) := {f = ck τk ϕ / (ck ) ∈ `2 (Z)}. k∈Z
Notice that V (ϕ) is (integer)shift invariant, this means that if f ∈ V (ϕ), then τk f ∈ V (ϕ), ∀ k ∈ Z. The following theorem is proved in [1]. 2
Theorem 1.4. Let ϕ ∈ Wc (L1 ) and ϕ satisfy (1.2). Then V (ϕ) ⊆ Wc (L2 ) and the following conditions are equivalent: (a) X is a stable set of sampling for V (ϕ) (b) There exists constants A1 , A2 > 0 such that A1 || c ||`2 (Z) ≤ || U c ||`2 (Z) ≤ A2 || c ||`2 (Z) ,
∀ c ∈ `2 (Z)
where U is an infinite matrix with entries Ujk = ϕ(xj − k), j, k ∈ Z. The condition (1.1) shows that the reconstruction is unique. The inequalities in theorem 1.4(b) show that the operator U is bounded and is also bounded below. Obviously U satisfies these conditions if we take ϕ(x) = sinπxπx and X is the class of integers or the perturbed version of integers. In fact this is the reformulation of the statements of Shannon’s theorem and Kadec’s theorem. To the best our knowledge, it has been shown that U is bounded and is also bounded below only when ϕ = sinπxπx or an nth order B-spline. In this paper, we show that if ϕ ∈ F1 and X is the perturbed version of integers, it turns out that U is invertible. This leads to the stable and unique reconstruction for all f ∈ V (ϕ). We then give a reconstruction procedure and illustrate it with some examples using Matlab. The invertibility and projection method of Laurent operator along with recurrence relation between the solutions of the matrix equation of the finite dimensional truncations of an operator help us to obtain an iterative reconstruction procedure. The questions of speed, storage, complexity of the algorithms are not investigated and hence no claims are made in this regard. We organize our paper as follows. In section 2, we discuss the invertibility of a Laurent operator and convergence of projection method for a Laurent operator. In section 3, we prove our main results. In section 4, we prove a recurrence relation between the solutions of the matrix equation of the finite dimensional truncations of an operator and discuss the reconstruction procedure. In section 5, we provide algorithm along with illustration.
2
Laurent Operators
In this section we briefly review some properties of Laurent operators and projection method applied to operator equation involving Laurent operators. A detailed treatment of these topics including proofs of theorems 2.2, 2.3, 2.4, 2.5 can be found in [9]. Definition 2.1. A bounded operator A on `2 (Z) is a Laurent operator if < Aek , ej > depends on the difference j − k only, where {en }∞ n=−∞ is the standard orthonormal basis 2 of ` (Z). 3
Let ‘a’ be a bounded complex valued function on [0, 1]. Then Z 1 a(t)e2πi(m−n)t dt = am−n , ∀ m, n ∈ Z amn = 0
which defines a Laurent operator A on `2 (Z). In this case we say that A is the Laurent operator defined by ‘a’. The function ‘a’ is called the symbol of A. The properties of the Laurent operator A and its symbol ‘a’ are closely related as shown by the following theorem. Theorem 2.2. Let A be the Laurent operator defined by the continuous function ‘a’. Then A is invertible iff a(t) 6= 0 for each 0 ≤ t ≤ 1. Theorem 2.3. Let T be a bounded and invertible operator on `2 (Z). If S is a bounded operator on `2 (Z) and ||S|| < ||T 1−1 || then T + S is invertible. Let {Pn } denote a sequence of projection in B(H) (the class of bounded operators on a separable Hilbert space H) such that Pn → I in the strong operator topology. Given T ∈ B(H), consider Tn = Pn T Pn . The method of obtaining solution x of the operator equation T x = y as a limit (in the norm topology) of xn satisfying Tn xn = Pn (y) is called a projection method. The projection method for T is said to converge if there exists an integer N0 such that for each y ∈ H and n ≥ N0 , there exists a unique solution xn satisfying Tn xn = Pn (y) and in addition the sequence {xn } converges to T −1 y. Let Π(Pn ) denote the class of invertible operators for which the projection method converges. Let H denote the separable Hilbert space having an orthonormal basis {..., e−n , e−n−1 , ..., e−1 , e0 , e1 , ..., en−1 , en , ...}. Let Hn denote the closed subspace of H generated by {e−n , ..., e0 , ..., en }. We shall need the following known theorems for proving the main results in the next section. Theorem 2.4. Suppose T ∈ Π(Pn ) and M = supn≥N (||(Pn T Pn )−1 || +||Pn ||). If S ∈ B(H) such that || S || < min{ ||T 1−1 || , M12 }, then T + S ∈ Π(Pn ). Theorem 2.5. Let T be an invertible Laurent operator with continuous symbol ω. Then the projection method converges for T iff the winding number of ω relative to zero is equal to zero.
3
The Main Results
Definition 3.1. Let F denote the class of all complex valued continuous functions ϕ defined on R satisfying the following conditions (i) ∃ c1 , c2 > 0 and α > 1 such that c2 c1 , | ϕ(x)| ˆ ≤ |ϕ(x)| ≤ 1 + |x|α 1 + |x|α where ϕˆ denotes the Fourier transform of ϕ. 4
, ∀ x ∈ R,
(ii) ϕ(x) ˆ 6= 0 on some unit interval [b, b + 1], for some b ∈ R. Example 3.2. It is easy to show that x2
• the class of Hermite functions hn where hn (x) = (−1)n e 2 2
2 dn (e−x ), dxn
n = 0, 1, ...
− x2
with h0 = e
• the function ϕ2n = χ[a,a+1] ∗ χ[a,a+1] ∗ · · · ∗ χ[a,a+1] , 2n − 1 fold convolution of a characteristic function on [a, a + 1], where a ∈ Z2 and n ≥ 1 are contained in F. Proposition 3.3. Let ϕ ∈ F. Then V (ϕ) is a closed subspace of L2 (R), {τk ϕ/k ∈ Z} is a Riesz basis for V (ϕ) and V (ϕ) ⊆ Wc (L2 ). P c2 2 ˆ ≤ Proof. Since ϕˆ satisfies |ϕ(x)| ˆ ≤ 1+|x| α , ∀ x ∈ R, it is easy to show that k∈Z |ϕ(x+k)| M, for some positive constant M . On the other hand as ϕ(x) ˆ 6= 0 on the interval [b, b + 1] which is compact and ϕˆ is continuous, hence there exists a real number m > 0 such that 2 |ϕ(x)| ˆ ≥ m, ∀ x ∈ [b, b + 1]. Thus, we get X 0 0 and B < ∞, which means that V (ϕ) is closed and {ϕ(· − k)} forms a Riesz basis for V (ϕ) (see [12] and [5] section 5.3). Let k ∈ Z and x ∈ [0, 1]. Then we have c1 ≤ |ϕ(x + k)| ≤ 1 + |x + k|α
(
c1 , for k < 0 1+|k+1|α c1 , for k ≥ 0 1+|k|α
which implies that X
ess sup{|ϕ(x + k)|/x ∈ [0, 1]} < ∞.
k∈Z
Hence ϕ ∈ F ⊆ Wc (L1 ) and using theorem 1.4, V (ϕ) ⊆ Wc (L2 ). Theorem 3.4. Suppose ϕ ∈ F also satisfies an additional condition namely X ˆe + k) 6= 0, ∀ t ∈ [0, 1] ϕ(t
(3.1)
k∈Z
where ϕ(x) e = ϕ(−x). Let Tϕ be the infinite matrix defined by (Tϕ )jk = ϕ(j − k) for j, k ∈ Z. Then Tϕ defines an invertible, Laurent operator on `2 (Z). 5
Proof. Let ϕ ∈ F. Consider the series a(t) =
X
ˆe + k). ϕ(t
k∈Z
ˆe is continuous This series converges uniformly on [0, 1] by Weierstrass M-test. Further τk ϕ P ˆe is convergent to a continuous function for each k ∈ Z. Then it follows that a = k∈Z τk ϕ on [0, 1] which also satisfies a(0) = a(1). Corresponding to this function a, we can define a Laurent operator Tϕ by Z (Tϕ )jk =
1
a(t)e−2πi(j−k)t dt = ϕ(j − k), j, k ∈ Z.
0
As ϕ is a continuous function satisfying 3.1 (b), the Poisson summation formula holds for ˆe viz ϕ X X ˆe + k) = ϕ(k)e2πikt , ∀ t ∈ [0, 1]. ϕ(t a(t) = k∈Z
k∈Z
In other words, ‘a’ is a symbol of Tϕ . Since a(t) 6= 0, ∀ t ∈ [0, 1] and a is continuous which defines Tϕ . By theorem 2.2, Tϕ is invertible. Definition 3.5. Let F1 denote the class of functions in F satisfying (3.1). x2
Example 3.6. The functions ϕ2n defined in Example 3.2 and h0 = e− 2 belong to F1 . Let xj = j, j ∈ Z. Then the following result can be proved in a straight forward manner using theorem 1.4 and 3.4. Corollary 3.7. If ϕ ∈ F1 , then Z is a stable set of sampling for V (ϕ). Notation: Let ϕ ∈ F1 and X = {xj ∈ R/j ∈ Z}. Let UX denote the operator on `2 (Z) given by the matrix αjk = ϕ(xj − k). Definition 3.8. Let δ > 0. Then a set Y = {yj ∈ R/j ∈ Z} is said to be in a δneighbourhood of X if |xj − yj | < δ, ∀ j ∈ Z. Next we prove the following theorem. c1 Theorem 3.9. Suppose ϕ is a continuous function on R such that |ϕ(x)| ≤ 1+|x| α for all x ∈ R, for some c1 > 0, α > 1. Then ∀ > 0, ∃ δ > 0 such that whenever Y is in a δ-neighbourhood of X, ||UX − UY || < .
Proof. Let ϕ be a function satisfying all the hypotheses. Then, ∃ M > 0 such that X |ϕ(x + k)| ≤ M < ∞, ∀ x ∈ R. k∈Z
6
Therefore, for given > 0, there exists N ∈ N such that X |ϕ(x + k)| < , ∀ x ∈ R. 4 |k|>N
Now, consider X
X
|ϕ(xj − k) − ϕ(yj − k)| =
|ϕ(xj − k) − ϕ(yj − k)|
|k|≤N
k∈Z
+
X
|ϕ(xj − k) − ϕ(yj − k)|
|k|>N
X
0 such that whenever |x − y| < δ, x, y ∈ R ⇒ |ϕ(x) − ϕ(y)| < r. Take any set Y such that it is a δ-neighbourhood of X which implies that X |ϕ(xj − k) − ϕ(yj − k)| ≤ , ∀ j ∈ N 2 |k|≤N
Now, we have X
|ϕ(xj − k) − ϕ(yj − k)| < , ∀ j ∈ N.
(3.2)
k∈Z
P
Notice that k∈Z |ϕ(x+k)| is uniformly bounded in x and |xj −yj | < δ, for j ∈ Z. Instead of considering the sum over k, we shall consider the sum over j and proceed as before, we shall show that X |ϕ(xj − k) − ϕ(yj − k)| < , ∀ k ∈ N. (3.3) j∈Z
√ 2 And wePknow that || A || ≤ αβ, P where A is a bounded operator on ` (Z), α = maxj∈Z k∈Z |Ajk | and β = maxk∈Z j∈Z |Ajk |. Thus, using (3.2) and (3.3), we get ||UX − UY || < .
Corollary 3.10. Let ϕ ∈ F1 and X be a stable set of sampling for V (ϕ) such that UX is invertible. Then ∃ δ > 0 such that whenever Y is in a δ-neighbourhood of X, UY is also invertible and Y is a stable set of sampling for V (ϕ).
7
Proof. Let Take = ||U 1−1 || . Then by theorem 3.9 and 2.3, UY is invertible. Thus the X required result follows from theorem 1.4. Corollary 3.11. Let ϕ ∈ F1 . Then ∃ δ > 0 such that if X is in a δ-neighbourhood of Z, then X is a stable set of sampling for V (ϕ) and UX is invertible. Proof. By theorem 3.4, T = UZ is invertible. Now applying Corollary 3.10, we get the required result. Corollary 3.12. Suppose ϕ ∈ F1 and UX ∈ Π(Pn ). Then ∃ δ > 0 such that whenever Y is in a δ-neighbourhood of X, UY ∈ Π(Pn ). Proof. Let M = supn≥N (||(Pn UX Pn )−1 || + ||Pn ||). Take < min{ ||U 1−1 || , M12 }. Then the X result follows from the theorem 3.9 and 2.4. Theorem 3.13. Let ϕ ∈ F1 . Define a(t) = ω(e2πit ) =
X
ϕ(n)e2πint .
n∈Z
Suppose ϕ satisfies one of the following conditions: (i)
P
n6=0
|ϕ(n)| ≤ |ϕ(0)|
or (ii) ω(e2πit ) is a (complex) constant multiple of a real valued function. Then the projection method converges for the operator T where Tjk = ϕ(j − k), j, k ∈ Z, that is T = UZ . Proof. Since ϕ ∈ F1 , it follows from theorem 3.4, that T is an invertible Laurent operator. P 2πit ˆ Also ω(e ) = a(t) = k∈Z ϕ(t ˜ + k) is a continuous symbol of T. Further if (i) holds, then X |ω(e2πit ) − ϕ(0)| = | ϕ(n)e2πint | n6=0
≤
X
|ϕ(n)|
n6=0
≤ |ϕ(0)|. Thus the curve ω lies inside the disc with centre ϕ(0) and radius of |ϕ(0)|. Notice that the circle passes through the origin and ω(e2πit ) 6= 0, ∀ t ∈ [0, 1], which means the curve 8
ω can not pass through the origin. Thus the winding number of ω relative to zero equals to zero. Hence the projection method for T converges. If we assume (ii), then ω(e2πit ) = kg(t) 6= 0, for all t ∈ [0, 1] for some k ∈ C. It is easy to see geometrically that ω is contained in as line segment not passing through the origin and hence the winding number of ω relative to zero is equal to zero. Hence the projection method for T converges in this case also. Examples 3.14. • The Gaussian function g(x) =
2
−αx e√ πα
, α > 0 satisfies both the conditions. 2
• We can choose ϕ to be the inverse Fourier transform of g where g(ξ) = kpn (ξ)e−ξ where k is constant and pn (ξ) = a0 + a2 ξ 2 + ... + a2n ξ 2n , a2m ≥ 0 or even more generally we can choose pn such that pn (ξ) ≥ 0, ∀ ξ ∈ R. • ϕ = χ[− 1 , 1 ] ∗ χ[− 1 , 1 ] ∗ χ[− 1 , 1 ] 2 2
2 2
2 2
• The projection method does not converge when ϕ = χ[a,a+1] ∗ χ[a,a+1] , where a ∈ Z \ {− 12 }. In fact, it is easy to show that ϕ ∈ F, supp(ϕ) ⊆ [2a, 2a + 2] and 2 ω(e2πit ) = ϕ(2a + 1)e2πi(2a+1)t , which in turn implies that the winding number of ω relative to zero is nonzero.
4
A Reconstruction Procedure
In this section, we wish to present an iterative reconstruction algorithm for reconstructing function f ∈ V (ϕ), ϕ ∈ F1 . In order to do this, we prove a theorem which is useful in computing the solution of a finite matrix equation recursively. First, we introduce the notations. Consider an (n + 2) × (n + 2) real matrix equation given by x1 y1 x2 y2 a11 a12 · · · a1,n+2 a21 a22 · · · a2,n+2 · = · (4.1) ··· · · · · an+2,1 an+2,2 · · · an+2,n+2 xn+2 yn+2 We first write this in terms of block marices. Let Ri = (ai2 ai3 ... ai,n+1 ), CiT = (a2i a3i ... an+1,i ), i = 1, ..., n + 2. Let XnT = (x2 x3 ... xn+1 ) and YnT = (y2 y3 ... yn+1 ). Let Tn+1 be the given matrix Tn+1 = (aij )(n+2)×(n+2) . a22 · · · a2,n+1 . Then the matrix equation (4.1) becomes ··· Let Tn be the matrix an+1,2 · · · an+1,n+1 Tn+1 X = Y. Now, we prove the following theorem. 9
Theorem 4.1. Suppose Tn and Tn+1 are invertible. Then the solution X of matrix equation Tn+1 X = Y is given by X T = (x1 , XnT , xn+2 ) where x1 = xn+2 =
(y1 − R1 Tn−1 Yn )α − (yn+2 − Rn+2 Tn−1 Yn )β (a11 − R1 Tn−1 C1 )α − β(an+2,1 − Rn+2 Tn−1 C1 ) (a11 − R1 Tn−1 C1 )(yn+2 − Rn+2 Tn−1 Yn ) − (y1 − R1 Tn−1 Yn )(an+2,1 − Rn+2 Tn−1 C1 ) (a11 − R1 Tn−1 C1 )α − β(an+2,1 − Rn+2 Tn−1 C1 )
Xn = Tn−1 Yn − Tn−1 C1 x1 − Tn−1 Cn+2 xn+2 α = (an+2,n+2 − Rn+2 Tn−1 Cn+2 ) and β = (a1,n+2 − R1 Tn−1 Cn+2 ).
Proof. The equation Tn+1 X = Y can be written as a system of equations in the following manner a11 x1 + R1 Xn + a1,n+2 xn+2 = y1 C1 x1 + Tn Xn + Cn+2 xn+2 = Yn an+2,1 x1 + Rn+2 Xn + an+2,n+2 xn+2 = yn+2
(4.2) (4.3) (4.4)
Since Tn is invertible, (4.3) can be solved to obtain Xn Xn = Tn−1 Yn − Tn−1 C1 x1 − Tn−1 Cn+2 xn+2
(4.5)
Substitute (4.5) in (4.2) and (4.4), we get a(n) x1 + b(n) xn+2 = u(n) c(n) x1 + d(n) xn+2 = v (n)
(4.6) (4.7)
where a(n) = a11 − R1 Tn−1 C1 , b(n) = a1,n+2 − R1 Tn−1 Cn+2 , c(n) = an+2,1 − Rn+2 Tn−1 C1 , d(n) = an+2,n+2 − Rn+2 Tn−1 Cn+2 , u(n) = y1 − R1 Tn−1 Yn , and v (n) = yn+2 − Rn+2 Tn−1 Yn . The above set of equations has a unique solution provided a(n) d(n) − b(n) c(n) 6= 0, we Tn+1 , which in turn will give required result. show that a(n) d(n) − b(n) c(n) = det det Tn Now consider, Tn+1 . This can be written as a product of two matrices A and B, where a11 R1 Tn−1 a1,n+2 1 0 0 I Cn+2 , B = 0 Tn 0 A = C1 −1 an+2,1 Rn+2 Tn an+2,n+2 0 0 1 det A =
det Tn+1 det B
=
det Tn+1 , det Tn
it is enough to show that a(n) d(n) − b(n) c(n) = det A.
If U (n) denotes a column vector of row ‘n’ and V (n) denote a row vector of column ‘n’ and In denotes the identity matrix of order ‘n’, then we have In U (n) (n) = d − V (n) U (n) (4.8) V d (n) U I n n (n) (n) and U ). (4.9) (n) = (−1) (d − V d V 10
Let R1 Tn−1 = (b12 ... b1,n+1 ) and Rn+2 Tn−1 = (bn+2,2 ... bn+2,n+1 ). Expand determinant of A with respect to first row. The following coefficients are calculated using (4.8) and (4.9). Coefficient of a11 in det A is an+2,n+2 −Rn+2 Tn−1 Cn+2 , coefficient of a1,n+2 in det A is (−1)n+1 [(−1)n (an+2,1 − Rn+2 Tn−1 C1 )] which is equal to Rn+2 Tn−1 C1 − an+2,1 and coefficient of b1j , j = 2, ..., n + 1 in det A is given by (−1)aj1 [an+2,n+2 − bn+2,n+1 an+1,n+2 − ... − bn+2,j+1 aj+1,n+2 −bn+2,j−1 aj−1,n+2 −...−bn+2,2 a2,n+2 ]−aj,n+2 [a21 bn+2,2 +...+aj−1,1 bn+2,j−1 + aj+1,1 bn+2,j+1 + ... + an+1,1 bn+2,n+1 − an+2,1 ]. Now, a(n) d(n) −b(n) c(n) = a11 (an+2,n+2 −Rn+2 Tn−1 Cn+2 )−R1 Tn−1 C1 (an+2,n+2 −Rn+2 Tn−1 Cn+2 )+ R1 Tn−1 Cn+2 (an+2,1 − Rn+2 Tn−1 C1 ) − a1,n+2 (an+2,1 − Rn+2 Tn−1 C1 ) By expressing a(n) d(n) − b(n) c(n) interms of a11 , b12 , ..., b1,n+1 , a1,n+2 and we can show by direct calculation that a(n) d(n) − b(n) c(n) = det A. Thus the conclusion follows. Now, assume that ϕ ∈ F1 and it satisfies the hypothesis of theorem 3.13. Then, by corollary 3.11 the infinite matrix U, given by Ujk = ϕ(xj − k), is invertible. By corollary 3.12, U ∈ Π(Pn ). Thus the solutions of U x = y can be obtained from its finite dimensional truncation Un xn = yn = y|R(PnP ) . In order to provide the reconstruction algorithm, we know that f ∈ V (ϕ) and f = k∈Z ck ϕ(· − k). Thus reconstructing f is equivalent to determining the sequence c = {cn } ∈ `2 (Z). We know the values of ϕ(xj − k) and we need to reconstruct ‘f ’ from these samples f (xj ). Now, Un xn = yn can be written as a−n,−n R−n a−n,n x−n y−n C−n Un−1 Cn Xn−1 = Yn−1 an,−n Rn an,n xn yn where ajk = ϕ(xj − k). Thus we obtain the following algorithm using theorem 4.1.
5
Algorithm
Let Zn = Un−1 Yn , Dn = Un−1 C−(n+1) , En = Un−1 Cn+1 , U00 = a00 and ei denote the vector whose i th element is one and other entries are zero. We further assume that every (n) Un is invertible and define Wi = Un−1 ei , i = −n, ..., n. Let N be a large positive integer. Step-1 Input the sample set X = {xj /j = −N, ..., N } and Y = (f (xj ))N j=−N . Step-2 Define Ujk = ϕ(xj − k), j, k = −N, ..., N Z0 =
C−1 C1 e0 f (x0 ) (0) ; D0 = ; E0 = ; W0 = ; U00 U00 U00 U00 11
Step-3 For n = 1, ..., N do the following Step-3.1 To calculate U −1 ei , for i = −n, ..., n. Take a(n) = a−n,−n − R−n Dn−1 , b(n) = a−n,n − R−n En−1 , c(n) = an,−n − Rn Dn−1 , d(n) = an,n − Rn En−1 . Step-3.2 For i = −n, ..., n. Let (n−1)
xn Xn−1 (n)
Wi
(n−1)
)d(n) − (ei (n) − Rn Wi )b(n) a(n) d(n) − b(n) c(n) (n−1) (n) (n−1) (n) (ei (n) − Rn Wi )a − (ei (−n) − R−n Wi )c = a(n) d(n) − b(n) c(n) (n−1) = Wi − Dn−1 x−n − En−1 xn
x−n =
(ei (−n) − R−n Wi
T = (x−n Xn−1 xn )T
Step-3.3 Define Zn = Dn = En =
n X i=−n n X i=−n n X
(n)
yi Wi
(n)
C−(n+1) (i)Wi (n)
Cn+1 (i)Wi
i=−n
Step-4 Let dN k = ZN (k), k = −n, ..., n. Reconstruct frN (x)
=
N X
dN k ϕ(x − k)
k=−N
and plot frN . P Let f (x) := k∈Z ck ϕ(x−k), x ∈ R be the function to be reconstructed and frN (x) := PN N k=−N dk ϕ(x − k), x ∈ R be the function reconstructed using the above algorithm with 2N + 1 sample values {f (xj ) : |j| ≤ N }. dN |k| ≤ N k Let c = (ck )k∈Z , dN = 0 otherwise N 2 Then c, d ∈ ` (Z) and X f (x) − frN (x) = (ck − dN k )ϕ(x − k), x ∈ R. k∈Z
Thus f − frN ∈ V (ϕ). The error in reconstruction can be defined as EN = ||f − frN ||2 and −frN ||2 . relative error by ||f ||f ||2 12
In practice, it may be difficult to estimate EN directly. Hence we make use of the following observations: 1. Since ϕ satisfies the conditions of Definition 1.3 and f − frN ∈ V (ϕ), we have A0 ||c − dN ||`2 (Z) ≤ ||f − frN ||2 ≤ B 0 ||c − dN ||`2 (Z) . 2. Since X is a stable set of sampling for V (ϕ) (see Definition 1.2), X A||f − frN ||22 ≤ |f (xj ) − frN (xj )|2 ≤ B||f − frN ||22 . xj ∈X
Hence in place of ||f − frN ||2 , we can use any one of the following estimates: ! 21 0 EN = ||c − dN ||`2 (Z) :=
X
2 |ck − dN k |
k∈Z
12
00 EN =
X
|f (xj ) − frN (xj )|2
xj ∈X
Again, in practice, we have to replace the infinite series in the above expressions by suitable finite (partial) sums. The above algorithm is implemented using Matlab and results are shown in figure 1, 2 and figure 3. We plot the original function f and the corresponding reconstructed function Pk=15 N −x2 fr in the following figures. In figure 1, we take ϕ(x) = e , f (x) = k=−15 ck ϕ(x − k), where c0 = 1, c3 = 2, c5 = − 23 , c9 = 1, c11 = −1, c13 = −1, c14 = − 12 , c−5 = − 21 , c−8 = 1, c−9 = 1, c−11 = 1, c−13 = 12 , c−15 = −1 and remaining ci = 0 with δ = 0.45. We also 0 00 and EN for various values of N in table 1. tabulate EN P 2 In figure 2, we take ϕ(x) = (4x2 − 2)e−x ,f (x) = k=15 k=−15 ck ϕ(x − k), where c0 = 1, c3 = 23 , c5 = − 21 , c9 = 32 , c11 = −1, c13 = −1, c14 = 1 12 , c−5 = − 12 , c−8 = 1, c−9 = 1, 0 c−11 = 1, c−13 = 12 , c−15 = −2 and remaining ci = 0 with δ = 0.1. We also tabulate EN 00 and EN for various values of N in table 2. P In figure 3, we take ϕ = χ[− 1 , 1 ] ∗ χ[− 1 , 1 ] ∗ χ[− 1 , 1 ] , f (x) = k=15 k=−15 ck ϕ(x − k), where 2 2 2 2 2 2 1 1 3 c0 = 1, c1 = 2 , c2 = − 2 , c3 = 2, c5 = − 2 , c6 = 1, c8 = 3, c9 = 1, c11 = −1, c13 = −1 c14 = − 12 , c−1 = −1, c−3 = 3, c−5 = − 12 , c−7 = −1, c−8 = 1, c−9 = 1, c−11 = 1, c−13 = 21 , 0 00 c−15 = −1 and remaining ci = 0 with δ = 0.47. We also tabulate EN and EN for various values of N in table 3.
13
In the following figures (a) denotes the original function and (b) denotes the reconstructed function.
14
00 . In all the above illustrations, we have taken |j| ≤ 1000 to calculate EN
N 5 10 15 20 30 50 100 300 600
0 EN 2.7386 2.3251 2.2365 ×10−15 2.9205 ×10−15 2.8164 ×10−15 2.9645 ×10−15 2.9645 ×10−15 2.9645 ×10−15 2.9645 ×10−15
00 EN 3.6458 2.3823 2.5428 ×10−15 3.7155 ×10−15 2.5112 ×10−15 2.1153 ×10−15 2.1331 ×10−15 2.1331 ×10−15 2.1331 ×10−15
N 5 10 15 20 30 50 100 300 600
Table 1
0 EN 4.3045 2.2617 2.4863 ×10−15 2.4918 ×10−15 2.4216 ×10−15 2.4698 ×10−15 2.4698 ×10−15 2.4698 ×10−15 2.4698 ×10−15
Table 3
15
00 EN 3.7464 1.6811 1.1985 ×10−15 1.4652 ×10−15 1.4530 ×10−15 1.4180 ×10−15 1.4180 ×10−15 1.4180 ×10−15 1.4180 ×10−15
N 5 10 15 20 30 50 100 300 500 700 900
0 EN 3.4278 3.3352 1.1132 ×10−14 1.9185 ×10−14 3.5195 ×10−14 6.6224 ×10−14 9.9978 ×10−14 1.4992 ×10−13 1.4919 ×10−13 1.4919 ×10−13 1.4919 ×10−13
00 EN 7.3028 5.6506 2.2548 ×10−14 4.2904 ×10−14 8.0473 ×10−14 1.4812 ×10−13 2.2335 ×10−13 3.0936 ×10−13 3.1109 ×10−13 3.1109 ×10−13 3.1109 ×10−13
Table 2 Remark 5.1. In this paper , we have shown that if U is invertible and if the projection method converges for U , then it is possible to obtain a stable reconstruction of a function f ∈ V (ϕ) from its samples. This again leads naturally to a question: “Which operators U have these properties?”. Sufficient conditions for this have been studied when U is a Toeplitz, block Toeplitz, Laurent, block Laurent operators in [4], [9] and tridiagonal operators in [3]. In future, if the above problem is solved for other class of operators, one can obtain reconstruction procedure for large class of functions. Acknowledgement The authors thank National Board for Higher Mathematics, Department of Atomic Energy No: 48/4/2004/R & D-II/2119 for the financial grant.
References [1] A. Aldroubi and K. Gr¨ochenig, Beurling-Landau-Type Theorems for Non-Uniform Sampling in Shift Invariant Spline Spaces, J. Fourier Anal. and Appl., 6 (1):93 - 103 (2000). [2] A. Aldroubi and K. Gr¨ochenig, Non-Uniform Sampling and Reconstruction in ShiftInvariant Spaces, SIAM Rev., 43 (4):585 - 620 (2001). [3] R. Balasubramanian, S.H. Kulkarni and R. Radha, Solution of a tridiagonal operator equation, Linear Algebra and its Applications, 414(1):389 - 405 (2006). [4] A.Bottcher and B.Silberman, Introduction to large truncated Toeplitz matrices, Springer, New York (1998). [5] I. Daubechies, Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics(SIAM), Philadelphia, PA (1992). 16
[6] H.G. Feichtinger, New results on regular and irregular sampling based on Wiener amalgams, Proc. Conf. Function Spaces (K. Jarosz, Ed.), Lecture Notes in Mathematics, No. 136:107 - 121, Springer-Verlag, New York (1991). [7] H.G. Feichtinger, Wiener amalgams over Euclidean spaces and some of their applications, Proc. Conf. Function Spaces (K. Jarosz, Ed.), Lecture Notes in Mathematics, No. 136:123 - 137, Springer-Verlag, New York (1991). [8] K. Gr¨ochenig and H. Schwab, Fast local reconstruction methods for nonuniform sampling in shift-invariant spaces, SIAM J. Matrix Anal. and Appl., 24 (4):899 - 913 (2003). [9] I. Gohberg, S. Goldberg and M.A. Kaashoek, Basic classes of linear operators, Brikhauser Verlang (2003). [10] M.I. Kadec, The exact value of Paley-Wiener constants, Soviet Math. Dokl., 5:559 561 (1964). [11] H. Landau, Necessary density conditions for sampling and interpolation of certain entire functions, Acta Math., 117 : 37 - 52 (1967). [12] S.G. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2 (R), Trans. Amer. Math. Soc., 315: 69 - 97 (1989). [13] C.E. Shannon, Communication in the presence of noice, Proc. IRE, 37(1):10-21 (1949). [14] L. Schumaker, Spline Functions:Basic Theory, Wiley-Interscience, Boston (1981). [15] E.M. Stein and G. Weiss, Introduction to Fourier Analysis on Euclidean spaces, Princeton University Press, Princeton, NJ (1971).
17
