1

H ∞-Optimal Fractional Delay Filters
Masaaki Nagahara, Member, IEEE, Yutaka Yamamoto, Fellow, IEEE

Abstract—Fractional delay filters are digital filters to delay
discrete-time signals by a fraction of the sampling period. Since
the delay is fractional, the intersample behavior of the original
analog signal becomes crucial. In contrast to the conventional
designs based on the Shannon sampling theorem with the bandlimiting hypothesis, the present paper proposes a new approachbased on the modern sampled-data H ∞ optimization that aims at
restoring the intersample behavior beyond the Nyquist frequency.
By using the lifting transform or continuous-time blocking thedesign problem is equivalently reduced to a discrete-time H ∞optimization, which can be effectively solved by numerical computation softwares. Moreover, a closed-form solution is obtained
under an assumption on the original analog signals. Design
examples are given to illustrate the advantage of the proposed
method.
Index Terms—Fractional delay filters, interpolation, sampleddata systems, H ∞ optimization, linear matrix inequality.

I. I NTRODUCTION
Fractional delay filters are digital filters to delay discretetime signals by a fractional amount of the sampling period.
Such filters have wide applications in signal processing, including sampling rate conversion [1], [2], [3], nonuniform
sampling [4], [5], wavelet transform [6], [7], digital modelingof musical instruments [8], [9], to name a few. For more
applications, see survey papers [10], [11], [12].
Conventionally, fractional delay filters are designed basedon the Shannon sampling theorem [13], [14] for strictly
bandlimited analog signals. Based on this theory, the optimalfilter coefficients are obtained by sampling a delayed sincfunction. This ideal low-pass filter is however not realizable
because of its non-causality and instability, and hence manystudies have focused their attention on approximating the ideal
filter by, for example, windowed sinc functions [15], [16],maximally-flat FIR approximation [17], [18], [19], [20], [21],all-pass approximation [22], [23], and minmax (Chebyshev)
optimization [24].
In particular, H 2 (or weighted least-squares) design has
been prevalent in the literature [10], [25], [26], [21]. Thismethod minimizes the H 2 norm of the weighted difference
between the ideal low-pass filter and a filter to be designed, andis based on the projection theorem in Hilbert space. There are,however, two major drawbacks in this conventional approach.
Copyright (c) 2012 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
M. Nagahara is with Department of Applied Analysis and Complex Dynamical Systems, Graduate School of Informatics, Kyoto University, Kyoto,606-8501, JAPAN (e-mail: nagahara@ieee.org). Mailing address: Kyoto
University, Yoshida Honmachi, Sakyo-ku, Kyoto 606-8501, Japan.
Y. Yamamoto is with Department of Applied Analysis and Complex
Dynamical Systems, Graduate School of Informatics, Kyoto University, Kyoto, 606-8501, JAPAN (e-mail: yy@i.kyoto-u.ac.jp). Mailing address:
Kyoto University, Yoshida Honmachi, Sakyo-ku, Kyoto 606-8501, Japan.

One is that due to the averaging nature of the H 2 criterion,
the obtained frequency response can have a sharp peak at a
certain frequency, thereby yielding a poor performance at thatfrequency, while still maintaining small H 2 error in the overallfrequency response. The other is that H 2 criterion can yielda truncated frequency response as an optimal approximant of
the ideal low-pass filter, which yields a distortion due to theGibbs phenomenon in the time domain. Furthermore, such adesign is mostly executed in the discrete-time domain, whichyields poor intersample response.
In view of these problems we employ sampled-data H ∞ optimization, recently introduced for signal processing by [27]1 .
This is based on sampled-data control theory [29] which
accounts for the mixed nature of continuous- and discretetime thereby enabling optimization of the intersample signals
via discrete-time controllers (filters). This also allows for
optimization according to the H ∞ norm, namely minimizing
the maximum of the error frequency response. This worstcase design is clearly desirable in that it does not have thedrawback due to the averaging property of the H 2 criterion.
Due to the nature of the H ∞ norm, however, this optimization
problem has been difficult to solve, but one can now utilizea standardized method to solve this class of problems [30],[29]. Furthermore, the obtained filter shows greater robustness
against unknown disturbances due to the nature of the uniform
attenuation of the error frequency response; see [27] fordetails. Based on this H ∞ optimization method, we formulate
the design of fractional delay filters as a sampled-data H ∞optimization problem2 .
In order to optimize the intersample behavior, we must
deal with both continuous- and discrete-time signals, andhence the overall system is not time-invariant. The key tosolving this problem is lifting, which is introduced in theearly studies of modern sampled-data control theory [33], [34],[35], [36], [37]. Indeed, continuous-time lifting gives an exact,
not approximated, time-invariant discrete-time model for asampled-data system, albeit with infinite-dimensional input
and output spaces. Hence the problem of the mixed time sets
is circumvented without approximation.
Lifting can also be interpreted as a continuous-time blocking
or polyphase decomposition. As in multirate signal processing [38], lifting makes it possible to capture continuoustime signals and systems in the discrete-time domain without
approximation; see Section III-A for details. The remaining
system becomes a time-invariant discrete-time system, albeit
with infinite-dimensional input and output spaces. In view
1 The approach dates back to [28], though.
2 This method was first proposed in our conference articles [31], [32]. The
present paper reorganizes these works with new results on the state-space
formulation (Proposition 1, Appendix A). Simulation results in Section IV
are also new.

