Causal Spline Interpolation by H∞ Optimization - 京都大学

H ∞ 最適化による因果的スプライン補間
Causal Spline Interpolation by H ∞ Optimization
永原正章 ∗ ，和田敏裕 ∗∗ ，山本裕 ∗∗∗
Masaaki NAGAHARA, Toshihiro WADA, Yutaka YAMAMOTO
京都大学大学院情報学研究科
Graduate School of Informatics, Kyoto University
∗
[email protected],
∗∗
[email protected],
∗∗∗
[email protected]
timization. We give a closed form solution to the
H ∞ 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 補間システムは非因果的となる．これは，Shannon
unstable zeros. For this purpose, we adopt H ∞ op- のサンプリング定理による信号補間と同様である [3]．
画像処理では，非因果性は問題とはならないので，ス
x(k)
プライン補間がよく用いられる．しかし，スプライ
ψ
c(k)
↑L
cL (k)
φL
xL (k)
ン補間を計測や音声処理などのリアルタイム処理に
用いる場合，非因果性は問題となる．
図 1: Spline interpolation
そこで，本論文では，非因果的なスプライン補間
を近似する因果的なシステムを設計する．この設計
:= φ(n/L)，
テムを設計する問題として定式化される．この問題 cL (n) := {(↑ L)c} (n) である．(1) と (2) を用いる
を解くために，本論文では，H ∞ 最適化を用いる．と，スプライン補間システムは，2 つのフィルタ ψ,
これにより，H ∞ 最適な安定逆システムが得られ， φL およびアップサンプラ ↑ L を用いて，図 1 のよ
因果的なスプライン補間システムが設計できる．さうに構成される．ここで，フィルタ ψ は
らに，フィルタを FIR と仮定することにより，設計
ψ ∗ φ = I.
(3)
問題は線形行列不等式として記述でき，MATLAB
等の数値計算ソフトウェアを用いて，効率的に最適を満たす φ の逆フィルタである．(3) は，遅れなし
解を計算することができる．
の完全再構成条件であるが，遅れ d ≥ 0 を許した完
問題は，不安定零点を持つシステムの安定な逆シス
と計算される [4]．ここで φL (n)
全再構成条件は，
2
スプライン補間
ψ ∗ φ = z −d
信号 x ∈ V N を考える．ここで V N は次数 N の
(4)
となる．
多項式スプラインの空間であり，次式で定義される
[4]．
3
∞
VN =
x(·) =
c(k)φ(· − k), c ∈
2
.
k=−∞
この定義式において φ(t) は次数 N の対称スプラ
イン

1,
β 0 (t) =
0,
3.1
H ∞ 最適化による因果的スプラ
イン補間
フィルタの分解による非因果的スプラ
イン補間
0 ≤ t ≤ 1,
N 次スプライン φ(t) は [0, N + 1) に台を持ち，
otherwise, そのサンプル値 φ(n) または φL (n) は FIR フィル
N +1
タで表現できる．例えば，N = 3 （キュービックス
である．ここで ‘∗’ は畳み込み積をあらわす．
プライン）の場合，
信号 x(t) ∈ V N をサンプル周期 T = 1 でサンプ
1
1 2
ルした x(n), n = 0, 1, 2, . . . は，
(5)
φ(z) = + z −1 + z −2 .
6 3
6
φ(t) = (β 0 ∗ · · · ∗ β 0 )(t),
∞
c(k)φ(n − k) = (c ∗ φ)(n)
x(n) =
(1)
である．完全再構成条件 (3) より，フィルタ ψ(z) を
計算すると，
k=−∞
と求められる．いっぽう，サンプル周期 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) =
1
2 1
z + + z −1 ,
6
3 6
まず，誤差システムを 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) = φ(z)−1 を次のように分解する．
ψ(z) =
1
1
+
−1 ,
−1
となる．|α1 | > 1 であるので，フィルタ ψ(z) は単
1 − αz
1 − αz
√
位円 D := {z ∈ C : |z| < 1} の外に極を持つ可能性
ここで， α = −2 + 3 である．|α| < 1 であるの
がある．簡単な計算により，フィルタ ψ(z) が安定
で，これは，安定な非因果的 IIR フィルタである．
である（すなわち，すべての極が D にある) ための
必要十分条件は，
ψ(z) = −
3.2
6α
1 − α2
6z 2
(z −d − E(z))
(z − α1 )(z − α2 )
H ∞ 最適化による因果的スプライン
補間
前節で述べたような非因果的な補間システムは，
リアルタイム処理には不向きである．そこで，本節
では，遅れを許した完全再構成条件 (4) を H ∞ ノ
ルムの意味で最適に近似する ψ(z) を求める．設計
問題は以下のように定式化される．
Problem 1 安定な伝達関数 φ(z) と遅れ定数 d ≥ 0
が与えられたとき，次の評価関数を最小化する因果
的かつ安定なフィルタ ψ(z) を求めよ．
J(ψ) = z −d − ψ(z)φ(z)
=
∞
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]．
複素関数の最大値原理より，
E
∞
= sup |E(z)| = sup |E(z)| ≥ |E(α1 )| = |α1 |−d
|z|=1
|z|≥1
であることがわかり，最小ノルムを持つ補間関数は
E(z) = α1−d , E ∞ = |α1 |−d であることがわかる．
これより，最適なフィルタ ψ(z) は，
ψ(z) =
6z 2
(z −d − α1−d )
(z − α1 )(z − α2 )
d−1
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) の状態空間表現は以下で与えられる．


