CAUSAL SPLINE INTERPOLATION BY H ∞ OPTIMIZATIONM. Nagahara and Y. Yamamoto
Graduate School of Informatics, Kyoto University, Kyoto, 606-8501, Japan
ABSTRACT
Spline interpolation systems generally contain non-causal ſlters, and hence such systems are difſcult to use for real-timeprocessing. Our objective is to design a causal system whichapproximates spline interpolation. This is formulated as aproblem of designing a stable inverse of a system with unstable zeros. For this purpose, we adopt H ∞ optimization.
We give a closed form solution to the H ∞ optimization in thecase of the cubic spline. For higher order splines, the optimalſlter can be effectively solved by a numerical computation.
We also show that the optimal FIR (Finite Impulse Response)
ſlter can be designed by an LMI (Linear Matrix Inequality),
which can also be effectively solved numerically. A design
example is presented to illustrate the result.
Index Terms— spline functions, interpolation, H ∞ optimization
1. INTRODUCTION
Signal interpolation has many applications; it is used for curve
ſtting, signal reconstruction, and sampling rate conversion including resolution conversion of digital images. Many methods have been proposed for signal interpolation such as polynomial splines [4, 5] and exponential splines [6]. Polynomial
splines are, in particular, widely used in image processing.
In polynomial (or exponential) spline interpolation, it isassumed that the original signal is a piece-wise polynomial(or exponential) function. Then, intersample values are computed via the Fourier coefſcients relative to the spline bases.
However, the coefſcients are computed by using the futuresamples as well as the present and the past ones, and hence
the interpolation system becomes non-causal. The same nature applies to the signal reconstruction by Shannon sampling
theorem [3]. In the case of image processing, non-causality
is not a restriction, and hence spline interpolation is widelyused in that ſeld. It is however difſcult to use the splines forreal-time processing such as instrumentation or audio/speechprocessing.
We therefore propose to design a causal system which approximates spline interpolation. This is formulated as a problem of designing a stable inverse of a system with unstable zeros. For this purpose, we adopt H ∞ optimization. By this, wecan obtain the H ∞ -optimal stable inverse, and hence causal

1­4244­0728­1/07/$20.00 ©2007 IEEEspline interpolation is obtained. Moreover, by assuming thatthe ſlter to be designed is an FIR (Finite Impulse Response)
ſlter, the optimization is reducible to an LMI (Linear MatrixInequality), which can be effectively solved by, for example,
standard MATLAB routines. In this article, we discuss polynomial spline interpolation. Exponential spline interpolation
can be discussed in the same way.
2. SPLINE INTERPOLATIONConsider x ∈ V N , where V N is the space of polynomial
splines of order N , which is deſned as [4],∞N2
V = x(·) =c(k)φ(· − k), c ∈  .
k=−∞In this equation, φ(t) is the symmetrical spline of order N ,that is,1, 0 ≤ t ≤ 1,000
φ(t) = (β ∗ · · · ∗ β )(t), β (t) =0, otherwise,
N +1

where ‘∗’ stands for convolution.
The sampled signal x(n), n = 0, 1, 2, . . . of x(t) ∈ V Nis given byx(n) =∞c(k)φ(n − k) = (c ∗ φ)(n).

(1)

k=−∞

On the other hand, the fast sampled signal xL (n) := x(n/L),
n = 0, 1, 2, . . . is given by [4]xL (n) = (cL ∗ φL ) (n),(2)

where φL (n) := φ(n/L) and cL (n) := {(↑ L)c} (n). By (1)and (2), the spline interpolation system is composed of twoſlters ψ and φL , and the upsampler ↑ L as shown in Fig. 1,
where ψ is a system such thatψ ∗ φ = I.

(3)

This is for the perfect reconstruction without any delay. If weallow a delay d ≥ 0 for reconstruction, the condition becomesIII ­ 1469

ψ ∗ φ = z −d .

(4)

ICASSP 2007

x(k)

c(k)

ψ↑LcL (k)φLxL (k)z −d1−φ(z)0Fig. 1. Spline interpolation
ψ(z)3. CAUSAL SPLINE INTERPOLATION BY H ∞OPTIMIZATION
3.1. Non-causal interpolation by decomposition

Fig. 2. H ∞ optimization