2of this infinite-dimensionality, we retain the term lifting to
avoid confusion. With such a representation, we show that ourdesign problem is reducible to a finite-dimensional discretetime H ∞ optimization without approximation. This type of
H ∞ optimization is easily solvable by standard softwares suchas MATLAB [39].
In some applications, digital filters with variable delay
responses (variable fractional delay filters [40], [10], [25],[26]) are desired. In this case, a filter should have a tunable
delay parameter, and hence a closed-form formula should bederived. In general, H ∞ optimal filters are difficult to solve
analytically. However, we provide a closed-form formula ofthe optimal filter with the delay variable as a parameter under
the assumption that the underlying frequency characteristic ofcontinuous-time input signals is governed by a low-pass filterof first order. While this assumption may appear somewhat
restrictive, it covers many typical cases and variations by somerobustness properties [31].
The paper is organized as follows. Section II defines fractional delay filters, and reviews a standard H 2 design method.
We then reformulate our design problem as a sampled-data
H ∞ optimization to overcome the difficulty due to the H 2
design. Section III gives a procedure to solve the sampled-data
H ∞ optimization problem based on the lifting transform. Section IV shows numerical examples to illustrate the superiorityof the proposed method.

Notation
Throughout this paper, we use the following notation. We
denote by L2 [0, ∞) and L2 [0, T ) the Lebesgue spaces consisting of all square integrable real functions on R+ := [0, ∞)and [0, T ), respectively. L2 [0, ∞) may be abbreviated as L2 .
By `2 we denote the set of all real-valued square summablesequences on Z+ := {0, 1, 2, . . .}. For a normed space X, wedenote by `2 (Z+ , X) the set of all sequences on Z+ taking
values in X with squared norms being summable. For normed
linear spaces X and Y , we denote by B(X, Y ) the set of
all bounded linear operators of X into Y . Rν and Rm×n
denote respectively the set of real vectors of size ν and realmatrices of size m×n. Finite-dimensional vectors and matrices
are denoted by bold letters, such as x or A, and infinitedimensional operators by calligraphic letters, such as B. The
transpose of a matrix A is denoted by A> . Symbols s andz are used for the variables of Laplace and Z transforms,
respectively. For a linear system F , its transfer function isdenoted by F̂ (z) (if F is discrete-time) or F̂ (s) (if F iscontinuous-time), and its impulse response√ by the lower-caseletter, f [n] or f (t). The imaginary unit −1 is denoted by j.
II. F RACTIONAL D ELAY F ILTERS
In this section, we review fractional delay filters withconventional design methods based on the Shannon sampling
theorem. Then, we reformulate the design problem as asampled-data H ∞ optimization problem.

v(t − D)v(t)(A) Delay
t

t

0

0 D
(C) Sampling(B) Samplingv(nT )

v(nT − D)(D) FDFt0t
0 D

Fig. 1. Fractional delay process: (A) a continuous-time signal v(t) (top left)
is delayed by D > 0. (B) the delayed signal v(t − D) is sampled at t = nT ,n = 0, 1, . . .. (C) the signal v(t) is sampled at t = nT , n = 0, 1, . . ..
(D) digital filtering (fractional delay filter, FDF) to produce (or estimate) the
sequence {v(nT − D)} from the sampled-data {v(nT )}.

A. Definition and standard design method
Consider a continuous-time signal v shown in Fig. 1 (topleft figure). Assume v(t) = 0 for t < 0 (i.e., it is a causal
signal). Delaying this signal by D > 0 gives the delayed
continuous-time signal v(t − D) shown in Fig. 1 (top-rightin Fig. 1). Then by sampling v(t − D) with sampling period
T , we obtain the discrete-time signal {v(nT − D)}n∈Z asshown in Fig. 1 (bottom-right in Fig. 1).
Next, let us consider the sampled signal {v(nT )}n∈Z of
the original analog signal v as shown in Fig. 1 (bottom-left in
Fig. 1). The objective of fractional delay filters is to reconstruct
or estimate the delayed sampled signal {v(nT − D)}n∈Z
directly from the sampled data {v(nT )}n∈Z when D is not
an integer multiple of T . We now define the ideal fractional
delay filter.
Definition 1: The ideal fractional delay filter Kid withdelay D > 0 is the mapping that produces {v(nT − D)}n∈Z
from {v(nT )}n∈Z , that is,
Kid : {v(nT )}n∈Z 7→ {v(nT − D)}n∈Z .
Assume for the moment that the original analog signal v isfully band-limited below the Nyquist frequency ΩN = π/T ,that is,v̂(jω) = 0, |ω| ≥ ΩN ,(1)
where v̂ is the Fourier transform of v. Then the impulse
response of the ideal fractional delay filter is obtained by [10]:sin π(n − D/T )= sinc(n − D/T ),π(n − D/T )(2)sin(πt)n = 0, ±1, ±2, . . . , sinc(t) :=.
πt
The frequency response of this ideal filter is given in thefrequency domain as
kid [n] =K̂id (ejωT ) = e−jωD ,ω ≤ ΩN .

(3)

Since the impulse response (2) does not vanish at n =−1, −2, . . . and is not absolutely summable, the ideal filter

3

(A)
w

W

v

(B)
ST

e−Ds
(C)
ST

ud

(D)
vd

K

ūd + ed−Fig. 2. Error system E for designing fractional delay filter K. (A)–(D)correspond to those in Fig. 1.

