Symbol approach in a signal-restoration problem involving block Toeplitz matrices∗ Vincenza Del Prete †, Fabio Di Benedetto ‡, Marco Donatelli §, Stefano Serra-Capizzano
¶
August 3, 2014
Abstract We consider a special type of signal restoration problem where some of the sampling data are not available. The formulation related to samples of the function and its derivative leads to a possibly large linear system associated to a nonsymmetric block Toeplitz matrix which can be equipped with a 2 × 2 matrix-valued symbol. The aim of the paper is to study the eigenvalues of the matrix. We first identify in detail the symbol and its analytical features. Then, by using recent results on the eigenvalue distribution of block Toeplitz matrix-sequences, we formally describe the cluster sets and the asymptotic spectral distribution of the matrix-sequences related to our problem. The localization areas, the extremal behavior, and the conditioning are only observed numerically, but their behavior is strongly related to the analytical properties of the symbol, even though a rigorous proof is still missing in the block case.
1
Introduction and description of the problem
We consider a special type of signal restoration problem, where a finite number of samples of a band limited signal is lost. A band limited signal is a function belonging to the space Bω of the square-integrable functions, whose Fourier transform is supported in [−ω, ω]. Functions belonging to this space can be represented by the Shannon series, which is an expansion in terms of the orthonormal basis of translates of the sinc function. The coefficients of the expansion, called sampling formula, are the samples of the function at a uniform grid on R, with “density ω/π”(Nyquist density). This theory has been extended by replacing the orthonormal basis with more general families, like Riesz bases or frames, formed by translates of a single or more functions and other sampling formulae have been obtained (see the one-channel (4) and the two-channel (15) formulas and reference [1]). Frames, unlike Riesz or orthonormal bases, are overcomplete and their redundancy allows the recovery of missing samples. The problem of the recovery of missing samples has been first investigated by Ferreira in [8] where the author shows that, under suitable oversampling assumptions, any finite set of samples can be recovered from the others. The recovery technique uses an oversampling one-channel formula (see (4)) and consists in solving a structured linear system (I − S)X = B where S is a positive definite matrix, whose eigenvalues lie in the interval (0, 1). If the missing samples are equidistant, this matrix is a scalar symmetric Toeplitz matrix associated to a scalar symbol (a characteristic function); thus, by standard localization and by distribution results [4], the eigenvalues lie in the convex hull of the symbol range and are clustered at the range. ∗ This
is a preprint of a paper published in J. Comput. Appl. Math., 272 (2014), pp. 399416. di Matematica, Via Dodecaneso 35, 16146 Genova (ITALY) (E-mail:
[email protected]) ‡ Dipartimento di Matematica, Via Dodecaneso 35, 16146 Genova (ITALY) (E-mail:
[email protected]) § Dipartimento di Scienza e alta Tecnologia, Via Valleggio 11, 22100 Como (ITALY)
[email protected]) ¶ Dipartimento di Scienza e alta Tecnologia, Via Valleggio 11, 22100 Como (ITALY)
[email protected]) † Dipartimento
1
(E-mail: (E-mail:
In a further step, Santos and Ferreira considered the case of the two-channel derivative oversampling formula [10], in which a band limited function is expanded in terms of the translates of two functions (generators) (see (15)) and the coefficients are the samples of the function and its derivative. The authors of [10] show that a finite number of missing samples either of the function or of the derivative can be recovered, solving in each case a nonsingular linear system. Their technique extends in a natural way to the case where the missing samples come both from the function and its derivative, leading to a system (I − S)X = B, where the matrix S (see (19)) is a 2 × 2 block matrix depending on the generators and on the position of the missing samples. The unknowns X are the missing samples of the function and of its derivative. In [1] Brianzi and Del Prete, using Ferreira’s technique, studied the stability of the matrix in the twochannel derivative case, especially when the missing samples are located in positions U = {m i1 , m i2 , . . . , m in }, where m is an integer (interleaving factor). The practical interest for studying these cases lies in the technique of interleaving the samples of a signal, prior to their transmission or archival; the advantage of this procedure is that the transmitted (or stored) information becomes less sensitive to the burst errors that typically affect contiguous set of samples [8]. The authors have obtained estimates of the minimum and maximum eigenvalues of the block submatrices of the matrix S, showing that in some cases the matrix reduces to a (block) lower triangular matrix. Moreover they performed several numerical experiments on the dependence of the eigenvalues on the parameters of the problem and on the reconstruction of the signal, also in the ill-conditioned case of contiguous missing samples. In this paper we consider the case of missing samples of the function and of its derivative at equidistant points. The formulation of the related signal restoration problem leads to a nonsymmetric block Toeplitz matrix, possibly large. Here we perform three steps. First we construct an explicit expression of the symbol f , by interpreting the entries of the matrix as Fourier coefficients of the symbol. As a second step, we study the related symbol and, as a third step, we use the information for giving a spectral characterization of the matrices S in terms of the different parameters. More specifically, concerning the second step, we prove that the double channel matrix is similar (via a permutation transform) to a standard block Toeplitz matrix Tn (f ) with 2-by-2 nonsymmetric matrix-valued symbol f (θ) associated with two main parameters: the oversampling parameter r ∈ (0, 1) (see (13)) and the interleaving factor m (positive integer), both depending on the problem. For any choice of the parameters the eigenvalues λ1 (f (θ)) and λ2 (f (θ)) of f (θ) have the following features: (i) their range is real (proved formally), (ii) their range consists of two intervals contained in (0, 1) (proved formally only for special sets of the parameters). For specific choices of the parameters it can also be seen that the range of λ1 (f (θ)) is made by two single points V1 , V2 and the range of λ2 (f (θ)) is made by two single points V3 , V4 with 0 < V1 ≤ V2 ≤ V3 ≤ V4 < 1 and where the Lebesgue measure of the θ-values such that λ(f (θ)) = Vj is a positive number Mj for every j = 1, 2, 3, 4. In addition this behavior of the symbol (though not proved rigorously for every possible choice of the parameters) seems to be general. On the other hand, the numerical tests show that 1. all the eigenvalues of Tn (f ) are contained in (V1 , V4 ); 2. the global spectrum of the sequence {Tn (f )} is distributed as the symbol f (see Theorem 2.3 and Theorem 5.1) and it is clustered to the set {V1 , V2 , V3 , V4 } and the number of the eigenvalues which are close within to Vj is equal to 2Mj n+o(n) (the quantity named o(n) seems to grow only logarithmically with n) where the Mj ’s are the Lebesgue measures mentioned above with M1 + M2 + M3 + M4 = 2π i.e. the total measure of the definition set of the symbol; 3. the minimal eigenvalue of Tn (f ) converges to V1 monotonically and exponentially with respect to n; 2
4. the maximal eigenvalue of Tn (f ) converges to V4 monotonically and exponentially with respect to n; 5. as the parameter m that characterizes the symbol f goes to infinity, the values V1 and V2 tend to r2 and the other two V3 and V4 tend to 2r − r2 . Now a brief discussion on these items is required since a lot of information can be extracted from them. As an example, if one puts together the second and the fifth items we conclude that, for m large enough, the symbol behaves as a step function, exactly as in the case of a one-channel problem. Whereas the statements contained in the second and the fifth items are proven rigorously (the second in Theorem 5.1 by using Szeg¨ olike results by Donatelli, Neytcheva, and Serra Capizzano [6], the fifth in the derivations before Subsection 5.1), the first, the third, and the fourth are only observed in practical simulations. However the above three statements carry a strong information that could be used theoretically: in fact the observed behavior in items 1, 3, 4 is typical of a Hermitian-valued symbol [16] so that one of the following alternatives has to be true: s1 (first alternative conjecture) We know that the symbol f (θ) is nonsymmetric, but it is symmetrizable i.e. there exists a matrix Q(θ) such that Q(θ)f (θ)Q(θ)−1 is symmetric. We would like to prove that the matrix Q(θ) could be chosen as a constant matrix independent of θ: however this seems unlikely. s2 (second alternative conjecture) There exist cases in which the statement s1 is not true, but the properties 1, 3, 4 are verified. In other words, the assumption of a constant transform can be weakened to prove properties 1, 3, 4. The statement s2 would be very interesting since it opens, from a theoretical viewpoint, a new research line and would pose a real challenge to the researchers working in the spectral theory of block Toeplitz sequences. The paper is organized in six sections. After this introduction, in Section 2 we recall the tools from the spectral Toeplitz machinery that are used in our problem, by giving also light on specific features of special interest in our context. In Section 3 and in Section 4 we prove that the involved matrices are similar via permutation matrices to (scalar or block) Toeplitz matrices with (scalar or matrix-valued) symbols and we derive an analytic expression of the symbol: in both cases such information is used to derive spectral results. In this direction of special interest is Subsection 4.3 where we report an analysis of the symbol, related to the projection phase in multigrid methods, which allows to derive interesting properties of its eigenvalues. The theoretical results are numerically validated in Section 5 and conclusions and open problems are discussed in the final Section 6.
2
Toeplitz sequences and symbols
This section contains the main theoretical results available in literature in their most general formulation (concerning multilevel matrices and multivariate L1 symbols). We point out that almost all the paper considers a simpler case, where the matrices have a scalar or block one-level structure and the related symbol is a piecewise continuous scalar or matrix-valued univariate function. The unexpert reader could skip this part and consider it back only for a deeper understanding of subsequent sections. Let Ms be the linear space of the complex s×s matrices and let f be a Ms -valued function of k variables, integrable on the k-cube Ik := (−π, π)k . We shall denote by L1 (s) the space L1 (I k , (2π)−k dx, Ms ) of the Ms -valued functions of k variables, absolutely integrable on Ik . The Fourier coefficients of f , given by Z b f (j) := − f (θ)e−i hj,θi dθ, j ∈ Zk (1) Ik
R R Pk where −Ik denotes the multiple integral (2π)−k Ik , i is the imaginary unit, and hj, θi = t=1 jt θt , are the entries of the k-level block Toeplitz matrices generated by f . More precisely, if n = (n1 , . . . , nk ) is a k-index 3
with positive entries, then Tn (f ) denotes the matrix of order sb n (throughout, we let n b := i X X h Tn (f ) = ··· Jn(j11 ) ⊗ · · · ⊗ Jn(jkk ) ⊗ fb(j1 , . . . , jk ). |j1 | C} = 0} . Finally we introduce the notion of Essential Numerical Range (EN R(f )) of f which is the set of all the complex numbers y so that for any positive we have µ{x ∈ Ik : ∃v so that kvk2 = 1 and v ∗ f (x)v ∈ B(y, )} > 0 where B(y, ) = {z ∈ C : |z − y| ≤ }. If s = 1 (scalar case) then σmin (f ) and σmax (f ) coincide with the essential infimum and the essential supremum of |f |. Here the word “essential” means up to zero measure sets: in general, for the sake of simplicity, we will indicate inf h and sup h instead of essinf h and esssup h for measurable functions h. Moreover in the case of s = 1, the quantity EN R(f ) coincides with ER(f ) i.e. with the standard essential range (for the notion of essential range of a scalar valued function see [9]). In the scalar case, the following Theorem 2.1 is essentially already in [2]. Theorem 2.1 [15] Suppose f ∈ L1 (s). Suppose also that there exists a straight lineSz in S the complex plane such that, denoting by H and H the two open half-planes such that C = H z H2 , it holds 1 2 1 T S EN R(f ) H1 = ∅ and 0 ∈ H1 z. Let d(z) ≥ 0 denote the distance of z from the origin (generalization of the notion of weakly sectoriality [15]). Then if σ is a singular value of Tn (f ), it holds σ ≥ d. In addition if ω is the rotation number so that z = {w ∈ C : Re w = d} and if the minimal eigenvalue of ωf (θ) + ω ¯ f ∗ (θ) is not identically equal to d, then the inequality is strict that is σ > d. Finally, we have supz∈S {d(z)} = maxz∈S {d(z)} = δ where S is the set of all separating lines and δ is defined as the distance of the complex zero from the convex hull of EN R(f ) (when s = 1 then the constant δ is the distance of the complex zero from Coh[ER(f )] where Coh[X] denotes the convex hull of a subset X of the complex plane). 4
Theorem 2.2 [15] Suppose f, g ∈ L1 (s). Consider R(f, g) = {λ ∈ C : P(f − λg)} where the proposition P(h) is true if the function h, h ∈ L1 (s), satisfies the assumptions of Theorem 2.1. Suppose also that P(g) holds. Then, for any n, the eigenvalues of Tn−1 (g)Tn (f ) belong to [R(f, g)]c i.e. to the complementary set of R(f, g). The assumptions of Theorem 2.1 are trivially satisfied when EN R(f ) ⊂ R (it suffices to choose e.g. z = {w ∈ C : Im w = 0} and H1 = {w ∈ C : Im w > 0}), but this will not happen in the two-channel case. However, consider the case g = I (the s × s identity matrix): except for pathological situations, EN R(f − λI) = EN R(f ) − λ and therefore Theorem 2.2 says that λ is in a localization region for the spectrum whenever a translation of EN R(f ) by λ is not well separated from the origin (in the sense of Theorem 2.1). In several cases, this means that the spectrum is localized in terms of the convex hull of EN R(f ). It is quite immediate to observe that the tools given in Theorem 2.1 and in Theorem 2.2 are difficult to use in pratice. The following observation gives a series of indications on how to use the results (especially in our setting). Remark 2.1 In Theorem 2.2, if f˜ is similar to f via a constant transformation and g˜ is similar to g via the same constant transformation, then Tn−1 (g)Tn (f ) is similar to Tn−1 (˜ g )Tn (f˜) and therefore for any n, the c −1 ˜ eigenvalues of Tn (g)Tn (f ) belong to [R(f , g˜)] as well. As a consequence if F denotes the set of all the pairs (f˜, g˜) satisfying the previous assumptions, then for any n, the eigenvalues of Tn−1 (g)Tn (f ) belong to T ˜ ˜)]c i.e. to the intersection of all complementary sets of R(f˜, g˜). (f˜,˜ g )∈F [R(f , g We will be interested in the particular case where g = I and f is similar to a Hermitian-valued symbol via a constant transform; then all the eigenvalues of Tn (f ) are real and belong to the interval [λmin (f ), λmax (f )]. Furthermore, if min1≤j≤s λj (f (θ)) is not constantly equal to λmin (f ), then the localization interval becomes open on the left (λmin (f ), λmax (f )]. Analogously, if max1≤j≤s λj (f (θ)) is not constantly equal to λmax (f ), then the localization interval becomes open on the right [λmin (f ), λmax (f )), so that if both min1≤j≤s λj (f (θ)) and max1≤j≤s λj (f (θ)) are not constant the localization result is refined by replacing the closed interval with the open interval (λmin (f ), λmax (f )). As we will see this is the case occurring for the symbols appearing in the one-channel and two-channel problems, even if there is theoretical gap since we have been not able to prove the existence of a constant similarity transform which makes the symbol Hermitian (see items s1, s2 and the related discussion in the Introduction). As a final remark, if µ{θ ∈ Ik : min1≤j≤s λj (f (θ)) = λmin (f )} > 0 then the monotone convergence (from above) of the minimal eigenvalue of Tn (f ) to λmin (f ) is exponentially fast; in perfect analogy, if µ{θ ∈ Ik : max1≤j≤s λj (f (θ)) = λmax (f )} > 0 then the monotone convergence (from below) of the maximal eigenvalue of Tn (f ) to λmax (f ) is exponentially fast. It is worth stressing that the latter situation occurs, both for the minimal and the maximal eigenvalues of Tn (f ), for our one-channel and two-channel problems.
2.2
Asymptotic results for the spectrum of non-Hermitian block Toeplitz sequences
Among the several asymptotic results for block Toeplitz sequences and their generalizations, we recall a theorem which is especially suited in our context since it gives the Szeg¨o formula and the cluster sets for the matrices coming from the one-channel and two-channel problems. 5
Theorem 2.3 S [6] Let f ∈ L∞ (s) with eigenvalues λj (f ), j = 1, . . . , s, and let Λn be the spectrum of Tn (f ). s Define R(f ) = j=1 R(λj (f )) as the union of the essential ranges of the eigenvalues of f . If R(f ) has empty interior part and does not disconnect the complex plane, then for every function F ∈ C0 (C) continuous with bounded support, the following asymptotic formula Z 1 X 1 F (λ) = − lim tr(F (f (θ))) dθ (3) n→∞ n bs Ik s λ∈Λn
holds. The latter relation is indicated in short using the notation {Tn (f )} ∼λ (f, Ik ). We emphasize that in our problems k = 1 and s = 1 (one-channel problem) or s = 2 (two-channel problem). In both the cases the set R(f ) is made by only qs real numbers in the interval (0, 1), at least for a wide choice of parameters, where q1 = 2 and q2 = 4. Therefore the assumptions of the above theorem are trivially fulfilled and the distribution result holds so that {Tn (f )} ∼λ (f, I1 ). We conclude this preliminary section of tools by giving a visual interpretation of the above distributional Theorem 2.3: we focus our attention to the case where k = 1 and f (θ) is piecewise continuous as in the application at hand. In fact, a practical reading of formula (3), for every test function F continuous (s) (1) with bounded support, is that the eigenvalues of Tn (f ) can be partitioned in s sets Λn , . . . , Λn , each of cardinality n (since k = 1). A suitable ordering of the pairs 2πj (l) (l) −π + , λj , j = 1, . . . , n, λj ∈ Λ(l) n , n reconstructs, approximately for n large enough, the function θ → λl (f (θ)) for every l = 1, . . . , s. In other words, the pairs 2πj (l) 2πj 2πj −π + , λj ≈ −π + , λl f −π + n n n are approximate sampling values of the function λl (f (θ)).
3
One-channel case: the symbol and the spectral analysis
First, we briefly summarize the general one-channel case and recall some stability results obtained in [8]. In the following we shall denote by χ[−T,T ] the characteristic function of the interval [−T, T ]. The Shannon one-channel oversampling formula for a signal F in Bω is F(x) =
ωt0 X sin(ω(x − nt0 )) F(nt0 ) , π ω(x − nt0 )
(4)
n∈Z
where t0 is less then the Nyquist density π/ω. We shall denote by r the oversampling parameter r=
ωt0 . π
Note that if r = 1, the family of the translates of the function sinc(x) = sin(πx)/πx becomes an orthonormal basis. Let U = {l1 , l2 , . . . , ln } be a set of integers and {F(lk t0 ) : k = 1, . . . , n} 6
(5)
be the set of missing samples. By evaluating (4) in lk t0 , k = 1, . . . , n and separating the known samples from the unknown ones, we obtain the following system (I − S)X = B,
(6)
where I is the identity n × n matrix, S(j, k) = r sinc(r(lj − lk )),
j, k = 1, . . . , n,
X = (F(lk t0 ))nk=1 , B = (b1 , b2 , . . . , bn )T , and bj = r
X
F(kt0 )sinc(r(lj − lk )),
j = 1, . . . , n.
k∈ /U
The matrix S is symmetric and positive definite. In [8], the author observes that when r is close to 1 the recovery is impossible. This is in agreement with the theory, since for r = 1 the family of translates of the sinc function is not a frame but an orthonormal basis and S is the identity matrix. In addition, the author shows that all its eigenvalues are in (0, 1) for every choice of the set U , by observing that S is always a principal submatrix of a scalar Toeplitz prolate matrix of sufficiently large order. Recall that the n × n prolate matrix is defined as M := (rsinc[r(j − k)])nj,k=1 . (7) An interesting situation occurs when U consists of equidistant points: let m be the interleaving factor and set U = {m, 2m, . . . , nm}. In this case the matrix S S := (rsinc[r(j − k)m])nj,k=1
(8)
keeps the Toeplitz structure and it can be seen as a generalization of M . For instance, in [8] the quadratic forms associated to the matrices M and S are expressed as follows Z Z 1 r 1 r ˜ 2 ∗ 2 ∗ v Mv = |φ(t)| dt, v Sv = |φ(t)| dt, (9) 2 −r 2 −r Pn−1 ˜ = Pn−1 vk eiπmkt . Note that φ(t) ˜ = φ(mt). where v = (v1 , v2 , . . . vn )T , φ(t) = k=0 vk eiπkt and φ(t) k=0 The prolate matrix M is studied in detail in [20] and is generated, in the sense of the previous section, by a scalar symbol represented by the characteristic function χ[−rπ,rπ] . This fact is related to the identity 2w + 2
+∞ X sin 2kπw k=1
kπ
cos kθ = χ[−2πw,2πw] (θ),
for |θ| = 6 2πw and w ∈ (0, 1/2)
(10)
and to the integral expression in (9) for v ∗ M v. Given the strong analogy with the formula associated to S in (9), we are tempted to conclude that its symbol too could be the characteristic function of some interval. Unfortunately things are not so straightforward. Indeed, if we try to compute directly the scalar symbol X f (θ) = r sinc(kmr)eikθ (11) k∈Z
of S by making a plain use of (10), we should assume w = rm/2 which does not belong to (0, 1/2). In reality we should proceed differently by defining L := brmc,
L0 := rm − L (fractionary part of rm, belonging to [0, 1)), 7
w := L0 /2.
With such a setting, the quantity w satisfies the required bounds. Denote by g(θ) the 2π-periodic function, that in [−π, π] is equal to χ[−2πw,2πw] (θ). Then the symbol is f (θ) =
g(θ + Lπ) + L . m
In order to obtain the explicit expression of f , we have to discuss the parity of L. When L is even, we have f (θ) = feven (θ) = {χ[−πL0 ,πL0 ] + L}/m. For odd L, we obtain f (θ) = fodd (θ) = {L + 1 − χ[−π(1−L0 ),π(1−L0 )] }/m so that fodd (θ + π) = feven (θ).
(12)
In both cases, since f is nonnegative and not identically zero, the matrix S is symmetric positive definite. (n) (n) (n) According to Remark 2.1, all the eigenvalues λ1 ≤ λ2 ≤ · · · ≤ λn of the Hermitian matrix S = Tn (f ) belong to the open interval (L/m, (L + 1)/m) = (minf, maxf ) where the interval is open since minf < maxf (all this derives from a classical integral representation of the Rayleigh quotients, for which we refer to [4]). Some spectrum localization results for S have been proved in [8] in a direct way, without the use of a symbol approach; but in the case of equidistant points, this new tool allows us to say much more. For instance, by using the classical distribution results of Szeg¨o, the graph of f gives a faithful representation of the eigenvalues of S which shows “jumps” as a function of the parameter r (which determines the value of the integer L). Moreover, we have exponential convergence of the extremal eigenvalues to the extremal values of f (since f attains minimum and maximum on a nontrivial interval). More precisely, by [12, Example 2 at page 123 and Section 3.1], since f − minf has a zero of infinite order (indeed f − minf vanishes on a non trivial interval), for any k, q fixed positive integers, we have lim
n→∞
(n) λk − minf nq = 0.
The same argument applied to maxf −f leads to the dual conclusion that, for any k, q fixed positive integers, one finds (n) lim maxf − λn−k nq = 0. n→∞
In other words the two points minf and maxf , which represent the range of the function, are exponential attractors for the extremal eigenvalues. Moreover all the spectrum is clustered at the range (it is a direct consequence of Theorem 2.3) and more precisely nL0 ± o(n) eigenvalues belong to ((L + 1)/m − , (L + 1)/m) and n(1 − L0 ) ± o(n) eigenvalues belong to (L/m, L/m + ), for any > 0 with the two o(n) terms depending on . Hence the two points of the range, minf and maxf , are the only sub-clusters of the spectra of the matrix sequence (see [4, p. 143]). A general conclusion is that the essential information on the asymptotic and localization of the spectra is totally independent on the parity of L: however this conclusion is not surprising and in fact it is expected and already written in the theoretical results of Section 2.1 and in relation (12), where it is evident that the range of the two functions fodd and feven is identical (one is the shift of the other by half a period). On the other hand the information on the domain where the minimum or the maximum is attained is important for the asymptotic eigenvector behavior (see the work in [21] for a global clustering study and the work in [3] for a single eigenvector asymptotic analysis): the eigenvectors related to the sub-cluster at (L + 1)/m are associated with the low frequencies for L even and with the high frequencies for L odd. The opposite behavior occurs for the eigenvectors related to the sub-cluster at L/m. 8
4
Two-channel case: the symbol and the spectral analysis
In [10] the authors obtained a two-channel derivative formula (see (15) below) for functions in Bω , where the sampling frequency in each channel is greater than the Nyquist frequency ω/2π (first order derivative oversampling formula). We shall denote by Fb the Fourier transform of an absolutely integrable function F Z 1 b F(t)e−itξ dt. F(ξ) =√ 2π Let t0 be a positive number less than 2π/ω and r be the ratio r=
ωt0 . 2π
(13)
Denote by ϕ e1 , ϕ e2 the functions whose Fourier transforms are 2
r r c ϕ e1 (ξ) = 1 − |ξ| χ[−ω,ω] (ξ) ω ω
r c ϕ e2 (ξ) = i 2 sign(ξ)χ[−ω,ω] (ξ). ω
Then any F ∈ Bω can be written according to the expansion √ X F(t) = 2π F(kt0 )ϕ e1 (t − kt0 ) − F 0 (kt0 )ϕ e2 (t − kt0 )
(14)
(15)
k∈Z
called derivative oversampling formula. We observe that 1/r = 2π/(ωt0 ) is the ratio between the sampling frequency and the Nyquist frequency and that r ∈ (0, 1). We shall be mainly interested in the case 1/2 < r < 1, since for 0 < r ≤ 1/2 it is possible to use one channel separately for the function and for its derivative [10]. Note that if r is close to 1, the frame is close to a Riesz basis, while if r is small, then the frame is very redundant [1]. Let (16) U = {l1 , l2 . . . , ln } ⊂ Z the positions of the missing samples and let X = (F(l1 t0 ), . . . , F(ln t0 ), F 0 (l1 t0 ), . . . , F 0 (ln t0 ))
T
(17)
be a vector containing the corresponding set of missing samples. To find a system with these samples as unknowns, we proceed as in the one-channel case. We compute the derivative of both sides of (15), evaluate the expansions in `k t0 and separate the unknown samples from the known ones (for explicit computations see [1]). Thus we obtain the system of 2n equations in 2n unknowns (I − S)X = B ,
(18)
where I is now the 2n × 2n identity matrix and S = S(U , r) is the real matrix S11 S12 S= , S21 S22
(19)
where S11 , S12 , S21 , S22 are the n × n blocks whose entries are √ √ S11 (k, j) = 2π ϕ f1 ((lk − lj )t0 ) S12 (k, j) = − 2π ϕ f2 ((lk − lj )t0 ) S21 (k, j) =
√
√
0
2π ϕ f1 ((lk − lj )t0 )
(20) 0
S22 (k, j) = − 2π ϕ f2 ((lk − lj )t0 )
k, j = 1, . . . , n. Finally B is a 2n-vector containing only known samples [1, formula (3.6) p.196]. 9
4.1
The entries of the matrix for equidistant points
From now on, we are only interested in the case where the missing samples in the set (16) are equally spaced, so that the block Toeplitz structure of S is preserved. Let U = {m, 2m . . . , nm} where m is a positive integer; then, by (13), lk t0 = 2πkmr/ω. In order to exhibit the entries of S in (20), we need to evaluate the explicit expressions of ϕ e1 and ϕ e2 and their derivatives. The former can be obtained from (14) via inverse Fourier transformation: ωx i ωx 1 h 2r (1 − r)sinc( ) + r2 sinc2 ( ) ϕ e1 (x) = √ π 2π 2π
(21)
ωx 1 xr2 sinc2 ( ). ϕ e2 (x) = − √ 2π 2π
(22)
Observe that ϕ e01 (0) = 0. The expressions of ϕ e01 and ϕ e02 are 2r 1 h (1 − r) cos(ωx) − sinc(ωx) + ϕ e01 (x) = √ 2π x ω ωx ωx i + r sinc( x) cos( ) − sinc( ) 2 2 2 1 ωx h ωx ωx i ϕ e02 (x) = − √ r2 sinc( ) 2 cos( ) − sinc( ) . 2 2 2 2π
(23)
(24)
We observe that ϕ e1 and ϕ e02 are even functions whereas ϕ e2 and ϕ e01 are odd, thus by (20) it follows that S11 and S22 are real symmetric, whereas S12 and S21 are distinct real anti-symmetric matrices. The four submatrices are all of Toeplitz nature with entries S11 = (ci−j ) with ck = 2r(1 − r)sinc(2rmk) + r2 sinc2 (rmk); 3 S12 = (ci−j ) with ck = 2πrωmk sinc2 (rmk); ω S21 = (ci−j ) with ck = πmk [(1 − r) cos(2πrmk) + (2r − 1)sinc(2rmk) − rsinc2 (rmk)]; 2 S22 = (ci−j ) with ck = r sinc(rmk)[2 cos(πrmk) − sinc(rmk)].
4.2
The matrix-valued symbol
We define the auxiliary parameters L := b2rmc,
L0 := 2rm − L (fractionary part of 2rm, in the interval [0, 1)), ( πL0 for L even α := π(1 − L0 ) for L odd.
Observe that in both cases [−α, α] ⊆ [−π, π]. Denote by fij (θ), i, j = 1, 2 the symbol related to the block Sij of S and by f the matrix-valued symbol f11 (θ) f12 (θ) f (θ) = . (25) f21 (θ) f22 (θ) The Fourier coefficients of f are fb(k) =
c fc 11 (k) f12 (k) c c f21 (k) f22 (k) 10
! k∈Z
(26)
where for each i, j (fc ij (k)) is the sequence (ck ) giving the entries of Sij in Subsection 4.1. Of course S is not in the form prescribed in (2), but is similar via a permutation matrix to Tn (f ) with f as in (25) and according to (2). In formulae we have ΠSΠT = Tn (f ) with Π of size 2n being the permutation matrix
eT1 0 .. .
Π= eTn 0
0 eT1 .. .
0 eTn
and with ej , j = 1, . . . , n, being the n vectors of the canonical basis of Cn . In the subsections below we find the explicit expressions of fij which depend on the parity of L and are piecewise defined on the subintervals (−π, −α), (−α, 0), (0, α) and (α, π). We also observe that fij is even or odd according with the parity of the sequence of its Fourier coefficients ck , so that f11 and f22 are even, whereas f12 and f21 are odd. Thus it is sufficient to give their expressions only for θ ∈ (0, π). 4.2.1
Even L
f11 (θ) =
c11 +
1 m
−
L 2m2
−
θ 2πm2
f12 (θ) =
c11
if θ ∈ (α, π)
where c11 := 2r − r2 − (1 − r) Lm0 −
f21 (θ) =
iωc1 4πrm3 (θ
if θ ∈ (0, α) f22 (θ) =
− π)
c22 +
if θ ∈ (0, α)
if θ ∈ (α, π)
0
L+θ/π 2m2
if θ ∈ (0, α)
c22
if θ ∈ (α, π)
if θ ∈ (α, π)
where c1 := L2 − 2mL, c2 := c1 + 2(1 − 4r)m + 2L and c22 := r2 − r Lm0 + 4.2.2
ir mω
L20 4m2 ,
iωθ 4πrm3 (c2 + θ/π)
if θ ∈ (0, α)
L20 4m2 .
Odd L
f11 (θ) =
c˜11 +
1 2πm2 (−θ
c˜11 +
1−r m
iωθ 4πrm3 (c4 + θ/π)
iωc3 4πrm3 (θ
if θ ∈ (0, α) f12 (θ) =
where c˜11 := 2r − r2 − (1 − r) Lm0 −
f21 (θ) =
+ α)
− π)
if θ ∈ (α, π)
ir mω
if θ ∈ (0, α)
if θ ∈ (α, π)
0
(1−L0 )2 4m2 ,
if θ ∈ (0, α) f22 (θ) = if θ ∈ (α, π)
c˜22 +
1 2πm2 (θ
r m
c˜22 +
where c3 := (1 + L)(1 + L − 2m), c4 := L2 − 2mL − 1 and c˜22 := r2 − r Lm0 + 11
− α)
if θ ∈ (0, α) if θ ∈ (α, π)
(1−L0 )2 4m2
.
4.3
A simplified matrix-valued symbol, multigrid methods, and some representations of the true symbol
As well known, several spectral results on block Toeplitz sequences (including those summarized in Section 2) are based on integral formulas in terms of the symbol and specific trigonometric polynomials. The same is true for the matrix S described in (19) since Z π 1 p(θ)∗ f (θ)p(θ) dθ, (27) v ∗ Sv = 2π −π where v = xy with x = (x1 , x2 , . . . , xn )T and y = (y1 , y2 , . . . , yn )T in Cn , whereas p(θ) = x(θ) with x(θ) y(θ) and y(θ) trigonometric polynomials isometrically associated to the vectors x and y, respectively, as follows: x(θ) =
n X
xk eikθ ,
y(θ) =
k=1
n X
yk eikθ .
k=1
Below we find an interesting analog of the classical Toeplitz representation (27) in terms of a new 2 × 2 matrix-valued symbol A(t) having a very simplified expression. Let m be the interleaving factor in the positions U = {m, 2m . . . , nm} of the missing samples. Theorem 4.1 Let r ∈ (0, 1) and let S be the matrix in (19). Then the following integral formula Z r v ∗ Sv = P (t)∗ A(t)P (t) dt
(28)
−r
holds, where A(t) =
1 − |t| iωt(1 − |t|)/r
−irsgn(t)/ω |t|
t ∈ [−1, 1],
(29)
yk e2πimkt .
(30)
P (t) = (Φ(t), Ψ(t))T and Φ(t) =
n X
xk e2πimkt ,
Ψ(t) =
k=1
k=1
Proof.
n X
From the expressions of the entries of S in (20) with `k = mk, k = 1, . . . n, we have n n X X √ xk xj ϕ f1 (m(k − j)t0 )+ xj yk ϕ f2 (m(k − j)t0 )+ v ∗ Sv = 2π
+
k,j=1
k,j=1
n X
n X
0
xk yj ϕ f1 (m(k − j)t0 )+
k,j=1
0
yk yj ϕ f2 (m(k − j)t0 ) .
k,j=1
By inverse transforming in (14) and changing variables in the integral, we obtain Z √ Z r √ r r ϕ f1 (x) = 2π (1 − |t|)eixωt/r dt ϕ f2 (x) = 2π i sign(t)eixωt/r dt . ω −r −r From this, by using the identity gb0 (ξ) = iξ gb(ξ) we deduce similarly for the derivatives Z √ √ Z r ω r 0 0 ixωt/r t(1 − |t|)e dt ϕ f2 (x) = − 2π |t|eixωt/r dt . ϕ f1 (x) = 2πi r −r −r 12
Thus, by reminding that t0 = 2πr/ω by (13), the quadratic form can be written as Z Z r n n X X r r sign(t) v ∗ Sv = (1 − |t|) xj yk e2πim(k−j)t dt xk xj e2πim(k−j)t dt − i ω −r −r k,j=1
k,j=1
+i
ω r
Z
r
t(1 − |t|) −r
n X
xk yj e2πim(k−j)t dt +
Z
r
|t| −r
k,j=1
n X
yk yj e2πim(k−j)t dt .
k,j=1
Hence, by (30), the expression of the quadratic form becomes Z r Z r r 2 ∗ (1 − |t|)|Φ(t)| dt − i sign(t)Φ(t)Ψ(t) dt v Sv = ω −r −r (31) +i
ω r
Z
r
Z
r
t(1 − |t|)Φ(t)Ψ(t) dt + −r
|t||Ψ(t)|2 dt .
−r
Thus formula (28) is proved. When 1/2 < r < 1 (the main case of interest) the integration domain in (31) can be further restricted. Corollary 4.1 Let 1/2 < r < 1 and let S, Φ, Ψ, P be as in Theorem 4.1, then Z 1/2 ˜ v ∗ Sv = P (s)∗ A(s)P (s) ds −1/2
˜ where P = (Φ, Ψ)T ; the matrix A(s) is equal to A(s) in (29) if |s| ≤ 1 − r and is the identity matrix if 1 − r ≤ |s| < 1/2. Proof. Since 1/2 < r < 1, we may write the first integral in (31) as the sum of the integrals on the intervals [−r, −1/2], [−r, r] and [r, 1/2]. Then we change the variable in the first and the third integral and use the fact that the function Φ has period 1, obtaining Z r Z 1/2 Z 1/2 Z −1+r (1 − |t|)|Φ(t)|2 dt = s |Φ(s)|2 ds + (1 − |s|) |Φ(s)|2 ds − s |Φ(s)|2 ds. −r
−1/2
1−r
This can be written as follows Z r Z 2 (1 − |t|)|Φ(t)| dt = −r
−1/2
1/2
s χ[1−r,1/2] (s) + (1 − |s|) − s χ[−1/2,r−1] (s) |Φ(s)|2 ds .
−1/2
By performing similar calculations for the other three integrals in (31), we obtain Z 1/2 h ∗ v Sv = a11 (s)|Φ(s)|2 + a12 (s)Φ(s)Ψ(s) + a21 (s)Φ(s)Ψ(s) + a22 (s)|Φ(s)|2 ds −1/2
where ( (1 − |s|) if |s| ≤ 1 − r a11 (s) = 1 if 1 − r ≤ |s| < 1/2 ( −i ωr sign(s) if |s| ≤ 1 − r a12 (s) = 0 if 1 − r ≤ |s| < 1/2
( a22 (s) =
|s| if |s| ≤ 1 − r 1 if 1 − r ≤ |s| < 1/2
( i ωr s(1 − |s|) if |s| ≤ 1 − r a21 (s) = 0 if 1 − r ≤ |s| < 1/2.
13
This concludes the proof. The results of Theorem 4.1 and Corollary 4.1 easily extend to the case where the missing samples are not equally spaced. Indeed it suffices to replace the functions Φ and Ψ in (30) by Φ(t) =
n X
xk e2πi`k t ,
Ψ(t) =
k=1
n X
yk e2πi`k t .
(32)
k=1
Now we can discuss the spectral properties of the simplified symbol A(t) in (29). The first key observation is that, since A12 (t)A21 (t) = |t|(1 − |t|), A(t) has rank equal to 1 independently of t so that the spectrum is given by 0 and 1 (that is the trace of A(t)). Second we observe that A(t) is similar to the symmetric matrix-valued symbol Asym (t) =
p 1 − |t| |t|(1 − |t|)
p
|t|(1 − |t|) |t|
.
Indeed if we prove that A(t) is similar via a constant similarity transformation to a Hermitian symbol, then by Remark 2.1 we would have that all the eigenvalues of S are real, belonging to the interval [0, 1]. Unfortunately the similarity transformation matrix is not constant with respect to t so that we cannot apply the intersection statement in Remark 2.1. This negative fact is structural in the sense that for every constant invertible matrix Q, a trivial computation shows that QA(t)Q−1 is not always Hermitian: however this fact does not imply that the same is true for f (θ). Moreover, we observe that, by a change of variable, the integral formula (28) may be written as 1 π
v ∗ Sv =
Z
πr
P (θ)∗ A(θ/π)P (θ) dθ,
(33)
−πr
where Φ(θ) :=
n X
xk ei2mkθ ,
Ψ(θ) :=
k=1
n X
yk ei2mkθ
(34)
k=1
for θ ∈ [−π, π] and P (θ) = (Φ(θ), Ψ(θ))T . It is intriguing to compare (33) with the classical representation formula (27) relating a block Toeplitz matrix and its symbol, which we remind below: Z π 1 v ∗ Sv = p(θ)∗ f (θ)p(θ) dθ 2π −π (with p(θ) different from P (θ)). We can guess that the product F (θ) = 2A(θ/π) · χ[−πr,πr] (θ)
(35)
plays the role of a symbol (a simplified symbol), but with a much simpler expression when compared with the true symbol f (θ). Therefore we try to determine the relations between the spectral properties of the block Toeplitz sequence under study and the analytical/spectral features of the simplified symbol. First we put in relation the simplified symbol and the true symbol: to this end classical tools used when designing multigrid methods are really useful (see [7]: equation (4.7) in Proposition 4.1 where g := 2m, p identically 1 and Section 6.1 therein for a discussion on the Toeplitz case). In fact the integral formula (33), in matrix terms, can be read as ΠSΠT = Tn (f ) = U T T2nm (F )U 14
˜ ⊗ I2 is the “cutting” matrix of size 4nm × 2n, with U ˜ij = 1 if and only if j = (i − 1)2m + 1 where U = U and zero otherwise. Therefore we have f (θ) =
2m−1 2m−1 1 X θ + 2πj 1 X θ + 2πj θ + 2πj F = A · χ[−πr,πr] . 2m j=0 2m m j=0 2mπ 2m
(36)
We observe that the characteristic function present in the expression (35) for F , evaluated in the various sub-intervals given by the change of variable θ 7→ θ+2πj 2m , becomes identically 1 or identically 0 most of the times so that, taking into account the definition domain of A(t), the symbol f (θ) in (36) becomes bL/2c 1 X θ + 2πj f (θ) = A m 2mπ j=−bL/2c θ − Lπ −A · χ[−π,−πL0 ] (θ) (37) 2mπ θ + Lπ · χ[πL0 ,π] (θ) −A 2mπ if L is even or f (θ)
bL/2c 1 X θ + 2πj = A m 2mπ j=−bL/2c θ + (L + 1)π +A · χ[−π,−π(1−L0 )] (θ) 2mπ θ − (L + 1)π · χ[π(1−L0 ),π] (θ) +A 2mπ
(38)
if L is odd. Formulae (37) and (38) can be considered as an alternative representation of f in place of that given in Subsection 4.2.
5
Numerical results (vs Theoretical results)
First of all we concentrate our attention on the distribution results (and a fortiori on the cluster sets). To this end the first task is to check the assumptions of Theorem 2.3: the latter means that the range R(f ) (the union of the ranges of λ1 (f (θ)) and λ2 (f (θ))) must not disconnect the complex plane. The main difficulty relies on the several parameters involved in the analytical expression of the symbol f (θ). Nevertheless the next proposition answers in the positive to our question, as it will be explicitly emphasized in Theorem 5.1. Proposition 5.1 Let αi , i = 1, . . . , l be real numbers belonging to [0, 1] with l ≤ 2l0 , both l, l0 being positive Pl integer numbers. Let si , i = 1, . . . , l be a collection of signs i.e. si ∈ {±1}, i = 1, . . . , l, let ν = i=1 si , and finally let γ be a nonzero complex value. For every i = 1, . . . , l define 1 − αi γ −1 si Ai = γsi (1 − αi )αi αi and A=
l 1X Ai . l0 i=1
15
(39)
If ν ≥ 0 then the eigenvalues of A, namely λ1 (A) and λ2 (A), are both real nonnegative with l l ≤ λ2 (A) ≤ 0 ≤ 2. 2l0 l
0 ≤ λ1 (A) ≤
(40)
Moreover for ν = 0 we have λ1 (A)
=
λ2 (A)
=
l−
Pl
i=1
αi
l0 Pl i=1 αi . l0
,
If ν < 0 then the eigenvalues of A can take either the same form as in (40) or the form λ1,2 (A) =
l ± iβ, 2l0
l ≥ β > 0. 2l0
(41)
Proof. It is enough to observe that every Ai can be seen as a sampling of the simplified symbol A(t) on a given point ti : therefore each Ai is a rank-one matrix with trace 1. Hence, by (39), the trace of A is given by l tr(A) = 0 , l while a simple check (again by (39)) shows that 0
det(l A) = (l −
l X
αi )
i=1
l X
l X αi − ν (1 − αi )αi
i=1
i=1
with
q 1 2 (42) λ1,2 (A) = tr(A) ± tr (A) − 4det(A) . 2 Pl Pl Pl Now assume ν ≥ 0. Since i=1 αi ∈ [0, l], it follows that (l − i=1 αi ) i=1 αi ∈ [0, l2 /4] with minimum Pl Pl Pl attained when i=1 αi = 0 or i=1 αi = l and maximum attained for i=1 αi = l/2. Since l X (1 − αi )αi ≥ 0 i=1
it directly follows that det(l0 A) ≤ l2 /4. In order to find a lower bound, we perform a direct computation so that, taking into account that ν ≤ l we have det(l0 A)
=
(l −
l X
αi )
l X
i=1
≥ (l −
l X
i=1
αi )
l X
i=1
= l =
αi − ν
i=1 kek22 kak22
(1 − αi )αi
i=1
αi − l
i=1
l X (αi )2 −
l X
l X
l X
(1 − αi )αi
i=1 !2
αi
i=1 T 2
− (e a) ≥ 0,
where e is the vector of length l of all ones, a is the vector whose i-th p component is αi , and where in the last line we have used the Cauchy-Schwartz inequality. Therefore tr2 (l0 A) − 4det(l0 A) ∈ [0, l] and the estimates in (40) directly follow. 16
With analogous computation the case of ν ∈ [−l, 0] leads to the estimate det(l0 A) ∈ [0, l2 /2] so that either (40) is true or (41) is true and the proof is complete. In the sequel, we use the notation Sn to emphasize that the matrix S has inner blocks of size n which can arbitrarily grow. Theorem 5.1 The eigenvalues of the sequence {Sn } are distributed as λ1 (f (θ)) and λ2 (f (θ)) that is {Sn } ∼λ (f, I1 ) according to the notations in Theorem 2.3. Proof. As already proved Sn is similar to Tn (f ) and, according to (36), the symbol f (θ) is such that, for any fixed θ, the matrix f (θ) is of the form (39) with l0 = m, l ≤ 2m and γ = iω/r. Therefore, by Proposition 5.1 the set R(f ) is a subset of the cross K in the complex plane given by [ l K = {z ∈ C : Im(z) = 0, Re(z) ∈ [0, 2]} z ∈ C : Re(z) = , Im(z) ∈ [−1, 1] . 2m Since K does not disconnect C, the same does R(f ) and the desired result follows from Theorem 2.3. Therefore, according to the discussion in Subsection 2.2, n eigenvalues of S = Sn behave as the function θ 7→ λ1 (f (θ)) and the remaining n eigenvalues of S = Sn behave as the function θ 7→ λ2 (f (θ)) or equivalently −π +
2πj (k) , λj n
≈
−π +
2πj 2πj , λk f −π + n n
are approximate sampling values of the function λk (f (θ)), k = 1, 2. At this point it is clear that our goal has to be the understanding as much as possible of the expression of λk (f (θ)), k = 1, 2. The information is in fact contained in the representations (37) and (38), and in those of Subsection 4.2. In connection with Proposition 5.1, a careful analysis shows that l ∈ {L, L + 1},
ν ∈ {0, ±1}.
In fact, in formula (36), the sub-domains
2j − 1 2j + 1 Ij = , 2m 2m
obtained from the affine transform θ 7→ θ+2πj 2m , θ ∈ [−π, π], are such that χ[−πr,πr] (θ) Ij is either identically zero or identically one most of the times. As a consequence, when ν = 0, the symbol is lower triangular so that the eigenvalues of f (θ) are the diagonal entries. i.e. L0 L2 + 02 , m 4m L0 L2 2 2r − r − (1 − r) − 02 , m 4m
λ1 (f (θ))
= r2 − r
λ2 (f (θ))
=
17
for L even and θ ∈ [−π, π]\[−α, α] or r (1 − L0 )2 , + m 4m2 (1 − r) (1 − L0 )2 2r − r2 − (1 − L0 ) , − m 4m2
λ1 (f (θ))
= r2 + (1 − L0 )
λ2 (f (θ))
=
for L odd and θ ∈ [−π, π]\[−α, α]. In the remaining domain, both for even and odd L, an explicit computation shows that r 1 , |λ1 (f (θ)) − r2 | ≤ +O m m3 1−r 1 2 |λ2 (f (θ)) − (2r − r )| ≤ +O m m3
(43) (44)
so that R(f ) is made by 4 points V1 , V2 , V3 , V4 , where for m tending to infinity we have V1 = V2 = r2 and V3 = V4 = 2r −r2 , i.e. two parabola symmetrically placed with respect to the first bisectrix with r2 ≤ 2r −r2 for r ∈ [0, 1] and equality only at r = 0, r = 1 (we recall in any case that the parameter r lies in the physical domain (0, 1)). Concerning the localization results we do not try to apply Theorem 2.2, because of its intrinsic difficulty. Instead we follow the discussion in Remark 2.1: indeed, we have been not able to prove the existence of a constant similarity transform which makes the symbol Hermitian and this leaves the door open for further investigations, in the sense of items s1 and s2 and of the related discussion in the Introduction. One of the desired properties (required by the algorithms used for treating the considered signalrestoration problem with missing data) is the fact that all the eigenvalues of S = Sn are nonnegative and strictly less than one. Indeed the asymptotic values r2 and 2r − r2 have the desired property and the same occurs for every n for the eigenvalues of S = Sn , when m is large. The problem is nontrivial when m is moderate, but a finer analysis, taking into account that L0 depends on r, shows that the two eigenvalues of the symbol have always essential infimum and essential supremum in the interval (0, 1).
5.1
Numerical evidences
The following set of numerical tests confirms that the spectrum of S = Sn is well described by the range of the two functions λ1 (f (θ)) and λ2 (f (θ)) according to Theorem 5.1. To this end the eigenvalues of the 2 × 2 matrix valued symbol f (θ) are sampled over 1000 equi-spaced points on the natural domain I1 = [−π, π]. The tests are performed using MatLab in double precision. Figures 1–3 show the eigenvalues of S = Sn and the values of λ1 (f (θ)) and λ2 (f (θ)), for several combinations of the parameters. The imaginary part of the computed eigenvalues has a size comparable to the machine precision and therefore we can safely claim that the eigenvalues of Sn are purely real. Furthermore, the agreement among the eigenvalues and R(f ) is impressive, giving evidence of the result described in Theorem 5.1 (see the discussion after the theorem for a ‘visual’ interpretation of the result that is perfectly visible in our tests). Now we focus our attention on the behavior of the extremal eigenvalues. Indeed it is proved that lim λk (Sn ) = r2 ,
m,n→∞
lim λ2n−k (Sn ) = 2r − r2
m,n→∞
(45)
for every fixed k independent of n (see also [1]), where we have assumed the ordering λ1 (Sn ) ≤ λ2 (Sn ) ≤ · · · ≤ λ2n (Sn ). Figure 4 shows that a similar behavior can be observed for the range of the two eigenvalues of f , for large m (they approach the limits in (45), symmetric around the line y = r, r ∈ [0, 1]). More specifically in Figure 4 (c) we see that for r = 0.8, n = 50 and large m the maximal eigenvalue of S = Sn differs from 18
−15
5
−15
x 10
8
−15
x 10
2
6
0
−5 0.4
0.5
0.6
0.7
0.8
0.9
1
x 10
1.5
4
1
2
0.5
0
0
−2
−0.5
−4
−1
−6
−1.5
−8 0.4
0.5
m=5
0.6
0.7
0.8
0.9
1
−2 0.4
0.5
0.6
m = 12
0.7
0.8
0.9
1
m = 21
Figure 1: The ◦ depict the eigenvalues of S for n = 100 and the × depict the range of the eigenvalues of f for r = 0.7.
−16
3
−16
x 10
4
−15
x 10
2
3
2
1.5
2 1
x 10
1
1 0.5
0
0 0 −1
−1
−0.5
−2 −2
−1
−3
−3 0.2
0.4
0.6
0.8
1
−4 0.5
0.6
r = 0.6
0.7
0.8
0.9
1
−1.5 0.7
0.8
r = 0.75
0.9
1
1.1
r = 0.9
Figure 2: The ◦ depict the eigenvalues of S for n = 100 and the × depict the range of the eigenvalues of f for m = 8
−16
1
3
−15
x 10
4
x 10
3
2 0.5
2 1 1
0
0 0 −1 −1
−0.5 −2 −1 0.4
0.5
0.6
0.7
n = 10
0.8
0.9
1
−3 0.4
−2 0.5
0.6
0.7
n = 100
0.8
0.9
1
−3 0.4
0.5
0.6
0.7
0.8
0.9
1
n = 300
Figure 3: For r = 0.7 and m = 8, × depict the range of the eigenvalues of f and the ◦ depict the eigenvalues of S for different values of n.
19
−12
1.2
10
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
−13
10
−14
10
0 −0.2 0
−15
10
0 20
40
60
80
100
−0.2 0
−16
20
(a)
40
60
80
100
10
0
20
40
(b)
60
80
100
(c)
Figure 4: For r = 0.8 and n = 50 versus m = 1, . . . , 100: (a) ◦ are λmax (S) and × are λmin (S), (b) ◦ are λmax (f ) and × are λmin (f ) (introduced in Definition 2.1), (c) the dashed curve depicts |λmax (S) − λmax (f )|/|λmax (S)| while the solid curve depicts |λmin (S) − λmin (f )|/|λmin (S)|. 1.1
1.1
1
1
1
0.9
0.9
0.9
0.8
0.8
0.8
0.7
0.7
0.7
0.6
0.6
0.6
0.5
0.5
0.5
0.4 0
10
20
(a)
30
40
0.4 0
1.1
50
100
(b)
150
200
0.4 0
100
200
300
400
500
600
(c)
Figure 5: For r = 0.8 and m = 3 varying n, ◦ are the eigenvalues of S while × denote R(f ): (a) n = 20, (b) n = 100, (c) n = 300. λmax (f ) for less than 10−14 , while λmin (f ) and the minimal eigenvalue of S = Sn are equal up to the machine precision (for m = 1 the error is only relatively big since all the involved quantities are of the order of the machine precision, see Figures 4 (a)–(b)). Figure 5 shows that the eigenvalues of S are all contained in R(f ) with the only exception of (O(log(n))) that are far away from the 4 points forming R(f ), but still in the envelope of the range (λmin (f ), λmax (f )); see [11] and [4, p. 143] for the study of a similar situation in the case of a scalar-valued symbol. Figure 6 shows an interesting asymptotic behavior of R(f ) as the parameter m tends to infinity: in the limit two of the points forming R(f ) collapse from above to r2 and the other two collapse from below to 2r − r2 according to (43) and (44), respectively. The last test shows how fast the minimum eigenvalue of S approaches λmin (f ). Figure 7 shows that increasing n, λmin (S) converges exponentially to λmin (f ): see conjectures s1 and s2. We finally remark that the exponential convergence (with difference going to zero monotonically as e−n ) is the fastest rate we can have for block Toeplitz matrices, according to the result conjectured in [13] and formally proved by Tilli in [19].
6
Conclusions and open problems
We have considered a special type of signal restoration problem where some of the sampling data are not available. The formulation leads to a large nonsymmetric block Toeplitz matrix which can be associated 20
1.1
1.1
1.1
1
1
1
0.9
0.9
0.9
0.8
0.8
0.8
0.7
0.7
0.7
0.6
0.6
0.6
0.5 0
200
400
600
800
1000
0.5 0
200
(a)
400
600
800
1000
0.5 0
200
400
(b)
600
800
1000
40
50
(c)
Figure 6: R(f ) for r = 0.8 varying m: (a) m = 4, (b) m = 31, (c) m = 101.
0
0
10
−5
10
−5
10
−5
10
−10
−10
10
−15
10
−15
10
−15
10
−20
10
−20
10
−20
10
−25
0
10
−10
10
10
0
10
10
−25
10
20
30
r = 0.7, m = 9
40
50
10
0
−25
10
20
30
r = 0.9, m = 4
40
50
10
0
10
20
30
r = 0.8, m = 11
Figure 7: The dashed curve depicts the reference function e−n , while the solid curve depicts |λmin (S) − λmin (f )|/|λmin (f )|, versus n = 1, . . . , 50.
21
to a 2 × 2 matrix-valued symbol f (θ) that has been identified and studied. The use of recent results on the eigenvalue localization and distribution of block Toeplitz matrix-sequences, has been of great help in describing the cluster sets, the spectral distribution, and the localization areas/extremal behavior of the spectrum of the double window matrix, although the information is not yet complete. In any case, we believe that further developments in the asymptotic theory for nonsymmetric symbols will be possible, providing a very effective tool to predict the full dependence of the spectral behaviour on the parameters. In addition, the applicative impact of such future results could address the most appropriate choice of free parameters in the design of transmission devices. Furthermore, the theory is very general and so we are confident that our analysis can be extended along the following lines: • The case of arbitrary positions of missing samples (one difficulty is given by the lack of spectral results relating a non-symmetric matrix to its submatrices). The more specific case where the positions of the missing samples are approximately described by a nonlinear function will lead to Generalized Locally Toeplitz structures whose spectral theory is well developed (see [17, 14]). • The case of s > 2 channels (leading to s × s matrix-valued symbols). • The case of 2D or 3D images in place of a one-dimensional signal F (leading to 2- or 3-level Toeplitz structures with symbol depending on 2 or 3 variables). This new reconstruction problem presents some analogies with the inpainting problem, investigated with the help of structured matrices by R. Chan and collaborators in a series of papers (see e.g. [5]).
References [1] P. Brianzi, V. Del Prete, Stability of the Recovery of Missing Samples in Derivative Oversampling, Sampl. Theory in Signal and Image Process. 10 (2011) 191–210. [2] A. B¨ ottcher, S. Grudsky, On the condition numbers of large semi-definite Toeplitz matrices, Linear Algebra Appl. 279 (1998) 285–301. [3] A. B¨ ottcher, S. Grudsky, E. Ramirez de Arellano, On the asymptotic behavior of the eigenvectors of large banded Toeplitz Matrices, Math. Nachr. 279 (2006) 121–129. [4] A. B¨ ottcher, B. Silbermann, Introduction to Large Truncated Toeplitz Matrices, Springer, New York, 1999. [5] J.F. Cai, R.H. Chan, Z.W. Shen, A Framelet-Based Image Inpainting Algorithm, Appl. Comput. Harmon. Anal. 24 (2008) 131–149. [6] M. Donatelli, M. Neytcheva, S. Serra-Capizzano, Canonical eigenvalue distribution of multilevel block Toeplitz sequences with non-Hermitian symbols, Operator Theory: Adv. and Appl. 221 (2012) 273–295. [7] M. Donatelli, S. Serra-Capizzano, D. Sesana, Multigrid methods for Toeplitz linear systems with different size reduction, BIT 52 (2012) 305–327. [8] P.J.S.G. Ferreira, The stability of a procedure for the recovery of lost samples in band-limited signals, Signal Processing 40 (1994) 195–205. [9] W. Rudin, Real and Complex Analysis, McGraw-Hill, New York, 1974. [10] D.M.S. Santos, P.J.S.G. Ferreira, Reconstruction from Missing Function and Derivative Samples and Oversampled Filter Banks, Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP 3 (2004) 941–944. 22
[11] S. Serra-Capizzano, Preconditioning strategies for Hermitian Toeplitz systems with nondefinite generating functions, SIAM J. Matrix Anal. Appl. 17 (1996) 1007–1019. [12] S. Serra-Capizzano, On the extreme eigenvalues of Hermitian (block) Toeplitz matrices, Linear Algebra Appl. 270 (1998) 109–129. [13] S. Serra-Capizzano, How bad can positive definite Toeplitz matrices be?, Numer. Funct. Anal. Optim., 21 (2000) 255–261. [14] S. Serra-Capizzano, Generalized Locally Toeplitz sequences: spectral analysis and applications to discretized Partial Differential Equations, Linear Algebra Appl. 366 (2003) 371–402. [15] S. Serra-Capizzano, P. Tilli, Extreme singular values and eigenvalues of non Hermitian block Toeplitz matrices, J. Comput. Appl. Math. 108 (1999) 113–130. [16] S. Serra-Capizzano, Spectral and computational analysis of block Toeplitz matrices having nonnegative definite matrix-valued generating functions, BIT 39 (1999) 152–175. [17] P. Tilli, Locally Toeplitz sequences: spectral properties and applications, Linear Algebra Appl. 278 (1998) 91–120. [18] P. Tilli, A note on the spectral distribution of Toeplitz matrices, Linear Multilin. Algebra 45 (1998) 147–159. [19] P. Tilli, Universal bounds on the convergence rate of extreme Toeplitz eigenvalues, Linear Algebra Appl. 366 (2003) 403–416. [20] J.M. Varah, The prolate matrix, Linear Algebra Appl. 187 (1993) 269–278. [21] N.L. Zamarashkin, E.E. Tyrtyshnikov, On the distribution of eigenvectors of Toeplitz matrices with weakened requirements on the generating function, Russ. Math. Surv. 522 (1997) 1333–1334.
23