The N th-order spline φ(t) is supported in [0, N + 1), andhence the sampled signal φ(n) or φL (n) is represented as anFIR (ſnite impulse response) ſlter. For example, in the caseof N = 3 (cubic spline), we haveφ(z) =

1 2 −1 1 −2+ z + z .
6 36(5)

By (3), the ſlter ψ(z) is the inverse ψ = φ−1 and given by
ψ(z) = 6z 2 /(z 2 + 4z + 1). One of the poles of ψ(z) lies outof the unit circle, and hence the ſlter ψ(z) is unstable. Thesame thing is said of the other N th-order splines [5]. A practical way to implement the ſlter is to decompose ψ(z) into acascade of stable causal and anti-causal ſlters [5]. In the caseof the cubic spline, we ſrst shift the impulse response of (5)as φ(z) = (1/6)z + (2/3) + (1/6)z −1 , and then decompose
ψ(z) = φ(z)−1 as
16α
1−1 ,+1 − α2 1 − αz −1
1 − αz√where α = −2 + 3. Since |α| < 1, this is a stable andnon-causal IIR (inſnite impulse response) ſlter.
ψ(z) = −3.2. Causal interpolation by H ∞ optimization
In image processing, causality is not necessary, and the noncausal ſlter mentioned above is used widely in that ſeld. However, it is difſcult to use such non-causal ſlters for real-timeprocessing, for example, in instrumentation or audio/speechprocessing. We therefore propose designing a causal ſlter
ψ(z) which approximates the condition (4) of delayed perfectreconstruction. Our problem is formulated as follows.
Problem 1 Given a stable transfer function φ(z) and delayd ≥ 0, ſnd the causal and stable ſlter ψ(z) which minimizes
J(ψ) = z −d − ψ(z)φ(z)∞=sup |e−jdθ − ψ(ejθ )φ(ejθ )|.

(6)

3.3. H ∞ optimal cubic spline
The cubic spline (N = 3) is widely used because of its simple
structure. We here give the H ∞ optimal ſlter ψ in a closed
form in the case of the cubic spline.
Deſne E(z) := z −d −ψ(z)φ(z). Substituting (5) into this−d2equation, we have E(z)√ = z −ψ(z)(z−α√1 )(z−α2 )/(6z )where α1 = −2 − 3 and α2 = −2 + 3. This equation
gives ψ(z) = 6z 2 (z −d − E(z))/{(z − α1 )(z − α2 )}. Since|α1 | > 1, the ſlter ψ(z) may have a pole which lies out of theopen unit disc D := {z ∈ C : |z| < 1}. It is easily shownthat the ſlter ψ(z) is stable (i.e., all poles of ψ(z) lie in D) ifand only if(7)E(α1 ) = α1−d .
Then our problem is to ſnd a stable E(z) of minimum H ∞norm under the interpolation constraint (7). This is a kind of
Nevanlinna-Pick interpolation problem [9]. By the maximum
modulus principle, we haveE∞ = sup |E(z)| = sup |E(z)| ≥ |E(α1 )| = |α1 |−d .
|z|=1|z|≥1The minimum inſnity norm interpolating function is therefore the constant function E(z) = α1−d , E∞ = |α1 |−d . By
this, we obtain the optimal ψ(z)
ψ(z) =6z 2(z −d − α1−d )(z − α1 )(z − α2 )

d−1

6z 2
=− d dαd−1−k z k .
α1 z (z − α2 ) k=0 1(8)Remark 1 In the case of higher order splines (i.e., N ≥ 4),
the optimal ſlter can be obtained by Nevanlinna algorithm[9]. A closed form solution is however very complicatedwhen N ≥ 4. In that case, the numerical computation mentioned in the previous section is available.

θ∈[0,2π)This is a standard H ∞ optimization problem, and it can beeffectively solved by standard MATLAB routines (e.g., robustcontrol toolbox [1]) by using the design block diagram shownin Fig. 2.

3.4. FIR ſlter design via LMIThe H ∞ optimal ſlter is generally an IIR one. Since the ſlter ψ(z) to be designed is linearly dependent upon the error
system E(z) = z −d − ψ(z)φ(z). By this nature, the optimalIII ­ 1470

FIR ſlter with ſxed order is obtained by optimization with alinear matrix inequality (LMI).
We here design the H ∞ optimal ſlter ψ(z) as an FIR oneψ(z) =N

ak z −k .

k=0