is noncausal and unstable, and hence the ideal filter is not
physically realizable. Conventional designs thus aim at approximating the impulse response (2) or the frequency response (3)by a causal and stable filter. We here review in particular theH 2 optimization, also known as weighted least squares [10].
Define the weighted approximation error by
E2 := (Kid − K)Wd

(4)

where Wd is a weighting function and K is a filter tobe designed, which is assumed to be FIR (finite impulse
response). The H 2 design aims at finding the FIR coefficientsof the transfer function K̂(z) of K that minimize the H 2 normof the weighted error system E2 :kE2 k22 = k(Kid − K)Wd k22
Z ΩN hi
2
1K̂id (ejωT ) − K̂(ejωT ) Ŵd (ejωT ) dω.
=
ΩN 0
(5)
As pointed out in the Introduction, this H 2 design has
some drawbacks. One is that the designed filter K may yielda large peak in the error frequency response Ê2 (ejωT ) due
to the averaging nature of the H 2 norm (5). If an input
signal has a frequency component at around such a peak ofÊ2 (ejωT ), the error will become very large. The second isthat the perfect band-limiting assumption (1) implies that theH 2 suboptimal filter is given as an approximant of the ideal
low-pass filter [41], which induces large errors in the time
domain [27]. Moreover, real analog signals always containfrequency components beyond the Nyquist frequency, andhence (1) never holds exactly for real signals.
B. Reformulation of design problem
To simultaneously solve the two problems pointed out
above, we introduce sampled-data H ∞ optimization [27]. Thismethod has advantages as mentioned in Section I. To adapt
sampled-data H ∞ optimization for the design of fractional
delay filters, we reformulate the design problem, instead of
mimicking the “ideal” filter given in (2) or (3).
Let us consider the error system shown in Fig. 2. W isa stable continuous-time system with strictly proper transfer
function Ŵ (s) that defines the frequency-domain characteristic of the original analog signal v. More precisely, we assume
that the analog signal v is in the following subspace of L2 :
W L2 := v ∈ L2 : v = W w, w ∈ L2 .

Note that the signal subspace W L2 is much wider than thatof band-limited L2 signals [42].
The upper path of the diagram in Fig. 2 is the ideal
process of the fractional delay filter (the process (A) → (B)in Fig. 1); that is, the continuous-time signal v is delayedby the continuous-time delay denoted by e−Ds (we use thenotation e−Ds , the transfer function of the D-delay system,as the system itself), and then sampled by the ideal sampler
denoted by ST with period T > 0 to become an `2 signal3
ud := ST e−Ds v, or
ud [n] := ST e−Ds v [n] = v(nT − D), n ∈ Z+ .
On the other hand, the lower path represents the real process((C) → (D) in Fig. 1); that is, the continuous-time signal v is
directly sampled with the same period T to produce a discretetime signal vd ∈ `2 defined by
vd [n] := (ST v) [n] = v(nT ),n ∈ Z+ .

This signal is then filtered by a digital filter K to be designed,
and we obtain an estimation signal ūd = KST v ∈ `2 .
Put ed := ud − ūd (the difference between the ideal output
ud and the estimation ūd ), and let E denote the error systemfrom w ∈ L2 to ed ∈ `2 (see Fig. 2). Symbolically, E isrepresented by (cf. (4))

E = ST e−Ds − KST W.
(6)Then our problem is to find a digital filter K that minimizes
the H ∞ norm of the error system E.
Problem 1: Given a stable, strictly proper W (s), a delay
time D > 0, and a sampling period T > 0, find the digital
filter K that minimizes (cf. (5))

kEk∞ = ST e−Ds − KST W ∞=supk ST e−Ds − KST W wk`2 .
w∈L2 , kwk2 =1

Note that W , or its transfer function Ŵ (s), can be interpreted as a frequency-domain weighting function for the
optimization. This is comparable to Ŵd (z) in the discrete-time
H 2 design minimizing (5). The point to use continuous-timeŴ (s) is that one can model the frequency characteristic ofsignals beyond the Nyquist frequency. Also, the advantage ofusing the sampled-data setup here is that we can minimize thenorm of the overall transfer operator from continuous-time w
to the error ed . In the next section, we will show a procedure
to solve Problem 1 based on sampled-data control theory.
III. H ∞ D ESIGN OF F RACTIONAL D ELAY F ILTERS
The error system E in Fig. 2 contains both continuousand discrete-time signals, and hence the system is not timeinvariant; in fact, it is T -periodic [29]. In this section, we
introduce the continuous-time lifting technique [37], [29] toderive a norm-preserving transformation from E to a timeinvariant finite-dimensional discrete-time system. After this,
3 If Ŵ (s) is stable and strictly proper, the discrete-time signal ud =ST e−Ds v belongs to `2 . Otherwise, ST is not a bounded operator on L2 ;see [29, Section 9.3].

4

w


LW L−1v

ST L−1
Fig. 3.

vd

K(z)ūd + ed−Fig. 4.

"!

!

!"

ud

ST e−Ds L−1

Factorization of E.

Lifted error system E.

one can use a standard discrete-time H ∞ optimization implemented on a computer software such as MATLAB to obtainan optimal filter. We also give a closed-form solution of the
optimization under an assumption on W (s).

Let {A, B, C} be a minimal realization [43] of Ŵ (s):dx(t)= Ax(t) + Bw(t), v(t) = Cx(t), t ∈ R+ ,dt

(7)

