H ∞ 最適化による因果的スプライン補間
Causal Spline Interpolation by H ∞ Optimization
永原正章 ∗ ，和田敏裕 ∗∗ ，山本裕 ∗∗∗Masaaki NAGAHARA, Toshihiro WADA, Yutaka YAMAMOTO
京都大学大学院情報学研究科
Graduate School of Informatics, Kyoto University∗

nagahara@ieee.org, ∗∗ twada@acs.i.kyoto-u.ac.jp, ∗∗∗ yy@i.kyoto-u.ac.jp

timization. We give a closed form solution to theH ∞ optimization in the case of the cubic spline.
スプライン補間システムには通常，非因果的なフィ For higher order splines, the optimal filter can beルタが用いられる．したがってリアルタイム処理が effectively solved by a numerical computation. We必要な場合，このままでは使いにくい．そこで本論 also show that the design of optimal FIR (Finite
文では，スプライン補間を近似する因果的なシステ Impulse Response) filter is described by an LMIムを設計する．この設計は，不安定零点を持つシス (Linear Matrix Inequality), which can also be efテムの安定な逆システムを求める問題として定式化 fectively solved numerically. A design example isされる．この問題を解くために，本論文では H ∞ 最 presented to illustrate the result.
適化を用いる．まずキュービックスプラインに対する H ∞ 最適化問題の解析解を与える．また高次のスプラインに対しては，数値計算により最適解が容易
1 はじめにに求まることを示す．つぎに，H ∞ 最適 FIR (有限インパルス応答）フィルタの設計問題が線形行列不信号補間には，カーブフィッティングや信号復元，等式として記述できることを示す．最後に設計例を 画像の解像度変換など様々な応用がある．特に，多用いて，提案手法の有効性を示す．
項式スプライン補間 [4, 5] は画像処理において広くあらまし用いられている．また，指数関数スプライン補間 [6]Abstractなど新しい補間方法も提案され，理論，応用ともに盛んに研究されている．

Spline interpolation systems generally contain
多項式スプラインや指数関数スプラインによる補non-causal filters, and hence such systems are diffi- 間では，元の信号が区分的な多項式（指数関数）であ
cult to use for real-time processing. Our objective ると仮定し，スプライン基底による Fourier 係数かis to design a causal system which approximates ら補間点の値を求める．しかし，それらの係数は未来
spline interpolation. This is formulated as a prob- のサンプル値を用いなければ求まらず，したがって
lem of designing a stable inverse of a system with 補間システムは非因果的となる．これは，Shannonunstable zeros. For this purpose, we adopt H ∞ op- のサンプリング定理による信号補間と同様である [3]．

画像処理では，非因果性は問題とはならないので，スx(k)プライン補間がよく用いられる．しかし，スプライ

ψ

c(k)

↑LcL (k)φLxL (k)ン補間を計測や音声処理などのリアルタイム処理に用いる場合，非因果性は問題となる．
図 1: Spline interpolation

そこで，本論文では，非因果的なスプライン補間を近似する因果的なシステムを設計する．この設計:= φ(n/L)，
テムを設計する問題として定式化される．この問題 cL (n) := {(↑ L)c} (n) である．(1) と (2) を用いるを解くために，本論文では，H ∞ 最適化を用いる． と，スプライン補間システムは，2 つのフィルタ ψ,
これにより，H ∞ 最適な安定逆システムが得られ， φL およびアップサンプラ ↑ L を用いて，図 1 のよ因果的なスプライン補間システムが設計できる．さ うに構成される．ここで，フィルタ ψ は
問題は，不安定零点を持つシステムの安定な逆シスと 計 算 さ れ る [4]．こ こ で φL (n)らに，フィルタを FIR と仮定することにより，設計ψ ∗ φ = I.

問題は線形行列不等式として記述でき，MATLAB(3)等の数値計算ソフトウェアを用いて，効率的に最適を満たす φ の逆フィルタである．(3) は，遅れなし解を計算することができる．の完全再構成条件であるが，遅れ d ≥ 0 を許した完全再構成条件は，2スプライン補間ψ ∗ φ = z −d

信号 x ∈ V N を考える．ここで V N は次数 N の(4)となる．

多項式スプラインの空間であり，次式で定義される[4]．
(VN =x(·) =∞X)c(k)φ(· − k), c ∈ `23H ∞ 最適化による因果的スプライン補間

3.1

フィルタの分解による非因果的スプライン補間.

k=−∞

この定義式において φ(t) は次数 N の対称スプライン

1,β 0 (t) =0,0 ≤ t ≤ 1,N 次スプライン φ(t) は [0, N + 1) に台を持ち，otherwise, そのサンプル値 φ(n) または φL (n) は FIR フィル
N +1タで表現できる．例えば，N = 3 （キュービックスである．ここで ‘∗’ は畳み込み積をあらわす．
プライン）の場合，
信号 x(t) ∈ V N をサンプル周期 T = 1 でサンプ
11 2ルした x(n), n = 0, 1, 2, . . . は，(5)φ(z) = + z −1 + z −2 .
6 3
6∞Xc(k)φ(n − k) = (c ∗ φ)(n) (1) である．完全再構成条件 (3) より，フィルタ ψ(z) をx(n) =k=−∞計算すると，φ(t) = (β 0 ∗ · · · ∗ β 0 )(t),|{z}と求められる．いっぽう，サンプル周期 T = 1/L（L は自然数）でサンプルした xL (n) := x(n/L),

n = 0, 1, 2, . . . はxL (n) = (cL ∗ φL ) (n),ψ(z) =

6z 2
z 2 + 4z + 1となるが，ψ(z) の極の一つが単位円外にあり，ψ(z)(2) は不安定である．このような状況は，他の N 次ス

プラインでも同様である [5]．これを回避するための一つの方法は，ψ(z) を因果的なフィルタと反因果的なフィルタに分解する方法である [5]．例えばキュービックスプラインの場合，まず，(5) のインパルス応答を次のようにシフトし，φ(z) =12 1z + + z −1 ,63 6フィルタ ψ(z) = φ(z)−1 を次のように分解する．
µ¶
6α11ψ(z) = −+−1,1 − α2 1 − αz −11 − αzここで， α = −2 +√まず，誤差システムを E(z) := z −d − ψ(z)φ(z) とおく. (5) をこれに代入すると，(z − α1 )(z − α2 ),6z 2√√が得られる．ここで，α1 = −2 − 3，α2 = −2 + 3である．これより，
E(z) = z −d − ψ(z)

ψ(z) =6z 2(z −d − E(z))(z − α1 )(z − α2 )となる．|α1 | > 1 であるので，フィルタ ψ(z) は単位円 D := {z ∈ C : |z| < 1} の外に極を持つ可能性3 である．|α| < 1 であるのがある．簡単な計算により，フィルタ ψ(z) が安定で，これは，安定な非因果的 IIR フィルタである．
である（すなわち，すべての極が D にある) ための必要十分条件は，

3.2

H ∞ 最適化による因果的スプライン補間

前節で述べたような非因果的な補間システムは，リアルタイム処理には不向きである．そこで，本節では，遅れを許した完全再構成条件 (4) を H ∞ ノルムの意味で最適に近似する ψ(z) を求める．設計問題は以下のように定式化される．

Problem 1 安定な伝達関数 φ(z) と遅れ定数 d ≥ 0が与えられたとき，次の評価関数を最小化する因果的かつ安定なフィルタ ψ(z) を求めよ．
J(ψ) = kz −d − ψ(z)φ(z)k∞=sup |e−jdθ − ψ(ejθ )φ(ejθ )|.

(6)

θ∈[0,2π)

この問題は，標準的な H ∞ 最適化問題（モデルマッチング問題）であり，MATLAB などの数値計算ソフトウェアを用いて効率よく解を得ることができる[1]．

E(α1 ) = α1−d(7)であることがわかる．したがって，設計問題は，補間条件 (7) を満たす E(z) のなかで，H ∞ ノルムが最小の安定な E(z) を見つける問題に帰着する．この問題は，Nevanlinna-Pick 補間問題と呼ばれる [9]．
複素関数の最大値原理より，

kEk∞ = sup |E(z)| = sup |E(z)| ≥ |E(α1 )| = |α1 |−d|z|=1|z|≥1であることがわかり，最小ノルムを持つ補間関数はE(z) = α1−d , kEk∞ = |α1 |−d であることがわかる．
これより，最適なフィルタ ψ(z) は，
ψ(z) =6z 2(z −d − α1−d )(z − α1 )(z − α2 )

d−1X
6z 2
αd−1−k z k=− d dα1 z (z − α2 ) k=0 1(8)と求められる．

H ∞ 最適キュービックスプライン

Remark 1 高次スプラインの場合，最適なフィル
タは Nevanlinna アルゴリズムを用いることによっ

ここでは，実装が簡単であるため最もよく使われて求められる [9]．しかし，その数式表現は非常に複るキュービックスプライン（N = 3）について，そ雑である．このような場合は，前節で述べたようなの最適解を解析的に与える．

数値計算による最適解の計算が有効である．

3.3

3.4

線形行列不等式による FIR フィルタ ψ(z)φ(z) の状態空間表現は以下で与えられる．

設計AB C0−B D
これまで述べた手法により求まる H 最適フィ E(z) = ルタは一般に IIR フィルタである．ここで，設計さ∞れるフィルタ ψ(z) は誤差システム E(z) = z −d −ψ(z)φ(z) に線形に依存している．この性質を用いると，次数を固定した H ∞ 最適 FIR フィルタは線形行列不等式による最適化により求めることができる．
設計されるフィルタ H ∞ を FIR とする．すな
わち，

ψ(z) =NX"=:

ψψ

0Aφ000AdCψ (α)AC(α)φψ

Dψ (α)Cφ Cd#B(z).
D(α)φ−BφBd (z)

−Dψ (α)Dφ

これより，設計パラメータ α は行列 C(α) およびD(α) に線形に依存することがわかる．有界実補題
を用いることにより，設計問題は以下のように線形

ak z−k行列不等式で記述される [10]．

k=0とおく．この FIR フィルタの状態空間表現は次式で与えられる．

0..
.
..
.
..
.
0...

...

...

aN

...

...

. . . a1#ψ(z) = 


"
=:

1..

.

0..

.

..

.

...
..
.
..
.
..
.

Aψ

Bψ

Cψ (α)Dψ (α)

h0..
.
010

0.. 
. 0  (z)0 
1 a0Proposition 1 実数 γ > 0 が与えられたとき，不等式 kE(z)k∞ < γ が成り立つ必要十分条件は，正定行列 P > 0 が存在して，
 T
A PA − P
AT P BC(α)T−γI + B T P B D(α)T  < 0 BT P A
C(α)

D(α)

−γI(9)が成り立つことである．

Remark 2 誤差システム E(z) の零点はC(α)(zi I − A)−1 B + D(α) = 0,(z).

i = 1, 2, . . . , `.

により設定できる．これは線形行列等式であり，線iaN . . . a1 a0 である．また， 形行列不等式 (9) と組み合わせることにより，零点状態空間表現において，次式で定義される Doyle の 条件を加味した最適解を容易に得ることもできる．
記法 [11] を用いた．
ここで，α :="ABCD#4(z) := C(zI − A)−1 B + D.

SN 比による性能の解析
前節において，遅れを許した完全再構成条件 (4)

ここで，設計パラメータ α は行列 Cψ (α) およびを近似する H ∞ 最適なフィルタ φ(z) の設計法を示Dψ (α) に線形に依存していることに注意する．フィ した．この節では，図 1 に示す補間システムのトータルの性能を SN (Signal-to-Noise) 比を用いて解析
ルタ φ(z) および z −d の状態空間表現をそれぞれ"
φ(z) =:

Aφ

Bφ

Cφ

Dφ

#"
(z),z −d =:Ad

Bd

Cd0とおく．このとき，誤差システム E(z) = z#する．

(z) Proposition 2 フィルタ φ(z) および ψ(z) は因果的かつ安定とし，x ∈ `2 とする．xeL を近似補間シス−d− テムの出力信号とする．すなわち，xeL := φL (↑ L)ψx．

このとき，φ と L にだけ依存する実数 C > 0 が存在し，kz −dL xL − xeL k2≤ Ckz −d − ψ(z)φ(z)k∞kxk2(10)が成り立つ．

Proof. ψI を完全再構成条件 ψI ∗ φ = z −d を満たす理想フィルタとする．このとき，(1) より，ψI x =z −d c が成り立ち，

表 1: Coefficient ak of FIR filter ψ(z)k H ∞ optimal CLSD [7]
KWA [8]012

0.1152359-0.4614954
1.7307475

0.0991561-0.4599156
1.7215190

0.06049527-0.37739071
1.63379087

34-0.4614951
0.1152352-0.4599156
0.0991561-0.37739071
0.06049527

z −dL xL − x
eL = φL (↑ L)ψI x − φL (↑ L)ψx= φL (↑ L)z −d c − φL (↑ L)ψφc

ルタ (8) は= φL (↑ L)(z −d − ψφ)cとなる．{φ(· − k)}∞k=0 は Riesz 基底 [2] であるので，φ に依存し，c および x には依存しない実数 K > 0が存在して，kck2 ≤ Kkxk2 が成り立つ．これより，

ψ(z) =−6z 2 − 6α1 z − 6α12
α13 z(z − α2 )kz −dL xL − x
eL k2 ≤ kφL (↑ L)k∞ kz −d − ψφk∞ kck2となる．また比較のために，制約付き最小二乗法に≤ kφL (↑ L)k∞ kz −d − ψφk∞ Kkxk2 よる設計 (CLSD) [7] および Kaiser 窓関数による打= Ckz −d − ψφk∞ kxk2 ,ち切り近似 (KWA) [8] による設計も行う．表 1 にH ∞ 最適 FIR フィルタおよび CLSD と KWA により設計した FIR フィルタの係数を示す．また，図 2¤に誤差システムE(z) = z −3 − φ(z)ψ(z) の周波数応
この命題より，誤差システム E(z) = z −d −答を示す．この図より，H ∞ 最適 IIR フィルタは全φ(z)ψ(z) の H ∞ ノルムが十分小さければ，補間シ周波数でフラットな特性を持つことがわかる．H ∞ステムの SN 比は大きくなることがわかり，したがっ最適 FIR フィルタは CLSD によるフィルタと似たて H ∞ 最適化により理想的な（すなわち，非因果的応答を示している CLSD と KWA はともに，直流入な）補間システムの性能に近い近似システムが得ら力に対する誤差を 0 にする，すなわち誤差システムれることがわかる．
E(z) が z = 1 に零点を持つように設計されていることに注意する．その代償として，高周波側で誤差が成り立つ．ここで，C := KkφL (↑ L)k∞ である．

5

設計例が大きくなっていることがわかる．表 2 にそれぞれの設計における誤差システム E(z) の H ∞ ノルムを

この節では，前節までで示した H ∞ 最適化による示す．一般に H ∞ 設計は，最悪入力に対する応答を因果的スプライン補間システムの設計例を示す．ス

最適化するので，入力信号の不確かさに対してロバプラインの次数は 3 とする（すなわち，キュービッ ストであり，これが H ∞ 設計の利点である．いっぽクスプラインを考える）．再構成の遅れを d = 3 とし，Hう，CLSD は誤差の平均値を最小化する設計である．

∞最適な IIR フィルタを (8) により，また， どちらが良いかは，応用によって異なる．もし入力5 タップの FIR フィルタを線形行列不等式 (9) を用 信号に不確かさが存在するのなら，H ∞ 設計（最悪いて設計する．N = 3 の場合の H ∞ 最適 IIR フィ ケースの最適化）が適している．

       [2] G. Strang and T. Nguyen, Wavelets and FilterBanks, Wellesley-Cambridge Press, 1996.
[3] M. Unser, Sampling — 50 years after Shannon,
Proceedings of the IEEE, Vol. 88, No. 4, pp.

 
1234

569–587, 2000.

+. */0, + ,-*&' ()[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.

      !#"$%

図 2: Magnitude plot of E(z): H ∞ optimal IIR(solid), H ∞ optimal FIR (dash), CLSD [7] (dashdots), and KWA [8] (dots).
表 2: H ∞ norm of E(z)MethodkEk∞

H ∞ optimal IIR
H

6∞optimal FIR
CLSD [7]
KWA [8]0.019238
0.038597
0.053446
0.16348

むすび
本論文では，因果的な多項式スプライン補間シス[5] M. Unser, A. Aldroubi, and M. Eden, B-Spline
signal processing: Part-II — Efficient design
and applications, IEEE Trans. Signal Processing, Vol. 41, No. 2, pp. 834–848, 1993.
[6] M. Unser and T. Blu, Cardinal exponentialsplines: Part I — Theory and filtering algorithms, IEEE Trans. Signal Processing, Vol.
53, No. 4, pp. 1425–1438, 2005.
[7] M. Unser and M. Eden, FIR approximations ofinverse filters and perfect reconstruction filterbanks, Signal Processing, Vol. 36, pp. 163–174, 1994.
[8] B. Vrcelj and P. P. Vaidyanathan, Efficient implementation of all-digital interpolation, IEEETrans. Image Processing, Vol. 10, No. 11,
pp. 1639–1646, 2001.

テムを H ∞ 最適化により設計する手法を提案した． [9] J. L. Walsh, Interpolation and Approximation
by Rational Functions in the Complex DoH ∞ 最適化により，入力信号の不確かさに対してロバストな補間システムが設計できる．また，本論分で述べた手法は，指数関数スプライン [6] においても同様に扱うことができる．

参考文献

main, 5th ed., American Mathematical Society, 1969.
[10] Y. Yamamoto, B. D. O. Anderson, M. Nagahara, and Y. Koyanagi, Optimizing FIR approximation for discrete-time IIR filters, IEEESignal Processing Letters, Vol. 10, No. 9, pp.

273–276, 2003.
[1] G. Balas, R. Chiang, A. Packard, and M. Safonov, Robust Control Toolbox Version 3, The [11] K. Zhou, J. C. Doyle, and K. Glover, Robust
MathWorks, 2005.
and Optimal Control, Prentice Hall, 1995.