設計
A
B C
0
−B D


これまで述べた手法により求まる H 最適フィ E(z) = 


ルタは一般に IIR フィルタである．ここで，設計さ
∞
れるフィルタ ψ(z) は誤差システム E(z) = z −d −
ψ(z)φ(z) に線形に依存している．この性質を用いる
と，次数を固定した H ∞ 最適 FIR フィルタは線形
行列不等式による最適化により求めることができる．
設計されるフィルタ H ∞ を FIR とする．すな
わち，
ψ
φ
ψ
φ
Aφ
0
0
0
Ad
−Bφ
Bd
Cψ (α)
Dψ (α)Cφ
Cd
−Dψ (α)Dφ
A
B
C(α)
D(α)


 (z)


(z).
これより，設計パラメータ α は行列 C(α) および
D(α) に線形に依存することがわかる．有界実補題
を用いることにより，設計問題は以下のように線形
N
ψ(z) =
=:
ψ
0
ak z
−k
行列不等式で記述される [10]．
k=0
とおく．この FIR フィルタの状態空間表現は次式で
与えられる．

0
..
.
..
.
..
.
0
...
...
...
1
0
aN
...
...
...
a1






ψ(z) = 





=:
1
..
.
0
..
.
..
.
...
..
.
..
.
..
.
Aψ
Bψ
Cψ (α)
Dψ (α)
0
..
.
0

0
.. 

. 


0  (z)


0 


1 
a0
Proposition 1 実数 γ > 0 が与えられたとき，不
等式 E(z) ∞ < γ が成り立つ必要十分条件は，正
定行列 P > 0 が存在して，
 T

A PA − P
AT P B
C(α)T


−γI + B T P B D(α)T  < 0
 BT P A
C(α)
D(α)
(9)
が成り立つことである．
Remark 2 誤差システム E(z) の零点は
C(α)(zi I − A)−1 B + D(α) = 0,
(z).
−γI
i = 1, 2, . . . , .
により設定できる．これは線形行列等式であり，線
aN . . . a1 a0 である．また，形行列不等式 (9) と組み合わせることにより，零点
状態空間表現において，次式で定義される Doyle の条件を加味した最適解を容易に得ることもできる．
記法 [11] を用いた．
ここで，α :=
A
B
C
D
(z) := C(zI − A)−1 B + D.
4
SN 比による性能の解析
前節において，遅れを許した完全再構成条件 (4)
ここで，設計パラメータ α は行列 Cψ (α) および
を近似する H ∞ 最適なフィルタ φ(z) の設計法を示
Dψ (α) に線形に依存していることに注意する．フィした．この節では，図 1 に示す補間システムのトー
タルの性能を SN (Signal-to-Noise) 比を用いて解析
ルタ φ(z) および z −d の状態空間表現をそれぞれ
する．
Aφ
Bφ
Cφ
Dφ
z −d =:
Bd
Cd
0
(z) Proposition 2 フィルタ φ(z) および ψ(z) は因果
的かつ安定とし，x ∈ 2 とする．xL を近似補間シス
−d
とおく．このとき，誤差システム E(z) = z
− テムの出力信号とする．すなわち，xL := φL (↑ L)ψx．
φ(z) =:
(z),
Ad
このとき，φ と L にだけ依存する実数 C > 0 が存
在し，
z −dL xL − xL
x 2
2
≤ C z −d − ψ(z)φ(z)
(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]
0
1
2
0.1152359
-0.4614954
1.7307475
0.0991561
-0.4599156
1.7215190
0.06049527
-0.37739071
1.63379087
3
4
-0.4614951
0.1152352
-0.4599156
0.0991561
-0.37739071
0.06049527
z −dL xL − xL = φ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
が存在して， c 2 ≤ K x 2 が成り立つ．これより，
z −dL xL − xL
2
≤ φL (↑ L)
∞
z −d − ψφ
∞
c
2
≤ φL (↑ L)
∞
z −d − ψφ
∞K
x
=C z
−d
− ψφ
∞
ψ(z) =
−6z 2 − 6α1 z − 6α12
α13 z(z − α2 )
となる．また比較のために，制約付き最小二乗法に
2
x 2,
よる設計 (CLSD) [7] および Kaiser 窓関数による打
ち切り近似 (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 := K φL (↑ L)
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 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.
✰✳ ✯✴✵✱
☎ ✁✝ ✞ ✰ ✱✲✯
✫✬ ✭✮
[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)
Method
E ∞
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 exponential
splines: 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 of
inverse filters and perfect reconstruction filter
banks, Signal Processing, Vol. 36, pp. 163–
174, 1994.
[8] B. Vrcelj and P. P. Vaidyanathan, Efficient implementation of all-digital interpolation, IEEE
Trans. 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, IEEE
Signal 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.

Download Report