A state space representation of this FIR ſlter is given by⎤⎡
010 ... 0 0.. ⎥.
⎢ ..
..
..
..
⎢ .
.
.
. ..
. ⎥
⎥
⎢
⎢ ..
.
.
..
.. 0 0 ⎥⎥⎢ .
ψ(z) = ⎢
⎥ (z)
⎥
⎢ .
..
⎥⎢ ..
.
01⎥⎢
⎣ 0 ... ... ... 0 1 ⎦aN...

...

...

a1BDAφCφ

BφDφ(z),z −d =:AdCdBd0BAC(α) D(α)

z −dL xL − x
L = φL (↑ L)ψI x − φL (↑ L)ψx= φL (↑ L)z −d c − φL (↑ L)ψφc= φL (↑ L)(z −d − ψφ)c.
Since {φ(· − k)}∞k=0 is a Riesz basis [2], there exists a real
number K > 0 which depends on φ and is independent of cand x such that c2 ≤ Kx2 . Finally, we have(z).

(z).

By this, the parameter α to be designed is linearly dependent
on the matrices C(α) and D(α). By using the bounded reallemma, we can describe our design problem as an LMI [10].
Proposition 1 Let γ be a positive number. Then the inequality E(z)∞ < γ holds if and only if there exist a positive
deſnite matrix P > 0 such that⎤⎡ T
AT P B
C(α)TA PA − P⎣ BT P A−γI + B T P B D(α)T ⎦ < 0. (9)
C(α)
D(α)
−γI

L 2 ≤ φL (↑ L)∞ z −d − ψφ∞ c2
z −dL xL − x≤ φL (↑ L)∞ z −d − ψφ∞ Kx2= Cz −d − ψφ∞ x2 ,where C := KφL (↑ L)∞ .

By this proposition, we conclude that if the H ∞ normof the error system z −d − φ(z)ψ(z) is adequately small, theSNR of the interpolator can be increased, and hence H ∞ optimization provides a good approximation of the ideal (i.e.,non-causal) spline interpolation.
5. DESIGN EXAMPLE
We here present a design example of causal spline interpolation. We consider the spline of order N = 3 (cubic spline).
We here take the reconstruction delay d = 3 and design theH ∞ optimal IIR ſlter by (8) and FIR one with 5 taps by thelinear matrix inequality (9). In the case of the cubic spline,
the H ∞ optimal IIR ſlter (8) with d = 3 is given by
ψ(z) =Remark 2 Zeros of E(z) can be set by
C(α)(zi I − A)−1 B + D(α) = 0,(10)Proof. Let ψI be the ideal ſlter which satisſes ψI ∗ φ = z −d .
Then, by (1), we have ψI x = z −d c, and(z) := C(zI − A)−1 B + D.

Then, a state space representation of the error system E(z) :=z −d − ψ(z)φ(z) is given by⎤⎡AψBψ Cφ0−Bψ Dφ⎢⎥0Aφ0−Bφ⎥ (z)E(z) = ⎢⎦⎣00AdBdCψ (α) Dψ (α)Cφ Cd −Dψ (α)Dφ=:L 2
z −dL xL − x≤ Cz −d − ψ(z)φ(z)∞ .
x2a0Note that the parameter α to be designed is linearly dependent
on the matrices Cψ (α) and Dψ (α). Set state space representation of φ(z) and z −d respectively by
φ(z) =:In the previous section, we have proposed the H ∞ optimization design of the ſlter φ(z) which approximates the delayed
perfect reconstruction condition (4). In this section, we analyze the overall performance of the interpolation system shownin Fig. 1. We here show that the approximation of the equation (4) is proper for increasing the SNR (signal-to-noise ratio) of the interpolation system.
Proposition 2 Assume that φ(z) and ψ(z) are causal andL be the output of the approximatedstable, and x ∈ 2 . Let xinterpolation system, that is, xL := φL (↑ L)ψx. Then there
exist a real number C > 0 which depends only on φ and Lsuch thatAψBψ=:(z),Cψ (α) Dψ (α)
where α := aN . . . a1 a0 , and we use the notationby Doyle [11]:AC

4. SNR PERFORMANCE ANALYSIS

i = 1, 2, . . . , .

This is a linear matrix equation. The LMI (9) combined withthis linear constraint is also easily solved.