where x(t) ∈ Rν is the state variable (ν is a positive integer).
We assume A ∈ Rν×ν , B ∈ Rν×1 , C ∈ R1×ν , andx(0) = 0. Let D = mT + d where m ∈ Z+ and d is a real
number such that 0 ≤ d < T . First, we introduce the lifting
operator L [37], [29] that transforms a continuous-time signal
in L2 [0, ∞) to an `2 sequence of functions in L2 [0, T ). Apply
L to the continuous-time signals w and v, and put w
e := Lw,
ve := Lv. By this, the error system in Fig. 2, is transformed
into a time-invariant discrete-time system E shown in Fig. 3.
Since the operator L gives an isometry between L2 [0, ∞) and`˜2 := `2 (Z+ , L2 [0, T )), we havesupw∈e `˜2 , kwk
e 2˜=1
`

kE wk
e `2

(8)

The following proposition is fundamental to the sampleddata optimization in (8).
Proposition 1: A state-space realization of the lifted error
system E is given byBξ[n + 1] = Ad ξ[n] +
w[n],
e
0
(9)
ed [n] = C 1 ξ[n] − ūd [n], vd [n] = C 2 ξ[n],ūd [n] = (k ∗ vd )[n],where ∗ stands for convolution, and the pertinent operators
Ad , B, C 1 and C 2 are given as follows: First, B is a linear(infinite-dimensional) operator defined by
2

ν+1

B : L [0, T ] → R
,


B w
e
w
e 7→ B w
e= 1
B2 w
e
 Z T

A(T −τ )



e
B w(τ
e )dτ


0

=
Z T −dA(T −d−τ )
Ce
B w(τ
e )dτ
0

C 1 := [0, 0, C m ] ∈ R1×(ν+1+m) ,C 2 := [C, 0, 0] ∈ R1×(ν+1+m) ,A. Lifted model of sampled-data error system

kEk∞ = kEk∞ :=The matrices Ad , C 1 , and C 2 in (9) are defined by


eAT
0
0
Ad :=  CeA(T −d)00  ∈ R(ν+1+m)×(ν+1+m) ,0
B m Am(10)

where Am , B m , and C m are state-space realization matricesof the discrete-time delay z −m .
Proof: See Appendix A.
The state-space representation (9) then gives the transfer
function of the lifted system E as
Ê(z) = Ĝ1 (z) − K̂(z)Ĝ2 (z),whereĜi (z) := C i (zI − Ad )
Put

−1



B
0



,

(11)

i = 1, 2.

Ê0 (z) := C 1 − K̂(z)C 2 (zI − Ad )−1 .

(12)Note that E0 is a finite-dimensional discrete-time system. Thenthe lifted system E(z) in (11) can be factorized (see Fig. 4)asBÊ(z) = Ê0 (z)
.
(13)
0

