DOI 10.1007/s00498-015-0141-6 O R I G I NA L A RT I C L E
A sampled-data approach to optimal non-causal
downsampling
Hanumant Singh Shekhawat1 · Gjerrit Meinsma1
Received: 10 April 2013 / Accepted: 29 March 2015 / Published online: 26 June 2015 © Springer-Verlag London 2015
Abstract Downsampling is the process of reducing the sampling rate of a discrete signal. This paper describes how sampled data system theory can be used to design an L2or L∞optimal downsampler which reduces the sampling rate by a positive integer factor from a given fast sampler. This paper also describes the effect of noise on the optimal downsampling.
Keywords Causality· Linearity · Sampling · Shannon’s sampling theorem · System norms· Time-shift invariance
1 Introduction
In most of the signal processing applications, information is taken from analog sources and is utilized in an analog format. Due to the high quality of digital transmission, information is often converted into a discrete form at the front end from an analog source and converted back to analog at the rear end. Sometimes in applications like image, audio, and video, it is necessary to change the sampling rate of the discrete signal. There are several reasons for this sampling rate conversion, such as different clocks in two digital hardwares. The process of increasing the sampling rate is known as upsampling or interpolation whereas decreasing the sampling rate is known as
downsampling or decimation [5]. The increment or reduction factor of the sampling rate can be a positive integer or rational number. In this paper, we concentrate on
B
Hanumant Singh Shekhawat h.s.shekhawat@alumnus.utwente.nl Gjerrit Meinsmag.meinsma@utwente.nl
downsampling with a reduction factor M (M is a positive integer). Instead of imposing any structure, we just assume that the downsampler is linear and M-shift invariant.
One approach to the problem of reconstruction of the analog signal from the discrete signal uses the theory of sampled-data systems [4,8,12,13,25]. This theory was first applied in the signal reconstruction problem in 1996 by Khargonekar and Yamamoto [8] (in 1995, Chen and Francis [3] applied the sampled-data system theory to the signal reconstruction problem entirely in the discrete domain). A lot of work has been done in this area, see [25] for a detailed discussion and applications. See also [10]. The two distinctive features of the sampled-data approach are the use of signal generation models and the use of system norms for the measurement of reconstruction performance [11]. The sampled-data system theory can be applied to downsampling problems also [7]. The sampled-data setup for the downsampling problem is shown in Fig.1. Here, the continuous time signal y is sampled at interval hby a given sampler or A/D converterSh, then the sampling rate of the resulting discrete signal ¯yhis changed by a downsampler ¯Sh to a slower rate of 1h := Mh1, where M is a positive integer
known as the downsampling factor. The output ¯yhof the downsampler is converted back to the analog domain by a hold or D/A converterH. Any discrete filter in between ¯ShandH can be combined in either ¯ShorH without loss of generality. The continuous time signals y andv are modeled as the output of a known linear continuous time-invariant (LCTI) systemG driven by a known set of signals w (see Fig.1).G can be thought as a model of the frequency characteristics of the signals y andv. The choice ofG depends upon the application at hand. The main aim in downsampling problem is to reconstruct u as close as possible tov with design parameters ¯ShandH, given
G and Sh. In most of the applications,v is the same as y but for a generic treatment of the downsampling problemv may be taken different from y. Also normalization of the signalw allows us to write the reconstruction performance in terms of the norm of the systemGethat mapsw to e in Fig.1[1,2,11]. Typically L2(or H2) and L∞(or H∞) system norms are used for measuring the reconstruction performance [11]. The advantage of using sampled-data system theory for downsampling is that it utilizes the inter-sample information.
A general approach for downsampling in the signal processing literature is to some-how band-limit the input signal y to the Nyquist frequencyωn := π/h rad/sec (or
any multiple of the Nyquist frequency) to prevent aliasing errors. However, we lose a lot of information by bandlimiting our input signal y as most of the signals in prac-tice are rarely bandlimited. To utilize the information available maximally, we use sampled-data system theory to solve the downsampling problem. The downsampling problem from the sampled-data viewpoint is studied first by Ishii et al. [7], where
Ge e ¯yh ¯yh y υ u H S¯h Sh G
they use the H∞error criteria to design causal optimal downsamplers ¯Sh (givenH,
Sh andG). Nagahara and Yamamoto [14] designed computationally efficient H∞ -optimal downsamplers (givenH, Sh andG). Meinsma and Mirkin [12] solved the problem of designing non-casual downsamplers (givenH, Sh andG) and the prob-lem of designing non-casual holds (given ¯Sh,Sh andG). The problem of designing optimal non-causal downsamplers and holds simultaneously (givenShandG) is first studied in [12], but for a very limited class of the signal generatorsG using L2 and L∞system norm for measuring the reconstruction performance.
The main contribution of this paper is to extend the downsampling result of [12] for all LCTI signal generatorsG using sampled-data system theory. This leads to a generic system theoretical treatment of the downsampling problem. Specifically, we consider the problem of designing both a non-causal downsampler and a hold, given
Sh andG in this paper. We also quantify the signal reconstruction error and show its variation in the presence of noise.
The outline of this paper is as follows. First, we describe mathematical representa-tions of the sampled-data setup for the downsampling problem in the Sect.2. Basics of the lifting technique and norms are summarized in the Sect.3. An important the-orem known as Rank thethe-orem and mathematical details of downsampling problem are discussed in the Sects.4and5. Sections6and 7contain the main results of this paper, i.e. the solution of the downsampling problem with different error norms along with some examples and an expression for the error norm. The effect of noise on the downsampling problem is discussed in the Sect.8.
1.1 Notation
Linear discrete time invariant (LDTI) means a linear h-time shift-invariant system. A system in the time domain will be represented in calligraphic, e.g.G and the corre-sponding classic frequency response is represented in uppercase like G(jω). Signals will be denoted in lowercase letters like y and discrete signal is represented by bar on top, e.g. ¯y. A discrete signal ¯y with interval hmeans that the samples of the discrete signal are htime unit apart. Square brackets are used to denote the value of a discrete signal at a given integer, e.g. ¯y[k] whereas parenthesis is used to denote the value of an analog signal at a given time, e.g. y(t).
The kernel and image of an operator A are denoted as Ker A and Im A, respectively. Im(A) denotes the closure of the image of an operator A and span{x1, x2, . . .} denotes
the closure of the span of vectors x1, x2,. . .. · H S and · denote the Hilbert– Schmidt norm and the induced 2-norm of an operator, respectively. rank(F) denotes the rank of the operator F . SVD stands for singular value decomposition. Also the left singular vectors of an operator A mean the eigenvectors of A A∗and the right singular vectors of an operator A mean the eigenvectors of A∗A.
Z represents the set of integers and Z+(Z−) represents the set of positive (negative)
integers.N represents the set Z+∪{0}. For any M ∈ Z+,M represents the set of integers {0, 1, . . . , M − 1}.
Throughoutωk :=θ+2πkh for any givenθ ∈ [−π, π] and all k ∈ Z. We also define sinc(x) := sinπx(πx)and j:=√−1. ¯δ[i] is the discrete unit pulse. The indicator function 1S(t) is identity when t ∈ S and 0 when t /∈ S.
L2(B) is the Hilbert space of square integrable functions f : B → C where
B ⊆ R. The inner product in this space is x, y = Bx(τ)y(τ)∗ dτ. Whenever
B = R, we just write L2.
2(B, Cn) is the Hilbert space of sequences mapping B ⊆ Z to n-dimensional vectors inCn, i.e. 2(B, Cn) := { f : B → Cn|i∈B fi2Cn < ∞}. Here, f :=
{ fi}i∈B. The inner product in this space is x, y =
i∈B xi, yi Cn. We use the
shorthand2(B) to denote 2(B, Cn) if n is unspecified or unambiguous. Moreover, ifB = Z, we just write 2.
2 Sampled-data setup
In this section, we review standard mathematical descriptions of all components of the sampled-data setup given in Fig.1(for detail see [11]). The sampled-data setup for the downsampling problem is again shown in Fig.2whereG =
Gv
Gy
is an LCTI system andGvandGyare the partition ofG with respect to the signals v and y, respectively.
The samplerSh is a device which samples an analog signal y : R → C at every integer multiple of hand gives a discrete signal ¯yh : Z → C. We assume Sh to be linear and h-time shift invariant, i.e.
¯yh = Shy: ¯yh[n] =
∞
−∞ψ(nh
− s)y(s) ds (1)
whereψ(t) is known as the sampling function.
The downsampler ¯Shis a device which operates on a discrete signal ¯yh : Z → C with interval hand produces a discrete signal ¯yh : Z → Cr with interval h = Mh at the output. Here, M is called downsampling factor and r is a positive integer. We assume the downsampler to be linear and M shift invariant, i.e.
¯yh= ¯Sh¯yh : ¯yh[n] = k∈Z
¯χ[Mn − k] ¯yh[k], n ∈ Z. (2) Here, ¯χ : Z → Cr×1is known as the discrete sampling function.
The holdH is a device which converts a discrete signal ¯yh : Z → Cr with interval
h back to an analog signal u : R → C. Here, r is a positive integer. We assume H to
be linear and h-time shift invariant, i.e.
F Sy y e ¯yh ¯yh h u H S¯h S Gv Gy -υ
u= H ¯yh: u(t) = n∈Z
φ(t − nh) ¯yh[n], t ∈ R. (3)
Theφ : R → C1×r is known as the hold function.
Since the output of the downsampler is the input of the hold, the output dimension of the downsampler and the input dimension of the hold must be same.
In this paper, we use interpolator to mean the following.
Definition 1 (Interpolator) An interpolatorF is a device which converts a discrete signal ¯yh : Z → C with interval hback to an analog signal u: R → C.
Note that the interpolator’s input is a discrete signal with interval hwhereas the hold’s input is a discrete signal with interval h. A downsampler followed by a hold (i.e.H ¯Sh) is an important component throughout this paper, so it deserves a special name: Definition 2 (Hybrid interpolator) An interpolatorF is a hybrid interpolator of order
r if it can be represented as a downsampler (with output dimension r ) followed by a
hold (with input dimension r ), i.e.F := H ¯Sh.
3 Lifting and transforms
Lifting was first introduced in [22] (see also [23]) for control and tracking applications. Lifting is also used for solving several signal reconstruction problems in [8,12,13] (see [4,11], for example, for more expositions). This section summarizes the lifting definition used in [11] and the lifting of a discrete signal with a sampling interval different from h [24]. We also review the meaning of lifted transfer functions and norms in this section. We start with lifting.
Definition 3 (Lifting) For a continuous time signal f : R → Cn, the lifted signal ˘f : Z → {[0, h) → Cn} is the sequence of functions { ˘f[k]} defined as
˘f[k](τ) := f (kh + τ), k ∈ Z, τ ∈ [0, h). Here,τ is inter-sample time parameter.
Sometimes it is beneficial to define ˘f[k](τ) := f (kh + τ) for arbitrary τ ∈ R. The lifted z-transform of an analog signal f is defined as
˘f(z; τ) := k∈Z
˘f[k](τ)z−k= k∈Z
f(kh + τ)z−k. (4)
We can obtain the lifted Fourier transform by substituting z= ejθin (4).
The following, which is known as the key lifting formula (see [11, §IV.A] for details), shows the bijection between the lifted Fourier transform and the classical Fourier transform:
˘f(ejθ; τ) = 1 h k∈Z F(jωk)ejωkτ (5a) F(jωk) = h 0 ˘f(e jθ; τ)e−jωkτ dτ (5b) whereωk := θ+2πkh . Here, F(jω) := ∞
−∞ f(t)e−jωtdt is the classical Fourier
trans-form and ˘f(ejθ; τ) represents the lifted Fourier transform of the continuous time signal
f : R → Cn.
Similar to lifting of a continuous time signal, lifting of a discrete signal with sam-pling interval h can be defined as follows (see e.g. [24] and the references cited therein).
Definition 4 (Discrete lifting) For any discrete signal ¯f : Z → Cnf sampled at
interval h = h/M, the lifting (with respect to h) f : Z → {M → Cnf} is the
sequence of discrete functions{ f[k]} defined as
f[k; m] = ¯f[Mk + m], k ∈ Z, m ∈ M. Here, nf ∈ Z+andM := {0, 1, . . . , M − 1}.
Sometimes, it is beneficial to define ˘f[k; m] := f [Mk + m] for arbitrary m ∈ Z.
The discrete time lifting is the same as the concept of polyphase decomposition (see, e.g. [21]). Similar to analog signals we can define the lifted z-transform for a discrete signal. The lifted z-transform f(z; m) of a discrete signal ¯f is defined as
f(z; m) := k∈Z
f[k; m]z−k = k∈Z
¯f[Mk + m]z−k (6)
where m∈ M. Similarly, the lifted discrete Fourier transform (LDFT) f(ejθ; m) of a
discrete signal ¯f is defined as
f(ejθ; m) := k∈Z f[k; m]e−jθk = k∈Z ¯f[Mk + m]e−jθk (7) where m∈ M.
Since lifting and z-transform are invertible processes, the inverse lifted z-transform can be defined as a combination of the inverse z-transform and the inverse of the lifting process.
The following is a key result (similar to [21, Eqs. (10), (11) and (6)]), which relates the lifted z-transform of a discrete signal f[n] with its classic z-transform.
Theorem 1 (Discrete key lifting formula) Consider a discrete signal ¯f : Z → C with time interval hand let ¯f(z) be its classic z-transform. If ¯f(z) exists then the lifted z-transform f(z; m) of ¯f[n] follows from its classic z-transform as
f(z; m) = 1 M M−1 i=0 ¯f(˜zi)˜zim m∈ M (8a) where˜zi = |z| 1 Mej Arg z M e j2πi
M (Arg is the principle argument function [19, chap3]) are
the M complex roots of the equation ˜zM = z for a given z. Conversely, the classic z-transform of ¯f follows from the lifted z-transform as
¯f(˜z) =M−1 m=0
f(˜zM; m)˜z−m (8b)
for a given˜z.
The proof of Theorem1is given in AppendixA. Let ¯f(ejν), ν ∈ [−π, π] represents
the classic discrete time Fourier transform of the signal ¯f in Theorem1[9, chap. 10]. Then, setting z= ejθin (8a) and˜z = ejνk (ν
k:= θ+2πkM , k ∈ M) in (8b), we have f(ejθ; m) = 1 M M−1 i=0 ¯f(ejωkh)ejωkhm = 1 M M−1 i=0 ¯f(ejνk)ejνkm (9a) ¯f(ejνk) = M−1 m=0 f(ejθ; m)e−jνkm (9b)
whereωk =θ+2πkh andθ ∈ [−π, π]. Note that at a given θ, the set {√1Mejωih m
}m∈M forms an orthonormal basis of spaceCM.
3.1 Transfer function in lifted frequency domain
In this section, we define transfer functions for h-time shift-invariant systems, down-samplers and holds in lifted frequency domain. For a more detailed discussion on (lifted) transfer function, we suggest [11].
Any linear system G that is h-time shift invariant in continuous time is by con-struction linear shift-invariant system in the lifted domain with respect to the discrete variable and hence can be written as a convolution
u= Gylifting⇒ ˘u[k] =
i
˘G[k − i] ˘y[i].
Here, ˘G[k − i] ˘y[i] for each i is a finite integral over inter-sample time [11]. Taking lifted Fourier transforms, we have
˘u(ejθ) = ˘G(ejθ) ˘y(ejθ),
where ˘G(ejθ) := i ˘G[i]e−jθi [11, §IV.B]. ˘G(ejθ) is known as the (lifted) transfer function ofG.
Now, the kernel of the operator ¯Shdefined in (2) in lifted frequency domain can be found as ¯yh(ejθ) = n∈Z ¯yh[n]e−jθn= n∈Z k∈Z χ[Mn − k] ¯yh[k]e−jθn = n∈Z l∈Z M−1 m=0 χ[Mn − Ml − m] ¯yh[Ml + m]e−jθn = M−1 m=0 χ(ejθ; −m) y h(ejθ; m) =: ´Sh(ejθ) yh(ejθ). (10)
Therefore, the kernel of the operator ¯Sh in the lifted frequency domain is given by χ(ejθ, −m). The (lifted) transfer function of ¯S
his denoted by ´Sh(ejθ).
Similarly, the (lifted) transfer function `H(ejθ) of the hold H defined in (3) is given by [11]
˘u(ejθ) := `H(ejθ) ¯y
h(ejθ) : ˘u(ejθ; τ) = ˘φ(ejθ; τ) ¯yh(ejθ).
3.2 Lifted norms and stability
In this section, we review the definition of some standard norms in sampled-data system theory. Norms defined in this section are very standard and discussed in great detail in [1,2,11]. In this section, systemG means any LDTI system (including LCTI system, sampler, downsampler and hold). The (lifted) transfer function is denoted by ˜G(ejθ).
L2is the space of LDTI systems with finite norm defined as [2,11]:
GL2 := 1 2πh π −π ˜G(e jθ)2 H Sdθ (11)
where.H Sstands for Hilbert–Schmidt norm of an operator [26, §8.1]. L∞is the space of LDTI systems with finite norm defined as [1,11]:
GL∞:= ess sup θ∈[−π,π] ˜G(e
jθ)
(12)
where ˜G(ejθ) = sup ˜x
2=1, ˜x∈˜S ˜G(e jθ)x
2. Here, ˜S means either L2[0, h) or Cn
or2(M, Cn) depending upon whether ˜x is lifted from an analog signal or a discrete signal (with interval h) or a discrete signal with interval h, respectively. The above definition of the L∞norm is equivalent to the induced norm interpretation given as GL∞ := supxSi =1GxSowhereSiandSocan be L2(R) or 2(Z) depending upon whether the signal in concern is analog or discrete, respectively [11, §V.C].
Definition 5 (Stability) An LDTI systemG : Si → So, whereSiandSocan be L2(R) or2(Z), is stable if G ∈ L∞.
The above definition says that an analog system is stable if it a bounded operator from L2(R) to L2(R), a discrete system (or downsampler) is stable if it is a bounded operator from 2(Z) to 2(Z), a sampler is stable if it is a bounded operator from
L2(R) to 2(Z), and a hold is stable if it is a bounded operator from 2(Z) to L2(R)
[11, theorem 6.1].
4 Rank theorem
In time domain, it is not so clear when a given interpolator is a hybrid interpolator of a given order r . However, it is somewhat clear in the (lifted) frequency domain as shown in this section. The proofs of the results in this section are given in AppendixA. First, we write the hybrid interpolatorF in the lifted frequency domain.
Lemma 1 The hybrid interpolatorF := H ¯Shin lifted frequency domain is an
oper-ator of the form
˘u = `F yh : ˘u(ejθ; τ) = M−1 m=0
˘ (ejθ; τ; m) y
h(ejθ; m)
where the kernel can be expressed in terms of the sampling and hold functions (see Sect.2) as
˘ (ejθ; τ; m) = ˘φ(ejθ; τ) χ(ejθ; −m).
(13) The above lemma is useful if we have an interpolator `F and want to represent it as
a cascade of a downsampler and a hold, i.e. as a hybrid interpolator.
Now we concentrate on obtaining a condition for an interpolator such that it is a hybrid interpolator. Consider an order r hybrid interpolator F shown in Fig.3. In
F Sy e y ¯y h ¯yh,1 ¯yh,0 ¯yh,r−1 u H1 S¯h, 1 Sh Gv Gy - H0 S¯h, 0 Hr−1 S¯h,r−1 υ
this case, the downsampling function ¯χ[n] of downsampler ¯Sh has r rows. Clearly, the lifted Fourier transform ¯yh(ejθ) of the output ¯yhof the downsampler has r rows. Therefore for allθ ∈ [−π, π], we have ¯yh(ejθ) ∈ Cr. Hence, rank ´Sh(ejθ) ≤ r for allθ. We can obtain such an upper bound for any hybrid interpolator F := H ¯Sh as shown in the following lemma.
Theorem 2 (Rank theorem) Given an interpolatorF and assume that the transfer
function `F(ejθ) exists. If F is a hybrid interpolator of order r, then rank `F(ejθ) ≤ r∀θ ∈ [−π, π].
Since the input of an interpolator is a discrete signal, we can expect that the rank of the interpolator in the lifted domain must have an upper bound (at a givenθ). Indeed. Lemma 2 IfF is an interpolator then rank `F(ejθ) ≤ M for all θ ∈ [−π, π].
IfF is a hybrid interpolator then rank `F(ejθ) ≤ M for all θ ∈ [−π, π] by Lemma2. Hence, increasing the order r of a hybrid interpolator beyond M is unnecessary.
5 Downsampling problem
In this section, we define the downsampling problem and simplify it as much as possible using lifting. The proofs of the results in this section are given in AppendixA. In this section,K represents either L2or L∞and · Krepresents the Hilbert–Schmidt (H S)
norm ifK = L2or the induced 2-norm ifK = L∞. Throughout this paper, we use short-handSy:= ShGy. Clearly,Syis a sampler sampling at interval h.
Now, we state the downsampling problem more precisely.
Problem 1 (Downsampling problem) GivenGv,Sy∈ L2∩ L∞, find a hybrid
interpo-latorF ∈ L∞of at most order r ∈ Z+such thatGv− FSyKis minimized.
Note that ifSy∈ L∞thenSy∈ L2as it is a sampler [11, prop.5.3]. Theorem2says
that an interpolatorF is a hybrid interpolator of order r only if rank( `F) ≤ r at all
θ ∈ [−π, π]. This result is used to solve Problem 1by translating it to an almost equivalent lifted frequency domain problem as shown below.
Theorem 3 GivenGv, Sy∈ L∞∩ L2. Define `Fopt(ejθ) at almost every θ ∈ [−π, π]
as
`Fopt(ejθ) := arg min `F(ejθ)
˘Gv(ejθ) − `F(ejθ) ´Sy(ejθ)K (14)
with constraint rank `F(ejθ) ≤ r. If interpolator Foptis well defined and stable, then it
minimizesGv− FSyKover all interpolatorsF with rank `F(ejθ) ≤ r at almost all
θ.
Remark 1 Note that we can always redefine Fopt in (14) by adding a null set in
chosen in such a way that rank `F(ejθ) ≤ r for all θ ∈ [−π, π]. This is a necessary
condition for a hybrid interpolator (see Theorem 2). Therefore, from now on we concentrate on interpolatorsF with constraint rank `F(ejθ) ≤ r at almost all θ ∈ [−π, π].
Remark 2 Note that the operators ´Sy(ejθ):L2[0, h) → CM and ˘Gv(ejθ) : L2[0, h) →
L2[0, h) at almost all θ ∈ [−π, π] are Hilbert–Schmidt (hence compact) operators as GvandSybelong to L2(see [20, lemma 2.4.5]).
If interpolatorFoptin (14) is also a (well defined) hybrid interpolator, then
Theo-rem3says that doing point-wise minimization in lifted Fourier domain is sufficient to solve Problem1. Therefore, to solve the downsampling problem, we consider the following problem.
Problem 2 GivenGv, Sy ∈ L2∩ L∞, find a well-defined and stable interpolatorF
with rank `F(ejθ) ≤ r such that ˘Gv(ejθ) − `F(ejθ) ´Sy(ejθ)Kis minimized at almost
eachθ ∈ [−π, π].
To solve Problem1, we first obtain a solution to Problem2using Theorem3. Then, we write the resulting optimal interpolator as a (well defined) hybrid interpolator. This is straight forward as we will see later.
In the rest of this paper, a transfer function ˘G(ejθ) at given θ is abbreviated as ˘G. We do the same for samplers, downsamplers and holds also. Unless necessary, we do this for all signals too.
6 L
2optimal downsampling
The L2optimal downsampling problem is the Problem1withK = L2. It can be solved by obtaining a solution of Problem2 with Hilbert–Schmidt norm. Section6.1 con-tains some basic results which are important in solving the L2optimal downsampling problem in Sect.6.2. The proofs of the results in this section are given in AppendixB. 6.1 Preliminaries
This section describes some basic results which are useful in obtaining a solution of Problem2with the Hilbert–Schmidt norm, i.e. obtaining
`Fopt:= arg min
`F ˘Gv− `F ´SyH S (15)
with constraint rank `Fopt ≤ r at each θ ∈ [−π, π]. Here, the minimization is
point-wise inθ. At almost all θ, the operators ˘Gv and ´Syare Hilbert–Schmidt operators
(see Remark2). Now, Problem2can be solved by considering a generic problem of obtaining
Fopt:= arg min F∈Fr
where A: H → H0and B : H → H1are Hilbert–Schmidt operators, F: H1→ H0
is a linear operator (not necessarily bounded) andFr denotes the set of all bounded linear operators mapping the spaceH1toH0and that have rank at most r . The spaces
H, H0andH1 are assumed to be separable Hilbert spaces. In our special case of
Problem2,H = H0= L2[0, h) and H1= CM, which are indeed separable.
Also, let P : H → H represent the orthogonal projection on (Ker B)⊥throughout this section. We also define the set
NG :=
{0, 1, . . . , rank G − 1} if rank G is finite
N if rank G is infinite (16)
for a linear operator G.
An algorithm for obtaining Fopt:= arg minF∈FrA − FBH Sis presented at the end of this section.
Since(Ker AP)⊥⊆ (Ker B)⊥, the orthogonal projection P helps us in identifying the space where F could play a role. Using B(I − P) = 0, we have the following result.
Lemma 3
A − FB2
H S= AN2H S+ AP − FB2H S. (17)
where AN := A(I − P).
From (17) it is clear that an F∈ Fr that minimizesA − FBH Salso minimizes AP − FBH S. If a rank-r F is such that F B cancels out the dominant r singular values in an SVD of A P, then this F is optimal. This fact is used in the following theorem to obtain Fopt.
Theorem 4 Suppose an SVD of A P is given by
A P x=
i∈NA P
σi x, ei fi, x ∈ H (18)
where{ei} and { fi} are orthonormal sequences in H and H0respectively, and{σi} is a
real non-increasing positive sequence. If Im B is closed then for a given r ≤ rank(AP),
Foptx:= r−1 i=0 σi B+x, ei fi, x ∈ H1 (19a)
is bounded and it minimizesAP − FBH Sover all F∈ Fr. Here B+: H1→ H is
the pseudo-inverse of the operator B [6, Definition 2.2]. Also, minF∈FrAP −FBH S
is given by
Equation (20) may be surprising because, even though Fr is not a subspace of the space of bounded linear operators, we still have a Pythagoras type result.
Remark 3 The image of a compact operator mapping two Hilbert spaces is closed if
and only if the operator is of finite rank (it follows from [17, Theorem 4.18] and [18, Lemma 1.10]). In the downsampling problems considered in this paper, the operator
B is a finite rank compact operator (hence, Im B is closed). For a general form of
Theorem4valid for both closed and not closed Im B, see [20, Theorem 4.9.2]. An example of not closed Im B (in other words, of infinite rank compact operator B) is discussed in [12, Theorem 5.1].
Now we can summarize the steps for obtaining Fopt:= arg minF∈Fr A − FBH Sin
the form of Algorithm1.
Algorithm 1 For Fopt:= arg minF∈Fr A − FBH S Require: A and B are Hilbert–Schmidt operators and B is of finite rank.
Obtain the orthogonal projection P on(Ker B)⊥ if A P= 0 then
Fopt= 0
else
Obtain SVD of A P. Obtain B+
Foptis given by (19a).
end if
CalculateA − FoptBH Susing (20)
return Fopt,A − FoptBH S
6.2 L2optimal downsampling solution
Recall thatSy:= ShGy. The Problem2with Hilbert–Schmidt norm can be solved by
application of Algorithm1. Among the important parameters to find in application of Algorithm1are the orthogonal projection ˘P on the space(Ker ´Sy)⊥, the pseudoinverse
´S+
y of ´Syand an SVD of the operator ˘Gv ˘P at almost every θ ∈ [−π, π]. To obtain ´S+y
and ˘P, we compute an SVD of ´Sy.
Lemma 4 GivenSy:= ShGy∈ L∞∩ L2. Then,
1. Syis a sampler with sampling functionψy(t) := ψ(t) ∗ gy(t). Here ∗ represents
the analog convolution operator and gyis the impulse response ofGy.
2. For almost all givenθ ∈ [−π, π], the lifted transfer function ´Syof the sampler
Syis a mapping from L2[0, h) to CM and it has an SVD (modulo ordering) of the
form ´Sy˘y = M h k∈M αk ˘y, ˘pk ek, (21)
where, for fixedθ, the αk ∈ C, ˘ek ∈ L2[0, h), ˘pk ∈ L2[0, h) and ek ∈ CM are defined as αk:= i∈Z |ψy(jωk+Mi)|2, ˘ek(ejθ; τ) := 1 √ he jωkτ, ˘pk(ejθ; σ ) := i∈Z ψ∗ y(jωk+Mi) αk ˘e k+Mi(ejθ; σ), ek(ejθ; m) := 1 √ Me jωkmh, m ∈ M. Hereτ, σ ∈ [0, h) and ωk= θ+2πkh , k ∈ Z.
Remark 4 Note that the operator ´Sy(ejθ) is a finite rank operator at almost each θ.
Hence, Im ´Sy(ejθ) is closed at almost each θ [6, p.37].
Remark 5 Sometimes direct calculation ofαk (defined in Lemma4) is not obvious. As an alternative, we can use the approach given below.
Let us defineψy(2)(t) := ψy(t) ∗ ψy(−t) where ∗ is the convolution operator and
the overline denotes the complex conjugation. Now, the classical Fourier transform of
ψy(2)(t) is |ψy(jω)|2. If the signalψy(2)(t) is ideally sampled at regular interval hthen
we have that ψy(2)(ejν) = 1 h i∈Z |ψy(jωi)| 2. (22)
Hereωi := ν+2πih for allν ∈ [−π, π] and ψy(2)(ejν) denotes the classical discrete
time Fourier transform of the discrete signalψy(2)[n] := ψy(2)(nh). The equivalence
in (22) is a standard result (see for example [16, chap 4]) and it can be also obtained by the key lifting formula (5) (by substitutingτ = 0).
Now substitutingν = θ+2πkM whereθ ∈ [−π, π] and k ∈ M, we have that
ωi = ν + 2πi h = θ + 2π(k + Mi) h = jωk+Mi. Here h= Mh. Now, ψ(2) y ejθ+2πkM = 1 h i∈Z ψy jθ + 2π(k + Mi h 2= 1 h i∈Z |ψy(jωk+Mi)|2. So, α2 k = hψy(2)(e jθ+2πkM ). (23)
Using (21), the pseudoinverse ´Sy+(of ´Syat a givenθ) mapping CMto L2[0, h) is given by ´S+ y yh = h M k∈M,αk=0 1 αk yh , ek ˘pk. (24) Having obtained an SVD of the operator ´Syin Lemma4, we can immediately write
down an orthonormal basis{ ˘pk}k∈M,αk=0of the space(Ker ´Sy)⊥, hence ˘P at almost
everyθ ∈ [−π, π] is
˘Px =
k∈M,αk=0
x, ˘pk ˘pk, x ∈ L2[0, h). (25)
The next step in Algorithm1is to obtain an SVD of the operator ˘Gv˘P at almost every
θ ∈ [−π, π].
Lemma 5 Given Gv, Sy ∈ L∞ ∩ L2. Using all the notations and conditions of
Lemma4and (25), an SVD (modulo ordering) of the operator ˘Gv ˘P : L2[0, h) →
L2[0, h) at almost all θ ∈ [−π, π] is given by
˘Gv ˘P ˘w = n∈M αn,σn=0 σn ˘w, ˘pn ˘qn (26) where σn= i∈Z|Gv(jωn+Mi)ψy(jωn+Mi)|2 αn , ˘qn= 1 σnαn i∈Z
Gv(jωn+Mi)ψy∗(jωn+Mi)˘en+Mi.
Remark 6 An approach similar to Remark5used for calculation ofαk can be used forσk.
Let (2)(ejν) denotes the classical discrete time Fourier transform of discrete signal
[n] := (nh) ∗ (−nh) where (t) is the classical inverse Fourier transform of
Gv(jω)ψy(jω). Here ∗ is the convolution operator and the overline denotes the complex
conjugation. Also define
Σν(ν) := ⎧ ⎨ ⎩ (2)(ejν) ψy(2)(ejν), ψ (2) y (ejν) > 0 0, ψy(2)(ejν) = 0. (27)
Using the approach used in Remark5, we can show that i∈Z |Gv(jωn+Mi)ψy(jωn+Mi)|2= h (2)(ej θ+2πk M )
Using (23) and the above equality, we can show that
σ2 k = Σν θ + 2πk M . (28)
Also, ordering of the singular valueσk can be easily done by looking at the function
Σν(ν). This is because σk ≥ σnat a givenθ if and only if
Σν θ + 2πk M ≥ Σν θ + 2πn M .
Now, according to Algorithm1, we have all the basic tools to obtain a solution of Problem2. For properly writing the results we need the dominant index set:
Definition 6 The dominant index set Dr(S) is defined as an ordered set of the indices of the dominant r elements (if well-defined) of a bounded sequenceS := {σ0, σ1, . . .}.
For example ifS := {10, 2, 9, 4, 1} then D3(S) = {0, 2, 3} and D4(S) = {0, 2, 3, 1}.
Note that r in Dr(S) cannot be greater than cardinality of the S. For the applications that we consider, the dominant index set is well defined.
The following theorem, which is one of the main results of this paper, describes a solution of Problem2(with the Hilbert–Schmidt norm).
Theorem 5 GivenGv, Sy := ShGy ∈ L∞∩ L2. LetS be the set of finite singular
values of the operator ˘Gv ˘P at almost every θ ∈ [−π, π] i.e. S := {σk}k∈M,αk=0
whereσkare defined in Lemma5. Given r ≤ M, define `Foptas
`Fopt yh := k∈Dr(S) Γk yh,√1M ek , (29) where Γk:=
i∈ZGv(jωk+Mi)ψy∗(jωk+Mi)ej(ωk+Mi)τ
i∈Z|ψy(jωk+Mi)|2 ,
(30)
wheneverθ ∈ A := {θ ∈ [−π, π] : ´Sy(ejθ) = 0} and τ ∈ [0, h). For θ ∈ [−π, π]\A,
we can take `Fopt = 0. If this interpolator Fopt is well-defined and stable, then it
minimizesGv− FSyL2 over all interpolatorsF with rank `F ≤ r at almost each
1 1 √ 2 -4π -3π -2π -π 0 π 2π 3π 4π Fig. 4 Gv(jω) of Example1
Also, the L2norm of the systemGv− FoptGyis given by
Gv− FoptGy2L2 = Gv2L2− FoptGy2L2 (31) whereFoptGy2L2 = 1 2πh Ak∈Dr(S)σk2dθ.
Note that the set Dr(S) can change with θ ∈ [−π, π].
Remark 7 To solve the L2downsampling problem1, we have to also write the optimal interpolatorFopt(given in Theorem5) as a hybrid interpolator i.e.Fopt:= ¯Sh,optHopt.
To this end, (13) is very useful. Define the lifted downsampler ¯Sh,optand holdHopt
such that their (lifted) transfer functions have kernel χopt(ejθ; −m) := M1 e−jωk0mh, . . . , e−jωkr−1mh T , (32a) ˘φopt(ejθ; τ) := Γk0, . . . , Γkr−1 (32b) at almost everyθ ∈ A (defined in Theorem5) respectively. Here{k0, k1, . . . , kr−1} :=
Dr(S). If the downsampler ¯Sh,opt and hold Hopt are well-defined then Fopt :=
Hopt¯Sh,opt is a hybrid interpolator. Note that we can always redefine optimal down-sampler and hold by adding a null set in Eq. (32a) and (32b). Whenever required, we use this fact (often without mention) in the examples of this paper for simplifying the results.
For any k ∈ M, the set of frequencies {θ/h + 2ωn(k + Mi)}i∈N are called the h
-aliased frequencies of the frequencyθ/h + 2ωnk. Theorem5clearly states that these
h-aliased frequencies play a crucial role in deciding the optimal downsampler (see the construction ofS in Theorem5) for any finite M. This fact is illustrated with an example given below.
Example 1 LetGvandGybe LCTI systems with classic frequency response Gv(jω) and Gy(jω) respectively. Also, let Gv(jω) = Gy(jω) and
Gv(jω) = 1[0,π)(|ω|) + 0.951[π,3π)(|ω|) + 1 √
21[3π,4π](|ω|).
Gv(jω) is shown in Fig.4. Also, suppose that the samplerSh has sampling function
ψ(t) whose classic frequency response is given by ψ(jω) = 1[−4π,4π](ω) and M = 2,
h = 1. The Nyquist frequency is ωn = π/h = π. It can be shown that αk∀k ∈ M (see (21)) are uniformly bounded for almost allθ ∈ [−π, π]. Therefore, Syis stable
¯χopt[n] 0.5 0 φopt(t) 1 0 ¯χopt(ejν) 1 0 -π -π 2 π2 π φopt(jω) 0.5 0 2ωN -2ωN 4ωN -4ωN
Fig. 5 Discrete sampling function¯χopt[n] of the L2optimal downsampler (top left) and its classical discrete
time Fourier transform¯χopt(ejν) (bottom left) in Example1. The top left figure also shows21cos(πt) sinc(t2)
(dotted). The figure also shows the hold functionφopt(t) of the L2optimal hold (top right) and its classical
Fourier transformφopt(jω) (bottom right) in Example1. Hereωn= πh
Now, the input of the samplerShis bandlimited toωB:= 4π rad/sec as Gv(jω) is bandlimited to 4π rad/sec. Since the sampling interval of Shi.e. h(=12) is greater than
the Nyquist interval (defined as ωπ
B) required for its input i.e. 1
4, there are h-aliased
frequencies at the output of samplerSh(see Fig.2) [9, § 5.1]. Calculating the singular valuesσi(see Lemma5) in the presence of h-aliased frequencies, we find that
σ0=
5/6 ≈ 0.913, σ1= 0.95
for almost allθ ∈ [−π, π]. Note that σ1 > σ0 for almost allθ ∈ [−π, π]. Now,
it follows from Theorem 5 that if r = 1 then the L2 optimal interpolator in lifted frequency domain is given by
`Foptx=e jω1τ+ejω−1τ 2 x,12ejω1mh
for almost allθ ∈ [−π, π]. Here m ∈ {0, 1}. Using (32) and the inverse lifted trans-form, we can write the discrete sampling function ¯χopt[n] and the hold function φopt(t)
of the optimal downsampler and hold as (see also Fig.5) ¯χopt[n] = 12sinc
n
2
(−1)n,
φopt(t) = cos(2ωnt) sinc(t).
Since the optimal downsampler and hold are well-defined,Foptis an order-1 hybrid
interpolator. Using (12), it can be shown that the above optimal downsampler and hold are stable. The classical discrete time Fourier transform of ¯χopt[n] and classical
¯χopt(ejν) = 1[π 2,π](ν) + 1[−π,−π2](ν) φopt(jω) = 1 2 1[ωn,3ωn](ω) + 1[−3ωn,−ωn](ω) .
Also, we find that
Gv− FoptSy2L2 = Gv2L2 − FoptSy2L2 = 1 2πh π −π 1+12+ 2 × 0.952− 0.952 dθ = 2.4025. andGv2L2 = 3.305.
If ¯K is an invertible discrete system between the downsampler ¯Sh,optand holdHopt
then the downsampler ¯K−1¯Sh,optand holdHopt ¯K are also optimal. This is because
Hopt¯Sh,opt = Hopt ¯K ¯K−1¯Sh,opt. We use this fact in the following example, with non-bandlimited signal models, for simplifying the results.
Example 2 LetGvandGybe LCTI systems with classic frequency response Gv(jω) and Gy(jω) respectively. Also, let Gv(jω) = Gy(jω) and
Gv(jω) = 1 1+ ajω
where a is a non-zero positive real number.|Gv(jω)| (with a = 0.2) is shown in Fig.6. Furthermore, suppose that the samplerShis ideal (i.e. sampling functionψ(t) = δ(t)). Therefore,ψ(jω) = 1. Also suppose M = 2 and h = 1. The Nyquist frequency is
ωn= π/h = π. Since Gvis not bandlimited, there are h-aliased frequencies at the output of samplerSh(see Fig.2) [9, §5.1].
Next, we use the method given in Remark 5to obtain αk defined in Lemma4. Clearly,
|ψy(jω)|2= |Gv(jω)|2= 1 1+ (aω)2.
Taking inverse Fourier, we have
ψy(2)(t) = 1 2ae −|t|a . ω → -4π -3π -2π -π 0 π 2π 3π 4π 1
Next, taking the classical discrete time Fourier transform of the discrete signal ψ2 y[n] := ψy2(nh), we have ψy(2)(ejν) = (R2− 1) 2a(1 + R2− 2R cos(ν)) where R:= e h
a andν ∈ [−π, π]. Using (23), for k∈ M, we have that
α2 k = hψy(2) ej(θ+2πkM ) .
It can be shown thatαk∀k ∈ M (see (21)) are uniformly bounded for almost all
θ ∈ [−π, π]. Therefore, Sy:= ShGyis stable (see (12)).
Next, we use the method given in Remark 6 to obtain σk defined in Lemma 5. Clearly,
|Gv(jω)ψy(jω)|2= |Gv(jω)|4= 1
(1 + (aω)2)2
Taking inverse Fourier, we have
(2)(t) = 1
4a2(a + |t|)e −|t|a
Now, taking the Fourier transform of the discrete signal (2)[n] := (2)(nh), we have
(2)(ejν) = a(R2− 1)(1 + R2− 2R cos(ν)) + 2Rh((1 + R2) cos(ν) − 2R)
4a2(1 + R2− 2R cos(ν))2
whereν ∈ [−π, π]. Now, Σν(ν) defined in (27) is given by
Σν(ν) = (2)(ejν) ψy(2)(ejν) =1 2 + Rh a(R2− 1)× cos(ν) − Y 1− Y cos(ν) (33)
where Y := 2R/(1 + R2). Using (28), we have that
σ2 k = Σν θ + 2πk M . By using derivative dΣν(ν) dν = − sin(ν) Rh(R2− 1) a(1 + R2− 2R cos(v))2
and the fact R> 1, it can be shown that the function Σν is monotonically decreasing for all ν ∈ [0, π] and monotonically increasing for all ν ∈ [−π, 0]. It means that
Now, it follows from Theorem5that if r = 1 then the L2optimal interpolator in lifted frequency domain is given by
`Foptx= Γ0 x,√1 M e0 ,
for almost allθ ∈ [−π, π].
Let ¯K be a discrete system with frequency responseα20. Clearly, ¯K has a stable
inverse (use (12) and R> 1). Using (32), the lifted downsampler ¯K−1¯Sh,optand hold
Hopt¯K are such that their (lifted) transfer functions have kernel
χK opt(ejθ; −m) := 1 Mα2 0 e−jω0mh =2a 1+ R2− 2R cos θ M h(R2− 1) e −jθm M, ˘φK opt(e jθ; τ) := i∈Z Gv(jωMi)ψy∗(jωMi)ej(ωMi)τ = i∈Z |Gv(jωMi)|2ej(ωMi)τ
for almost allθ ∈ [−π, π] respectively. Using the discrete key-lifting formula (9) and the key-lifting formula (5), the classical discrete time Fourier transform of the discrete sampling function ¯χoptK [n] of the optimal downsampler ¯K−1¯Sh,optand classical Fourier transform of the hold functionφoptK (t) of the optimal hold Hopt ¯K are given by (see
Fig.7) ¯ χK opt[n ] 0.36 0 φK opt(t) 1.38 0 ¯ χK opt(ej ) 0.81 0 -π -π 2 π2 π φK opt(jω) 1 0 2ωN -2ωN 4ωN -4ωN
Fig. 7 Discrete sampling function ¯χoptK[n] of the L2 optimal downsampler (top left) and its classical discrete time Fourier transform ¯χoptK(ejν) (bottom left) in Example2. The top left figure also shows
2a h(R2−1) (1 + R2) sinc t 2
− R sinc t+12 + sinc t−12 (dotted). The figure also shows the hold functionφoptK (t) of the L2optimal hold (top right) and its classical Fourier transformφoptK (jω) (bottom right) in Example2. Here a= 0.2 and ωn=πh
χK opt(e jν) = 2a h(R2− 1)(1 + R 2− 2R cos(ν))1 [−Mπ,Mπ](ν), φK opt(jω) = |Gv(jω)|2 i∈Z h1[−ωn,ωn](ω − 2Miωn) respectively.
Hereωn:= πh. Also, using the inverse lifted discrete Fourier transform and the classical inverse Fourier transform, we have that
¯χK opt[n]= 2a h(R2−1) (1+ R2) sinc n M − R sinc n+ 1 M +sinc n− 1 M φK opt(t) = i∈Z hψy(2)(t − ih) sinc i M .
See Fig.7. Since the optimal downsampler and hold are well-defined,Foptis an
order-1 hybrid interpolator. Using (12), it can be shown that the above optimal downsampler and hold are stable.
Finally, for a= 0.2, h = 1, M = 2, we find that
Gv− FoptSy2L2 = Gv2L2 − FoptSy2L2 = 2.5 − 0.6134 = 1.8866.
7 L
∞optimal downsampling
The L∞downsampling problem is the Problem1withK = L∞. It can be solved by obtaining a solution of Problem2 with induced 2-norm. Section7.1contains some basic results, which are important in solving the L∞optimal downsampling problem in Sect.7.2. The proofs of results in this section are given in AppendixC.
7.1 Preliminaries
This section describes some basic results which are useful in obtaining a solution of Problem2with induced 2-norm i.e. obtaining
`Fopt:= arg min
`F ˘Gv− `F ´Sy
with constraint rank `Fopt≤ r at each θ ∈ [−π, π]. Here the minimization is point-wise
inθ. At almost all θ, the operators ˘Gvand ´Syare compact operators (see Remark2).
Now, Problem2can be solved by considering a generic problem of obtaining
Fopt:= arg min F∈Fr
A − FB (34)
where A: H → H0and B : H → H1are compact operators, F : H1 → H0is a
operators mapping the spaceH1toH0and that have rank at most r . The spacesH,
H0andH1are assumed separable Hilbert spaces. Also, P : H → H represents the
orthogonal projection on(Ker B)⊥, and we define AN and Tγ as
AN := A(I − P), Tγ := I − γ−2ANA∗N ∀γ > AN,
where · represents the induced 2-norm of the operator throughout this section. NG for an operator G has the same meaning as in Sect.6.1(see (16)).
An algorithm for obtaining Fopt:= arg minF∈FrA − FB is presented at the end
of this section.
Similar to (17), the orthogonal projection P provides some clue about the lower bound ofA − FB for all bounded F. Indeed,
Lemma 6 infFA − FB ≥ AN.
The following is a standard but important step in L∞optimization [12]: Lemma 7 Ifγ > AN then
A − FB ≤ γ ⇔ T−12
γ (AP − FB) ≤ γ
It is also clear from Lemma7that ifγ > AN then the singular values of T− 1 2 γ A P
(or singular values of(AP)∗Tγ−1A P) provide a clue about optimal Foptat the given
γ .
We start with an SVD of A A∗ and other simplifications that are later useful in obtaining singular values of(AP)∗Tγ−1A P for a givenγ .
Lemma 8 Let an SVD of A A∗: H0→ H0be given by
A A∗x=
k∈NA
α2
k x, vk vk,
where{αk} is a non-increasing positive sequence. Define for a given fi, i ∈ NA P(see (18)) the subspace
Vi := span{vk|k ∈ NA& vk, fi = 0} (35)
and the orthogonal projection Pi onto the subspace Vi. Then fi ∈ Vi and an SVD (modulo ordering) of PiA A∗Pi : H0→ Vi exists and it has the form
PiA A∗Pix=
k∈NA, vk, fi =0
α2
k x, vk vk. (36)
In Lemma8, there is some degree of freedom in selection of the singular vectorsvk of A A∗. However, once the singular vectorsvk of A A∗are fixed then the spaceVi is fixed.
Corollary 1 The spacesVi (defined in Lemma8) are invariant under the operator
A A∗.
The above corollary follows from the fact that thevkare eigenvectors of A A∗. Without loss of generality an SVD of PiA A∗Pi can be written by rearranging (36) as PiA A∗Pix= n∈NPi A α2 in x, vin vin, x ∈ H0 (37)
where{α2in}n∈Nare in descending order,α2inare positive singular values of A A∗and
vinare singular vectors of A A∗such that vin, fi = 0.
In general, obtaining singular values of(AP)∗Tγ−1A P for a givenγ can be very
cumbersome but in some special cases it is relatively easy. One of the special case happens if the following assumption is satisfied.
Assumption A1: The subspaces{Vi}i∈NA Pdefined in (35) are mutually orthogonal. The above is a strong assumption. However in case of the downsampling this is true as we will see later.
Since fi ∈ Vi (see Lemma8), therefore if AssumptionA1is satisfied then fi ⊥
Vj, j = i. Hence, the orthogonal projection Pi on the space Vi satisfies Pifj = ¯δ[i − j] fj.
Under AssumptionA1, an SVD of the operator(AP)∗T− 1 2
γ A P can be obtained,
which eventually helps in obtaining Fopt, as follows:
Lemma 9 Suppose an SVD of A P exists and is given by (18). Assume that for each i ,
an SVD of PiA A∗Pi : H0→ Viexists and has the form given in (37). Ifγ > AN and
AssumptionA1is satisfied for the spaces{Vi}i∈NA P, then an SVD of(AP)∗Tγ−1A P
exists and it is given by (modulo ordering) (AP)∗Tγ−1A P x= i∈NA P η2 i(γ ) x, ei ei where η2 i(γ ) = ⎛ ⎝γ−2+ ! σ2 i k | fi, vi k |2 1− γ−2αi k2 "−1⎞ ⎠ −1 , (38)
and ei and fiare as defined in (18).
Even after obtaining an SVD of(AP)∗T−
1 2
γ A P, our original problem of obtaining
arg minF∈FrA − FB is still two-fold as both minF∈FrA − FB and Foptare still
Lemma 10 Suppose that an SVD of A P is given by (18). Given r ≤ rank AP, define
the setF that consists of all rank-r F of the form
F x= i∈C σi B+x, ei fi, x ∈ H1 (39)
whereC is any set of different r non-negative integers in the set NA Pand eiand fiare
as defined in (18). For convenience, we define F0∈ F as that F which has C = Dr(S) (see Definition6) whereS := {αi 0}i∈NA P. If Im B is closed and AssumptionA1is
satisfied for the spaces{Vi}i∈NA P, then F0minimizesA − FB over all F ∈ F and
A − F0B = max i/∈Dr(S),i∈N
{αi 0, AN} . (40)
In the following theorem, we show that F0defined in Lemma10minimizesA −
F B not only over the set F but also over Fr (see (34)) under certain conditions. Theorem 6 For a given r ≤ rank AP, let the set Dr(S) be as in Lemma10. Also
suppose that an SVD of A P has the form given in (18). If Im B is closed and
Assump-tionA1is satisfied for the spaces{Vi}i∈N, then F0(defined in Lemma10) minimizes
A − FB over all F ∈ Fr (see (34)).
The case of non-closed Im B in Lemma10and the above theorem is discussed in [20, Lemma 4.10.6, Theorem 4.10.7]. Note that Im B is closed if and only if B is a finite rank compact operator (see Remark3). Now, we can summarize all the steps for obtaining Fopt:= arg minF∈FrA − FB in a form of Algorithm2.
Algorithm 2 For Fopt:= arg minF∈Fr A − FB Require: A and B are compact operators and B is a finite rank operator.
Obtain the orthogonal projection P on(Ker B)⊥ Obtain SVD of A A∗and A P, and ei, fidefined in (18).
ObtainVidefined in Lemma8andS defined in Lemma10
if theVi’s satisfy AssumptionA1 then
Obtain B+
Fopt= F0(see Lemma10).
else Exit end if
CalculateA − FoptB using (40)
return Fopt,A − FoptB
7.2 L∞optimal downsampling solution
Problem2with induced 2-norm can be solved by application of Algorithm2. As in Sect.6, we defineSy := ShGy. The orthogonal projection operator ˘P on the space
(Ker ´Sy)⊥and an SVD of the operator ˘Gv ˘P at almost every θ ∈ [−π, π], required in Algorithm2, are given by (25) and Lemma5, respectively.
The next thing in Algorithm2is to find the invariant spacesVibased on an SVD of ˘Gv˘G∗vand left singular vectors˘qiof ˘Gv ˘P (see Lemma5). An SVD (modulo ordering) is given by [12]
˘Gv ˘G∗v ˘w = i∈Z
|Gv(jωi)|2 ˘w, ˘ei ˘ei,
where ˘ei are defined in Lemma4. An SVD (modulo ordering) of the operator ˘Gv ˘P :
L2[0, h) → L2[0, h) at almost all θ ∈ [−π, π] is given by (26). Hence, theVidefined in Lemma8in this case is given as
Vi := span{˘ei+Mk}k∈Qi, (41)
where Qi := {k ∈ N : ˘qi, ˘ei+Mk = 0} and ˘qi are defined in Lemma 5. Now, AssumptionA1is satisfied in the downsampling problem as shown in the following
lemma.
Lemma 11 Vi defined in (41) are mutually orthogonal.
By the above Lemma theViare (mutually) orthogonal (so AssumptionA1is satisfied),
therefore we can proceed further according to Algorithm2. Suppose ˘Pi represents the orthogonal projection onto the spaceVi at eachθ ∈ [−π, π]. Since, the ˘ek are eigenvectors of the operator ˘Gv˘G∗v, we have for every ˘w ∈ L2[0, h) that
˘Pi ˘Gv˘G∗v ˘Pi˘w = k∈Qi
|Gv(jωi+Mk)|2 ˘w, ˘ei+Mk ˘ei+Mk.
Now, we give a solution of Problem2with induced 2-norm using Theorem6. Theorem 7 Given Gv, Sy := ShGy ∈ L∞∩ L2. Let ˘P be given by (25), Qi be
as in (41) and ek be as in Lemma 4. Define|Gvmax,i| := maxl∈Qi{|Gv(jωi+Ml)|}, ˘GN := ˘Gv(I − ˘P), and S := {|Gvmax,i|}i∈M at almost everyθ ∈ [−π, π]. Given
r≤ M, define `Foptas `Fopt yh := k∈Dr(S) Γk yh,√1 M ek , (42)
where Dr(S) is defined in Definition6 and Γk is defined in (30), whenever θ ∈ A := {θ ∈ [−π, π] : ´Sy(ejθ) = 0} and τ ∈ [0, h). For θ ∈ [−π, π]\A, we can
take `Fopt = 0. If this interpolator Foptis well defined and stable, then it minimizes
Gv− FSyL∞ over all interpolatorsF with rank `F ≤ r at each θ ∈ [−π, π]. The
optimal norm is given by
Gv− FoptSyL∞ = ess sup
where ˘Gv− `Fopt´Sy = maxi/∈Dr(S),i∈M
|Gvmax,i|, ˘GN
.
The quantity ess supθ∈[−π,π] ˘GN is known as Parrott lower bound [15]. Note that
in (42) the setS is defined as a set of all singular values of ˘Gv˘P at a given θ, whereas in (29) the setS is defined as a set of maximal singular values of ˘Pi ˘Gv.
The interpolator `Foptin (42) can be written as a hybrid interpolator using the method
described in Remark7.
Example 3 Consider Example1but here we find the L∞optimal downsampling solu-tion instead of the L2optimal downsampling solution.
Since M = 2, we have M = {0, 1}. For i ∈ M, qi defined in Lemma5at a given
θ ∈ [−π, π] becomes
˘q0= √1
5(2˘e0+ ˘e21[−π,0](θ) + ˘e−21[0,π](θ)),
˘q1= √12(˘e−1+ ˘e1),
where˘ekare defined in Lemma4. Now, for almost allθ ∈ [−π, π] and i ∈ M, Gvmax,i defined in Theorem7becomes
|Gvmax,0| = 1, |Gvmax,1| = 0.95.
Clearly, D1(S) = {0} as defined in Theorem7. Using Theorem7, if r = 1 then an
L∞optimal interpolator in lifted frequency domain is given by
`Foptx= ⎧ ⎨ ⎩ 2 3e jω0τ +1 3e jω2τ x,1 2e jω0mh
, for almost all θ ∈ [−π, 0]
2 3e jω0τ +1 3e jω−2τ x,1 2e jω0mh
, for almost all θ ∈ [0, π]
where m∈ {0, 1}. Using (13) and the inverse lifted transform, we can write the discrete sampling function ¯χopt[n] and the hold function φopt(t) of the optimal downsampler
and hold as ¯χopt[n] = 12sinc(n2), φopt(t) = 2 3sinc(t) + 1 3cos 7πt 2 sinc t2 .
See Fig.8. Since the optimal downsampler and hold are well defined,Foptis an order-1
hybrid interpolator. Using (12), it can be shown that the above optimal downsampler and hold are stable.
The classical discrete time Fourier transform of ¯χopt[n] and classical Fourier
trans-formφopt(t) are given by
¯χopt(ejν) = 1[−π 2,π2](ν) φopt(jω) = 2 31[−ωn,ωn](ω) + 1 3 1[3ωn,4ωn](ω) + 1[−4ωn,−3ωn](ω) .
¯χopt[n] 0.5 0 φopt(t) 1 0 ¯χopt(ejν) 1 0 -π -π 2 π2 π φopt(jω) 2 3 1 3 0 2ωN -2ωN 4ωN -4ωN
Fig. 8 Discrete sampling function ¯χopt[n] of the L∞optimal downsampler (top left) and its classical
discrete time Fourier transform¯χopt(ejν) (bottom left) in Example3. The top left figure also shows12sinc(t2)
(dotted). The figure also shows the hold functionφopt(t) of the L∞optimal hold (top right) and its classical
Fourier transformφopt(jω) (bottom right) in Example3. Hereωn= πh
Since Gv(jω) has finite support, we can obtain an equivalent finite-dimensional matrix representation of the operator ˘GNfor norm calculation (see [26, §7.6]). This
yields ˘GN = 0.95 for all θ and
Gv− FoptSyL∞= ess sup
θ∈[−π,π]max(0.95, ˘GN) = 0.95.
Comparing with Example1, we can say that the L2and L∞downsampling problems may have different solutions.
Now, we take an example of non-bandlimited signal models.
Example 4 We consider again the system of Example2, but now we find the L∞
optimal downsampling solution instead of the L2optimal downsampling solution. For everyθ ∈ [−π, π] and i ∈ M, we have
|Gvmax,0| ≥ |Gvmax,1|.
This is because|Gv(jω)| is monotonically decreasing in [0, ∞) and monotonically increasing in (−∞, 0]. Hence, D1(S) = {0} for all θ ∈ [−π, π] in Theorem 7.
Therefore, the optimal downsampler and hold are as found in Example2.
To calculate the error, we use the following simple but useful technique. Since
Gv(jω) has infinite support, we can obtain an equivalent infinite-dimensional matrix representation of the operator ˘GN at eachθ ∈ [−π, π] (see [26, §7.6]). We denote
the infinite size matrix by A(θ) at a given θ. Now, we can easily make a sequence of matrices An(θ) such that A(θ) = limn→∞An(θ). Next,
lim
n→∞|A(θ) − An(θ)| ≤ limn→∞A(θ) − An(θ) = 0.
This impliesA(θ) = limn→∞An(θ). Since Gv(jω) is continuous at all ω ∈ R, every entry of A(θ) is continuous at all θ ∈ [−π, π] (considering the one-sided continuity at the end points). This implies that iterating over all θ ∈ [−π, π] for sufficiently large n, we can obtain ess supθ∈[−π,π] ˘GN. For a = 0.2, h = 1 and
M = 2, we obtain
Gv− FoptSyL∞ = max( ess sup
θ∈[−π,π]|Gvmax,1|, ess supθ∈[−π,π] ˘GN) = max(0.8467, 0.6227) = 0.8467.
8 Downsampling in the presence of noise
In this section, we will see the effect of colored noise on the L2and L∞downsampling problem. Noise analysis in this section is just an application of the theory we discussed in the previous sections; therefore, the proofs in this section are omitted. For a detailed discussion, see [20, §4.12]. The setup for noise analysis is shown in Fig.9. Here, ¯wc : Z → C is colored noise modeled by the h-time-invariant system ¯Gnwith input
signal ¯wn : Z → C which is white noise. The h-time-invariant system ¯Gn ∈ L∞is
defined as
¯wc= ¯Gn¯wn: ¯wc[k] =
i∈Z
¯gn[k − i] ¯wn[i] (44)
where ¯gnis the impulse response of the discrete system ¯Gn. We denote the z-transform
of ¯gnby ¯Gn(z). The SVD (modulo ordering) of the lifted operator ˘Gn(see (44)) at
almost everyθ ∈ [−π, π] is given by ˘Gn wn= M−1 k=0 ¯Gn(ejωkh ) wn, ek ek (45)
where the ekare defined in Lemma4. This SVD can be obtained by lifting ¯Gnand (9)
(see [20, lemma 4.11.1]).
For noise analysis, using the sampled-data system theory, we redraw Fig.9as Fig.
10[12]. Figure10is similar to Fig.1except forGv,Syand input signal. Now, the error
systemGein the presence of noise is defined as the mapping from
w ¯wn T to e: Ge:= GA− FGB whereGA:= Gv0 andGB:= Sy ¯Gn .
Our aim is to obtain a stable hybrid interpolatorF := H ¯Shof a given order r such thatGe2orGe∞is minimized.