−6z 2 − 6α1 z − 6α12.
α13 z(z − α2 )For comparison, we also design 5-tap FIR ſlter by the constrained least square design (CLSD) [7] and the Kaiser windowed approximation (KWA) [8]. Table 1 shows the coefſcients of the H ∞ optimal FIR ſlter, the ſlters by CLSD andIII ­ 1471

k01234

Table 2. H ∞ norm of E(z)MethodE∞H ∞ optimal IIR 0.019238
H ∞ optimal FIR 0.038597
CLSD [7]0.053446
KWA [8]0.16348

Table 1. Coefſcient ak of FIR ſlter ψ(z)H ∞ optimal
CLSD [7]
KWA [8]0.1152359
0.0991561
0.06049527-0.4614954 -0.4599156 -0.37739071
1.7307475
1.7215190
1.63379087-0.4614951 -0.4599156 -0.37739071
0.1152352
0.0991561
0.06049527













higher order optimal ſlter can be effectively solved. We havealso shown that the H ∞ optimal FIR ſlter can be designedby an LMI. We have shown a design example to illustrate theresult.
7. REFERENCES




1234

[1] G. Balas, R. Chiang, A. Packard, and M. Safonov, Robust Control Toolbox Version 3, The MathWorks, 2005.
*/0+.
,



+,-*&'()[2] G. Strang and T. Nguyen, Wavelets and Filter Banks,Wellesley-Cambridge Press, 1996.


!

"



[3] M. Unser, Sampling — 50 years after Shannon, Proceedings of the IEEE, Vol. 88, No. 4, pp. 569–587, 2000.

%

Fig. 3. Magnitude plot of E(z): H ∞ optimal IIR (solid),
H ∞ optimal FIR (dash), CLSD [7] (dash-dots), and KWA [8](dots).
by KWA. Fig. 3 shows the magnitude of the frequency response of E(z) = z −3 − φ(z)ψ(z). By this ſgure, we cansee that the H ∞ optimal IIR ſlter has an allpass characteristic. The H ∞ optimal FIR ſlter shows almost the same as
the CLSD ſlter except for the zero frequency. This is because
CLSD aims at exact inversion for DC signals. At the price ofthat, the CLSD ſlter shows much error in the high frequency.
The KWA ſlter shows the same nature. Table 2 shows theH ∞ norm of the error system E(z). In general, the purposeof H ∞ design is to minimize the error in the worst case. Itfollows that the system optimized by H ∞ design is robust
against uncertainty of input signals. This is efſciency of H ∞design. On the other hand, CLSD is a design to minimize themean value of the error. Which is better depends on its application. If we do not have much information of the input, H ∞design (worst case optimization) is more suitable.
6. CONCLUSION
In this article, a design of causal interpolation with polynomial splines has been proposed. The design is formulated as
H ∞ optimization. In the case of the cubic spline, the optimal solution is given in a closed form. By using MATLAB,[4] M. Unser, A. Aldroubi, and M. Eden, B-Spline signal
processing: Part-I — Theory, IEEE Trans. Signal Processing, Vol. 41, No. 2, pp. 821–833, 1993.
[5] M. Unser, A. Aldroubi, and M. Eden, B-Spline signal
processing: Part-II — Efſcient design and applications,
IEEE Trans. Signal Processing, Vol. 41, No. 2, pp. 834–848, 1993.
[6] M. Unser and T. Blu, Cardinal exponential splines: PartI — Theory and ſltering algorithms, IEEE Trans. Signal
Processing, Vol. 53, No. 4, pp. 1425–1438, 2005.
[7] M. Unser and M. Eden, FIR approximations of inverseſlters and perfect reconstruction ſlter banks, SignalProcessing, Vol. 36, pp. 163–174, 1994.
[8] B. Vrcelj and P. P. Vaidyanathan, Efſcient implementation of all-digital interpolation, IEEE Trans. Image Processing, Vol. 10, No. 11, pp. 1639–1646, 2001.
[9] J. L. Walsh, Interpolation and Approximation by Rational Functions in the Complex Domain, 5th ed., American Mathematical Society, 1969.
[10] Y. Yamamoto, B. D. O. Anderson, M. Nagahara, and Y.
Koyanagi, Optimizing FIR approximation for discretetime IIR ſlters, IEEE Signal Processing Letters, Vol. 10,
No. 9, pp. 273–276, 2003.
[11] K. Zhou, J. C. Doyle, and K. Glover, Robust and Optimal Control, Prentice Hall, 1995.

III ­ 1472