B. Norm-equivalent finite-dimensional system
The lifted system E given in (9), or its transfer function Ê in(11), involves an infinite-dimensional operator B : L2 [0, T ) →Rν+1 . Introducing the dual operator [44] B ∗ : Rν+1 →L2 [0, T ) of B, and composing this with B, we can obtaina norm-equivalent finite dimensional system of the infinitedimensional system E.
The dual operator B ∗ of B is given by [44]>B ∗ = B1∗ B2∗ , B1∗ (θ) := B > eA (T −θ) ,>B2∗ (θ) := 1[0,T −d) (θ)B > eA (T −d−θ) C > ,θ ∈ [0, T ),where 1[0,T −d) is the characteristic function of the interval[0, T − d), that is,(1,θ ∈ [0, T − d),1[0,T −d) (θ) :=0,otherwise.
Then we have the following lemma:

5

!"!

Fig. 5.

!

!

"!

Discrete-time system Ed .

Lemma 1: The operator BB ∗ is a positive semi-definite
matrix given by
B B ∗ B1 B2∗
BB ∗ = 1 1∗B2 B1 B2 B2∗(14)
M (T )
eAd M (T − d)C >=>CM (T − d)eA d CM (T − d)C >where M (·) is defined by
Z t>
M (t) :=eAθ BB > eA θ dθ ∈ Rν×ν ,0

t ≥ 0.

Proof: First, the equality in (8) and E ∈ B(L2 , `2 ) give
kEk∞ = kEk∞ < ∞. Using the factorization (13), we have

 2
B
2
2
kEk∞ = kEk∞ = E00 ∞ ∗
B
B
= E0
E0∗
0
0
∞


>
Bd
Bd
2
= kEd k∞ .
= E0
E0∗
0
0
∞

Thus the sampled-data H ∞ optimization (Problem 1) isequivalently transformed to discrete-time H ∞ optimization.
A MATLAB code for the H ∞ -optimal fractional delay filter
is available on the web at [46]. Moreover, if we assume thatthe filter K(z) is an FIR filter, the design is reduced to aconvex optimization with a linear matrix inequality. See [32],[47] for details.

Proof: We first prove B1 B2∗ = eAd M (T − d)C > .
For every u ∈ R, we haveZ T
∗B1 B2 u =eA(T −θ) B B2∗ (θ)u dθ
0
Z T


>
=
eA(T −θ) B 1[0,T −d) (θ)B > eA (T −d−θ) C > u dθ0
Z T −d>Ad=eeA(T −d−θ) BB > eA (T −d−θ) dθ C > u

C. Closed-form solution under a first-order assumption
Assume that the weighting function Ŵ (s) is a first-order
low-pass filter with cutoff frequency ωc > 0:
ωc
.
(15)
Ŵ (s) =s + ωcUnder this assumption, a closed-form solution for the optimalfilter is obtained [31], [32]:Theorem 2: Assume that Ŵ (s) is given by (15). Then the
optimal filter K̂(z) is given by

Similarly, we can prove the equalities B1 B1∗ = M (T ) andB2 B2∗ = CM (T − d)C > .
Remark 1: The matrix M (t) can be computed via thematrix exponential formula [45]:

where

0

= eAd M (T − d)C > u.

M (t) = F >22 (t)F 12 (t),



 
F 11 (t) F 12 (t)−A BB >:= expt .
0
F 22 (t)0A>K̂(z) = a0 (d)z −m + a1 (d)z −m−1 ,a0 (d) :=(16)sinh (ωc (T − d)), a1 (d) := e−ωc T eωc d − a0 (d) .
sinh(ωc T )Moreover, the optimal value of kEk∞ is given bysωc sinh(ωc d) sinh(ωc (T − d))kEk∞ =.
sinh(ωc T )(17)

By this formula, we can easily compute the matrices M (T )and M (T − d) in (14) without performing a numericalintegration.
From Lemma 1, BB ∗ is a positive semi-definite matrix andhence there exists a matrix B d such that BB ∗ = B d B >d . With
matrix B d and discrete-time system E0 given in (12), define
a finite-dimensional discrete-time system by
Bd
Ed := E0.
0

Since the optimal filter K̂(z) in (16) is a function of thefractional delay d and the integer delay m, the filter can beused as a variable fractional delay filter [10].
Remark 2: Fix d > 0 and m ∈ Z+ arbitrarily. By definition,
we have T − d < T . It follows that as ωc → ∞, we havea0 (d) → 0, a1 (d) → 0, and kEk∞ → ∞. This meansthat if the original analog signals contain higher frequency
components (far beyond the Nyquist frequency), the worstcase input signal becomes more severe, and the H ∞ -optimalfilter becomes closer to 0.

See Fig. 5 for the block diagram of Ed . Then the discretetime system Ed is equivalent to the original sampled-data error
system E in Fig. 2 with respect to their H ∞ norm as describedin the following theorem:
Theorem 1: Assume that the sampled-data error system E
gives an operator belonging to B(L2 , `2 ), the set of all boundedlinear operators of L2 into `2 . Then the discrete-time system
Ed belongs to B(`2 , `2 ) and equivalent to E with respect totheir H ∞ norm, that is, kEk∞ = kEd k∞ .

IV. D ESIGN E XAMPLESWe here present design examples of fractional delay filters.
The design parameters are as follows: the sampling period
T = 1 (sec), the delay D = 5.5 (sec), that is, m = 5 and d =0.5. The frequency-domain characteristic of analog signals to
be sampled is modeled byωcŴ (s) =, ωc = 0.1.
s + ωc6Frequency Reponse Gain

Frequency Response−23.5

Magnitude (dB)

0
−10

−24

−20
−30
−40
−3
10

−2

10

−1

10
Frequency (rad/s)

0

10

1

10

Frequency Response
0

−25

−25.5

−20
Phase (deg)Magnitude (dB)

−24.5

−40
−26
−60
−80
−3

10

−2

10

−1

10
Frequency (rad/s)

0

10

1

10

Fig. 6. Bode plot of Ŵ (s) (solid) and its impulse-invariant discretizationŴd (z) (dash). The vertical line shows the Nyquist frequency.

−26.5
−3
10

30
amplitude

Magnitude (dB)

1

10

40

0
−10
−20
−30

20
10
0
−10
−20

−2

10

−1

10
Frequency (rad/s)

0

10

1

10

40

−200

30

−400

20

amplitude

0

−600

0

−1000

−10
−20
−2

−1

10
Frequency (rad/s)

0.5

1

1.5

2
time (sec)

2.5

3

3.5

4

3

3.5

4

10

−800

10

0

Time Response (H−2 design)Frequency Response

Phase (deg)010Time Response (H−inf design)Frequency Response

−1200
−3
10

−1

10
Frequency (rad/s)

Fig. 8. Frequency response gain of error system E in Fig. 2: sampled-data
H ∞ design (solid), and discrete-time H 2 design (dash).

10

−40
−3
10

−2

10

0

10

1

10

0

0.5

1

1.5

2
time (sec)

2.5

Fig. 7. Bode plot of filters: sampled-data H ∞ design (solid), and discretetime H 2 design (dash).

Fig. 9. Time response: sampled-data H ∞ design (above) and discrete-time
H 2 design (below) with sampling frequency 1000 (Hz).

Note that Ŵ (s) has the cutoff frequency ωc = 0.1 (rad/sec) ≈0.016 (Hz), which is below the Nyquist frequency π (rad/sec)= 0.5 (Hz).
We compare the sampled-data H ∞ optimal filter obtainedby Theorem 2 with conventional FIR filters designed bydiscrete-time H 2 optimization [10], which minimizes the costfunction (5). The weighting function Ŵd (z) in (4) or (5) is
chosen as the impulse-invariant discretization [48] of Ŵ (s).
Fig. 6 shows the Bode plots of Ŵ (s) and Ŵd (z).
The transfer function of the proposed filter is given byK̂(z) = z −5 0.4994 + 0.4994z −1 .

E shown in Fig. 2. This is because the conventional designs
cannot take into account the frequency response of the source
analog signals while the present method does.
To see the difference between the present filter and theconventional one, we show the time response against a piecewise regular signal produced by the MakeSignal functionof WaveLab [49] in Fig. 9. The present method is superior
to the conventional one that shows much ringing at edgesof the wave. To see the difference more finely, we show thereconstruction error in Fig. 10. The H 2 -optimal filter has muchlarger errors around edges of the signal than the proposed H ∞ optimal one. In fact, the L2 norm of the error is 1.34 × 10−2
for H ∞ design and 2.07×10−2 for H 2 design. This illustrates
the effectiveness of our method.

Fig. 7 shows the Bode plots of the designed filters. As
illustrated in Fig. 7, the H 2 optimal filter is closer to the ideal
filter (3) as expected, so that it appears better in the context
of the conventional design methodology.
However, the H 2 optimal filter exhibits much larger errors
in the high-frequency domain as shown in Fig. 8 that shows
the frequency response gain of the sampled-data error system

V. C ONCLUSION
We have presented a new method of designing fractional
delay filters via sampled-data H ∞ optimization. An advantage
here is that an optimal analog performance can be attained.

7

where B1 is defined in (10). On the other hand, from (7),
we have v(t) = Cx(t) for t ∈ R+ . Putting t1 := nT andt2 := nT + θ for n ∈ Z+ and θ ∈ [0, T ), we have

Error (H−inf design)4

amplitude3

2v(nT + θ) = Cx(nT + θ)

1

0

= Ce
0

0.5

1

1.5

2
time (sec)

2.5

3

3.5

Aθ

x(nT ) +Z θ

CeA(θ−τ ) Bw(nT + τ )dτ.

0

4

(19)

Error (H−2 design)
4

By this, we haveamplitude3

vd [n] = v(nT ) = Cx1 [n].

(20)

2

Next, from (19), we have
1

0

0

0.5

1

1.5

2
time (sec)

2.5

3

3.5

yd [n] = v(nT − d) = v(nT − T + T − d)4

Fig. 10. Absolute value of reconstruction error: sampled-data H ∞ design(above) and discrete-time H 2 design (below).

= CeA(T −d) x1 [n − 1]
Z T −d+CeA(θ−τ ) B w[ne − 1](τ )dτ.
0

Put x2 [n] := yd [n]. Then we have
The optimal design problem can be equivalently transformed
to discrete-time H ∞ optimization, which is easily executedby standard numerical optimization toolboxes. A closed-form
solution is given when the frequency distribution of the input
analog signal is modeled as a first-order low-pass filter. Design
examples show that the H ∞ -optimal filter exhibits a much
more satisfactory performance than the conventional H 2 optimal filter.

x2 [n + 1] = CeA(T −d) x1 [n] + B2 w[n],eyd [n] = x2 [n],where B2 is defined in (10). By the relation

ud [n] = yd [n − m] = z −m yd [n],and the state-space matrices Am , B m , and C m for m-stepdelay z −m , we havex3 [n + 1] = Am x3 [n] + B m yd [n],ud [n] = C m x3 [n].

ACKNOWLEDGEMENT
This research is supported in part by the JSPS Grant-in-Aid
for Scientific Research (B) No. 24360163, (C) No. 24560543,
and Grant-in-Aid for Exploratory Research No. 22656095.

(21)

(22)

Combining (18), (20), (21), and (22) all together gives the> >state-space representation (9) with ξ > := [x>1 , x2 , x3 ] .
R EFERENCESA PPENDIX A
P ROOF OF P ROPOSITION 1
From the relation (6), the lifted system E is described as(see also Fig. 3)
E = ST e−Ds W L−1 − KST W L−1= ST e−mT s e−ds W L−1 − KST W L−1= z −m ST e−ds W L−1 − KST W L−1= z −m yd − Kvd ,where yd := ST e−ds W L−1 w
e and vd := ST W L−1 w.
e Fromthe state-space representation of W in (7), for any t1 and t2
such that 0 ≤ t1 ≤ t2 < ∞, we have
Z t2
A(t2 −t1 )x(t2 ) = ex(t1 ) +eA(t2 −t) Bw(t)dt.
t1

Putting t1 := nT and t2 := (n + 1)T for n ∈ Z+ gives
Z T
x(nT + T ) = eAT x(nT ) +eA(T −τ ) Bw(nT + τ )dτ.
0

Define x1 [n] := x(nT ) and w[n]e := (Lw)[n]. Then we havex1 [n + 1] = eAT x1 [n] + B1 w[n],e(18)[1] T. A. Ramstad, “Digital methods for conversion between arbitrary
sampling frequencies,” IEEE Trans. Acoust., Speech, Signal Processing,vol. 32, no. 3, pp. 577–591, 1984.
[2] J. O. Smith and P. Gossett, “A flexible sampling-rate conversion
method,” in IEEE ICASSP’84, 1984, pp. 19.4.1–19.4.4.
[3] S. Park, G. Hillman, and R. Robles, “A novel structure for real-time
digital sample-rate converters with finite precision error analysis,” inIEEE ICASSP’91, 1991, pp. 3613–3616.
[4] H. Johansson and P. Löwenborg, “Reconstruction of nonuniformly
sampled bandlimited signals by means of digital fractional delay filters,”
IEEE Trans. Signal Processing, vol. 50, no. 11, pp. 2757 – 2767, Nov.
2002.
[5] R. Prendergast, B. Levy, and P. Hurst, “Reconstruction of band-limited
periodic nonuniformly sampled signals through multirate filter banks,”
IEEE Trans. Circuits Syst. I, vol. 51, no. 8, pp. 1612–1622, Aug. 2004.
[6] R. Yu, “Characterization and sampled-data design of dual-tree filterbanks for Hilbert transform pairs of wavelet bases,” IEEE Trans. Signal
Processing, vol. 55, no. 6, pp. 2458 –2471, Jun. 2007.
[7] B. Dumitrescu, “SDP approximation of a fractional delay and the designof dual-tree complex wavelet transform,” IEEE Trans. Signal Processing,vol. 56, no. 9, pp. 4255 –4262, Sep. 2008.
[8] V. V. H.-M. Lehtonen and T. I. Laakso, “Musical signal analysis using
fractional-delay inverse comb filters,” in Proc. of the 10th Int. Conf. on
Digital Audio Effects, Sep. 2007, pp. 261–268.
[9] V. Välimäki, J. Pakarinen, C. Erkut, and M. Karjalainen, “Discrete-timemodelling of musical instruments,” Rep. Prog. Phys, vol. 69, no. 1, pp.
1–78, 2006.
[10] T. I. Laakso, V. Välimäki, M. Karjalainen, and U. K. Laine, “Splitting
the unit delay,” IEEE Signal Processing Mag., vol. 13, pp. 30–60, 1996.
[11] V. Välimäki and T. I. Laakso, “Principles of fractional delay filters,” inIEEE ICASSP’00, 2000, pp. 3870–3873.

8[12] V. Välimäki and T. I. Laakso, “Fractional delay filters — design
and applications,” in Nonuniform Sampling.
New York: Kluwer
Academic/Plenum Publishers, 2001, pp. 835–885.
[13] C. E. Shannon, “Communication in the presence of noise,” Proc. IRE,
vol. 37, no. 1, pp. 10–21, 1949.
[14] M. Unser, “Sampling — 50 years after Shannon,” Proc. IEEE, vol. 88,
no. 4, 2000.
[15] G. D. Cain, A. Yardim, and P. Henry, “Offset windowing for FIR
fractional-sample delay,” in IEEE ICASSP’95, 1995, pp. 1276–1279.
[16] J. Selva, “An efficient structure for the design of variable fractional delay
filters based on the windowing method,” IEEE Trans. Signal Processing,vol. 56, no. 8, pp. 3770 –3775, Aug. 2008.
[17] E. Hermanowicz, “Explicit formulas for weighting coefficients of maximally flat tunable FIR delayers,” Electronics Letters, vol. 28, pp. 1936–
1937, 1992.
[18] S.-C. Pei and P.-H. Wang, “Closed-form design of maximally flat FIR
Hilbert transformers, differentiators, and fractional delayers by power
series expansion,” IEEE Trans. Circuits Syst. I, vol. 48, pp. 389–398,
2001.
[19] S. Samadi, O. Ahmad, and M. N. S. Swamy, “Results on maximally flat
fractional-delay systems,” IEEE Trans. Circuits Syst. I, vol. 51, no. 11,
pp. 2271–2286, 2004.
[20] H. Hachabiboglu, B. Gunel, and A. Kondoz, “Analysis of root displacement interpolation method for tunable allpass fractional-delay filters,”
IEEE Trans. Signal Processing, vol. 55, no. 10, pp. 4896 –4906, Oct.
2007.
[21] J.-J. Shyu and S.-C. Pei, “A generalized approach to the design ofvariable fractional-delay FIR digital filters,” Signal Processing, vol. 88,
no. 6, pp. 1428–1435, Jun. 2008.
[22] Z. Jing, “A new method for digital all-pass filter design,” IEEE Trans.
Acoust., Speech, Signal Processing, vol. 35, pp. 1557–1564, 1987.
[23] S.-C. Pei and P.-H. Wang, “Closed-form design of all-pass fractional
delay filters,” IEEE Signal Processing Lett., vol. 11, pp. 788–791, 2004.
[24] W. Putnam and J. Smith, “Design of fractional delay filters using
convex optimization,” in Applications of Signal Processing to Audio andAcoustics, 1997 IEEE ASSP Workshop on, Oct 1997.
[25] A. Tarczynski, G. D. Cain, E. Hermanowicz, and M. Rojewski, “WLS
design of variable frequency response FIR filters,” in Proc. IEEE Int.
Symp. Circuits Syst., 1997, pp. 2244–2247.
[26] T.-B. Deng and Y. Nakagawa, “SVD-based design and new structures forvariable fractional-delay digital filters,” IEEE Trans. Signal Processing,vol. 52, pp. 2513–2527, 1991.
[27] Y. Yamamoto, M. Nagahara, and P. P. Khargonekar, “Signal reconstruction via H ∞ sampled-data control theory — Beyond the Shannon
paradigm,” IEEE Trans. Signal Processing, vol. 60, no. 2, pp. 613–625,
2012.
[28] P. P. Khargonekar and Y. Yamamoto, “Delayed signal reconstruction
using sampled-data control,” in Proc. 35th IEEE CDC, 1996, pp. 1259–
1263.
[29] T. Chen and B. A. Francis, Optimal Sampled-data Control Systems.
Springer, 1995.
[30] J. C. Doyle, K. Glover, P. P. Khargonekar, and B. A. Francis, “Statespace solutions to standard H∞ and H2 control problems,” IEEE Trans.
Automat. Contr., vol. 34, pp. 831–847, 1989.
[31] M. Nagahara and Y. Yamamoto, “Optimal design of fractional delay
filters,” in Proc. of 35th Conf. on Decision and Control, 2003, pp. 6539–
6544.
[32] M. Nagahara and Y. Yamamoto, “Optimal design of fractional delay FIR
filters without band-limiting assumption,” in IEEE ICASSP’05, 2005, pp.
221–224.
[33] Y. Yamamoto, “New approach to sampled-data systems: a function spacemethod,” in Proc. 29th Conf. on Decision and Control, 1990, pp. 1881–
1887.
[34] B. Bamieh, J. B. Pearson, B. A. Francis, and A. Tannenbaum, “A lifting
technique for linear periodic systems with applications to sampled-data
control systems,” Syst. Control Lett., vol. 17, pp. 79–88, 1991.
[35] H. T. Toivonen, “Sampled-data control of continuous-time systems withan H∞ optimality criterion,” Automatica, vol. 28, pp. 45–54, 1992.
[36] B. Bamieh and J. B. Pearson, “A general framework for linear periodic
systems with applications to H∞ sampled-data control,” IEEE Trans.
Automat. Contr., vol. 37, pp. 418–435, 1992.
[37] Y. Yamamoto, “A function space approach to sampled-data control
systems and tracking problems,” IEEE Trans. Automat. Contr., vol. 39,
pp. 703–712, 1994.
[38] P. P. Vaidyanathan, Multirate Systems and Filter Banks. Prentice Hall,
1993.

[39] G. Balas, R. Chiang, A. Packard, and M. Safonov, Robust Control
Toolbox, Version 3. The Math Works, 2005.
[40] C. W. Farrow, “A continuously variable digital delay element,” in Proc.
IEEE Int. Symp. Circuits Syst., vol. 3, 1988, pp. 2641–2645.
[41] N. J. Fliege, Multirate Digital Signal Processing. New York: John
Wiley, 1994.
[42] M. Nagahara, M. Ogura, and Y. Yamamoto, “H ∞ design of periodicallynonuniform interpolation and decimation for non-band-limited signals,”
SICE Journal of Control, Measurement, and System Integration, vol. 4,
no. 5, pp. 341–348, 2011.
[43] W. J. Rugh, Linear Systems Theory. Prentice Hall, 1996.
[44] Y. Yamamoto, “On the state space and frequency domain characterization of H ∞ -norm of sampled-data systems,” Syst. Control Lett., vol. 21,
pp. 163–172, 1993.
[45] C. F. V. Loan, “Computing integrals involving the matrix exponential,”
IEEE Trans. Automat. Contr., vol. 23, pp. 395–404, 1994.
[46] http://www-ics.acs.i.kyoto-u.ac.jp/mat/fdf/.
[47] M. Nagahara, “Min-max design of FIR digital filters by semidefinite
programming,” in Applications of Digital Signal Processing. InTech,
Nov. 2011.
[48] L. Jackson, “A correction to impulse invariance,” IEEE Signal Processing Lett., vol. 7, no. 10, pp. 273–275, Oct. 2000.
[49] D. Donoho, A. Maleki, and M. Shaharam, “Wavelab 850,” 2006.

Masaaki Nagahara (S’00–M’03) received the
Bachelor’s degree in engineering from Kobe University in 1998, the Master’s degree and the Doctoral
degree in informatics from Kyoto University in 2000
and 2003. He is currently a Lecturer at Graduate
School of Informatics, Kyoto University. His research interests include digital signal processing anddigital control systems. He received Young Authors
Award in 1999 and Best Paper Award in 2012 fromthe Society of Instrument and Control Engineers(SICE), and Transition to Practice Award from IEEE
Control Systems Society in 2012. He is a member of ISCIE, SICE, and IEICE.

Yutaka Yamamoto (M’83–SM’93–F’98) received
the Ph.D. degree in mathematics from the University
of Florida, Gainesville, in 1978, under the guidanceof Professor Rudolf E. Kalman. He is currently aProfessor in the Department of Applied Analysis andComplex Dynamical Systems, Graduate School ofInformatics, Kyoto University, Kyoto, Japan.
His current research interests include the theoryof sampled data control systems, its application todigital signal processing, realization and robust control of distributed parameter systems and repetitive
control.
Dr. Yamamoto was the recipient of G. S. Axelby Outstanding Paper Awardof the IEEE Control Systems Society (CSS) in 1996, the Commendation for
Science and Technology by the Minister of Education, Culture, Sports, Scienceand Technology, Prize for Science and Technology, in Research Category in
2007, Distinguished Member Award of the Control Systems Society of theIEEE in 2009, Transition to Practice Award of the IEEE CSS in 2012, andseveral other awards from the Society of Instrument and Control Engineers(SICE) and the Institute of Systems, Control and Information Engineers(ISCIE). He has served as Senior Editor of the IEEE Trans. Automatic Control(TAC) for 2010–2011, and as an associate editor for several journals including
TAC, Automatica, Mathematics of Control, Signals and Systems, Systems andControl Letters. He was Chair of the Steering Committee of MTNS for 2006–
2008, and the General Chair of the MTNS 2006 in Kyoto.
He is currently President of the Control Systems Society (CSS) of the IEEE,and was a vice president for 2005–2008 of CSS, and President of ISCIE of
Japan for 2008–2009. He is a fellow of the IEEE and SICE.