数理解析研究所講究録
第 1869 巻 2013 年 1-14

1

連続信号に対するスパースサンプリング
山口大学・大学院医学系研究科

平林晃

Akira Hirabayashi
Graduate School of Medicine
Yamaguchi University
概要

連続信号に対する標本化と言えばこれまで，ナイキスト間隔で行われることが常識的となってきた．しかし，レーダーやソナー，エコーなどに現れる信号は広帯域でありながら，特定波形を平行移動したものの線形結合になっていることに着目すれば，ナイキスト間隔より遥かに広いスパースサンプリングが可能になる．また，区間ごとに多項式や指数関数で定義された信号は帯域制限を受けてすらいないが，同様のスパースサンプリングが可能である．本稿では，これらの信号を包含する不確定率有限信号とよばれる信号クラスを紹介し，無雑音と雑音を含む場合のそれぞれに対して，標本化と再構成の理論を概説する．
特に雑音を含む場合にこの問題が行列の構造制約付き低ランク近似 (Strucmred
Low-Rank Approximation) と呼ばれる問題に帰着され，画像処理や最適化理論
で注目集める分野に結びついていく事を示す．

1

ナイキストとスパース標本化
広く知られているように，信号 $s(t)$ にフーリエ変換$\hat{s}(\omega)=\int_{\infty}^{\infty}s(t)e^{-i\omega t}dt$が存在し，そのサポートに関して条件(1)$\hat{s}(\omega)=0(|\omega|\geq\omega_{。})$

が成立している場合，$\pi/\omega_{c}$

(2)

以下の間隔

$T$

で得られた標本値から信号 $s(t)$ を$s(t)= \frac{\omega_{c}T}{\pi}\sum_{k=-\infty}^{\infty}f(kT)$

によって完全再構成できる [1], [2]1

sinc

sinc

$\{\frac{\omega_{c}(t-kT)}{\pi}\}$

ここで，$(t)=\{\begin{array}{ll}\sin(\pi t)/(\pi t) (t\neq 0)1 (t=0)\end{array}$

1 この定理を戦後混乱期の日本で独自に発見した染谷勲の事績については文献 [3], [4] に詳しい．

2である．条件 (1) は信号処理の分野で「低域通過型の帯域制限」と呼ばれている．間隔$\pi/\omega_{c}$はナイキスト間隔と呼ばれ，アナログ信号の標本化に関するほぼ唯一の指導原理となっている．ここで周波数 。が小さな値であればナイキスト間隔は大きな値となり，それほど問題にはならない．しかし， が大きな値である場合，ナイ
キスト間隔は非常に小さな値となり，種々の問題を引き起こす．すなわち，膨大なデータ量やそれに伴う計算量，あるいは標本化にかかる時間や実装するハードウエ
アのコストである．これらの問題を回避する為に， が大きくナイキスト間隔が非常に小さい場合であっても，より広い間隔で得られた標本値から信号を完全再構成する研究が従来から行われてきた．例えば，信号 $s(t)$ が帯域通過型の帯域制限を受けている場合，即ち条件 (1) に加えて$\omega$$\omega_{c}$$\omega_{c}$(3)$\hat{s}(\omega)=0(|\omega|\leq\omega_{l}<\omega_{。})$も成立している場合を考える．条件 (3) を見過ごしてしまい，この信号を低域通過型の帯域制限を受けていると見なしてしまうと，完全再構成の為の標本間隔は非常に小さな値になってしまう．しかし，条件 (3) をきちんと活用すれば，標本間隔を飛躍的に拡大できるのである [5]. 見過ごすなどという事は考えがたいと思われるかもしれないが，工学の問題では起こり得るのである．実際，光計測の分野である白色光干渉法による表面凹凸形状測定では，観測対象信号は帯域通過型帯域制限を受けているにも関わらず，低域通過型と見なされ，非常に細かいナイキスト間隔で標本化が行われてきた [6].

これに対して筆者等は，帯域通過型標本化定理を同技術に適用し，当時として世界最高速の測定アルゴリズムを開発した [7]. 装置に用いる光学フィルターにも依存するが，測定精度を同程度に保ちながら，標本間隔を 6 倍から 14 倍に拡大する事に成功した．また当該技術の特許を日本，アメリカなどで取得し，そのアルゴリズムを搭載した計測装置を販売した [8]. このように，信号の性質を有効に活用すれば，ナイキスト間隔による密な標本化ではなく，より広い間隔による疎な標本化 (スパースサンプリング) が可能になるのである．

最近，このような標本化理論に対する研究が盛んに行われるようになってきた．不確定率有限信号 2 (Signals with Finite Rate of Innovation, FRI 信号) に対する標本化理論と呼ばれているものがそれである [9]. 不確定率有限信号の典型例は反射波の重畳で構成される信号である．レーダーやエコー，ソナーなどでは，対象物との距離を計測する為に電磁波や超音波などを送出し，対象物体から得られる反射波の遅延時間を測定する．送出波は一般に広帯域 ( が大) であるので，標準的手法である相関関数に基づく手法では，非常に狭いナイキスト間隔で得られた標本に基づいて遅延時間を推定する [10]. しかし，送出波形は既知である事を考慮すれば，未知$\omega_{c}$パラメータは遅延時間と減衰率のみであり，これらの情報を得る為だけにナイキス
ト間隔による密な標本化を行う事は冗長である．できることならば，未知パラメー
タである遅延時間と減衰率の出現頻度に近い間隔で標本化を行いたい．これによつ
2 訳語は東京工業大学小川英光名誉教授のご提案による．

3て，連立方程式の未知数と条件数が一致する場合のように，測定データから未知パ
ラメータを非線形な手法で計算できるのではないかと考えるのである．このことに着目した理論が不確定率有限信号に対する標本化理論である．本稿ではこの理論を概説する．

本稿は以下のように構成されている．第 2 節で不確定率有限信号を定義する．本稿では特に周期信号を論じる．第 3 節では，不確定率に非常に近い頻度で得られた
無雑音の標本値から FRI 信号を完全再構成できることを示す．再構成手法は非線形になっている．第 4 節では，標本値に雑音が含まれている場合に標準的に用いられている雑音抑制法とその問題点を説明し，著者等が提案している最尤推定による再構成手法を述べる．第 5 節は結論である．

2 不確定率有限信号式 (2) を一般化する事から始める．式 (1) が成立しているということは，信号が$\frac{\omega_{c}T}{\pi}$sinc$( \frac{\omega_{c}t}{\pi})$を平行移動した関数の線形結合で表現されることを意味している．式 (2)においては平行移動量が $kT$ であり，結合係数が標本値 $f(kT)$ になっている．これらを次のように一般化する．信号が既知の関数$\varphi(t)$を平行移動した関数の線形結合で表現されると分かっている場合を考える．ただし，平行移動量と係数は未知であ$k<k’$ のときを仮定
り，それぞれ および で表す．一般性を失う事無く，する．このとき信号 $s(t)$ は，$t_{k}$$t_{k}<t_{k’}$$c_{k}$(4)$s(t)= \sum_{k=-\infty}^{\infty}c_{k}\varphi(t-t_{k})$および同一のに含まれると表現されることになる．区間
個数の総和をで表し，不確定率 (rate of innovation) を$[t_{a},t_{b}]$$t_{k}$$C_{s}(t_{a},t_{b})$

(5)

$k$

に対する

$c_{k}$の$\rho$$\rho=\lim_{\tauarrow\infty}\frac{1}{\tau}C_{s}(-\tau/2, \tau/2)$と定義する．

Definition 1 [9] 式 (4) でパラメトリック表現される信号 $s(t)$ に対して，式 (5) で定$s(t)$ を不確定率有限信号 $(a$ signal with
義される不確定率 が有限であるとき，finiteinnovation)と呼ぶ．
rate of$\rho$

式 (5) は不確定率を大域的に定義しているが，局所的に定義する事もできる．すなわち，時刻(6)$t$を中心とする固定長$\tau$の区間を考え，$\rho_{\tau}(t)=\frac{1}{\tau}C_{s}(t-\tau/2, t+\tau/2)$

4
によって局所的不確定率を定義する．この場合は更に， に関する最大値$t$$\rho_{\max}(\mathcal{T})=\max_{t\in R}\rho_{\mathcal{T}}(t)$が重要となる．信号が周期 を有する場合には局所的不確定率が扱いやすい．
不確定率有限信号の中で特に重要なものが，ディラックのデルタ関数列$\tau$(7)$s(t)= \sum_{k=-\infty}^{\infty}c_{k}\delta(t-t_{k})$である．なぜなら，式 (4) で与えられる一般の不確定率有限信号は，式 (7) の $s(t)$ ととの畳込みで表現され，標本化・再構成もディラックのデルタ関数列 $s(t)$ に対する場合と同様に論じる事ができるからである．本稿でも，デルタ関数列の標本化．
再構成問題を論じる事とする．
本稿では更に簡単の為に周期信号を考える．周期 に含まれるデルタ関数の個数を $K(<\infty)$ で表せば，周期的ディラック列 $s(t)$ は，区間 $[0,\tau)$ の信号$\varphi(t)$$\tau$$s_{0}(t)= \sum_{k=0}^{K-1}c_{k}\delta(t-t_{k})(0\leq t_{0}<t_{1}<\ldots<t_{K-1}<\tau)$

を用いて
(8)

$s(t)= \sum_{k’=-\infty}^{\infty}s_{0}(t-k’\tau)$と表すことができる．本稿ではこの信号の標本化・再構成問題を論じる．当然のこ, ナイキスト間隔は無限小ととながら，この信号の帯域は無限大でありと対応する係数なる．一方，式 (8) の 1 周期あたりの自由パラメータは位置$(\omega_{c}=\infty)$$\{t_{k}\}_{k=0}^{K-1}$に依らずに $2K/\tau$ となる．この値に非常に近い頻度で得られた標本値からデルタ関数列が完全再構成されることを次節で示す．
$\{c_{k}\}_{k=0}^{K-1}$

であり，不確定率$\rho_{\tau}(t)$は$t$

なお，信号 $s(t)$ は周期的なのでフーリエ級数で表現され，その係数は(9)となる．ここで$\hat{d}_{p}=\frac{1}{\tau}\int_{0}^{\tau}s(t)e^{-i2p\pi t/\tau}dt=\frac{1}{\tau}\sum_{k=0}^{K-1}c_{k}u_{k}^{p}$$u_{k}=e^{-i2\pi t_{k/\mathcal{T}}}$である．このようにフーリエ係数がべき級数$u_{k}^{p}$の線形結合になっていることが重要である．

3 雑音が含まれない場合の標本化と再構成
周期的デルタ関数列の$t=t_{n}t\ovalbox{\tt\small REJECT}$

おける標本値$s(t_{n})$はほとんどが$0$になってしまい，

このような標本値を集めたところであまり意味が無い．そこで，標本化核$\psi(t)$を $nT$

5
だけ平行移動した関数 $\psi_{n}(t)=\psi(t-nT)$ と $s(t)$ との内積を標本値として用いる :(10)$d_{n}= \langle s, \psi_{n}\rangle=\int_{-\infty}^{\infty}s(t)\psi(t-nT)dt$$n=0\sim N-1$ および $T=\tau/N$ である．厳密に言えば，この標本値を一般
ここで，標本値と呼び， を理想標本値と呼んで区別する [11,12]. ただし本稿では，単に標本値と呼ぶことにする．標本化核の取り方には様々なものがある．例えば，$\psi(t)=$ Bsinc $(Bt)$ とおけば，信号処理で理想ローパスフィルターやアンチェイリアス
フイルターと呼ばれているものになる．この場合，完全再構成の為に条件$s(t_{n})$$\psi(t)$(11)$B \geq\rho=\frac{2K}{\tau}$を仮定する．この他にも，-spline [13] や-spline[14] のようにコンパクトサポー
トの関数を用いることや [15], 周波数領域で定義された Sum of Sinc と呼ばれる関数を用いることがある [16]. 標本化核の選択は本理論の文字通り中核となる部分で
$B$

$E$

ある．本稿では，その簡易性から $\psi(t)=$ Bsinc $(Bt)$ を用いる事とする．

本論文で論じる課題は，標本値$\{d_{n}\}_{n=0}^{N-1}$

から $2K$ 個の未知数を求(Prony の手法とも呼ば$\{t_{k}\}_{k=0}^{K-1}$

めることである．この問題はアニヒレーティングフィルタ
れる) を用いることによって解決できる [9]. 概略を述べる．
整数を で表す :. このとき，ボアソン総和式$B\tau/2$$P$と$\{c_{k}\}_{k=0}^{K-1}$

を超えない最大$P=\lfloor B\tau/2\rfloor$$B \sum_{k’=-\infty}^{\infty}$sinc$\{B(t+k’\tau)\}=\frac{1}{\tau}\sum_{p=-P}^{P}e^{-i2p\pi t/\tau}$より，$d_{n}= \sum_{p=-P}^{P}\hat{d}_{p}e^{i2pn\pi/N}$が成立する．この式は，式 (9) のフーリエ係数$\hat{d}_{p}$が逆離散フーリエ変換によってシンク核によるサンプル値に対応付けられることを意味している．逆変換によってからを求めることができるように，$d_{n}$$d_{n}$$\hat{d}_{p}$(12)$N\geq 2P+1$を仮定すれば，(13)$\hat{d}_{p}=\frac{1}{N}\sum_{n=0}^{N-1}d_{n}e^{-i2pn\pi/N}$を得る．第 番目の要素が$n$が$\hat{d}_{p-P}$(14)$d_{n}$

である $N$ 次元ベクトルを $d$ で表し，第 $p$ 番目の要素である $2P+1$ 次元ベクトルを$\hat{d}$で表せば，式 (13) は$\hat{d}=Fd$6と行列表現される．ここで， は式 (13) に応じて定義された離散フーリエ変換 (DFT)$F$

行列である．
系列$\{\hat{d}_{p}\}_{p=-P}^{P}$

は

$z$

変換が

$A(z)= \sum_{k=0}^{K}a_{k}z^{-k}=a_{0}\prod_{k=0}^{K-1}(1-u_{k}z^{-1})$

であるフィルタ $[a_{0},a_{1}, \ldots,a_{K}]$ によって零変換される :(15)$a_{0}\hat{d}_{p}+a_{1}\hat{d}_{p-1}+\ldots+a_{K}\hat{d}_{p-K}=0(p=K-P, \ldots, P)$

これらの式を連立して解くことによりフィルタ係数が求まる．具体的には，(16)$T=[\hat{d}_{K-P+1}\hat{d}_{K,.....\cdot-P}\hat{d}_{P} \hat{d}_{K-P-1}\hat{d}_{K,.\cdot..\cdot\cdot-P}\hat{d}_{P-1} ..\cdot.

と定義される $(2P-K+1)\cross(K+1)$ の行列$T$\hat{d}_{-P+1}\hat{d}_{P-K}\hat{d}_{-P})$を用いて表される方程式$Ta=0$(17)を満たすベクトル $a=(a_{0}, a_{1}, \ldots, a_{K})$ を求める．ここで， の定義より $P\leq B\tau/2<$$P+1$ であり，条件 (11) より $K\leq B\tau/2$ であるので，
$P$

(18)

$P\geq K$$P>K$ の場合にはが成立している．従って， は $P=K$ の場合に正方行列であり，縦に長い矩形行列になる．列数を考えれば のランクは最大で $K+1$ であるが，$T=USV^{*}$ と特異値となる :rank $(T)=K$ . 従って，が個であるので実際には分解すれば，特異値には 1 個の が含まれる事となり，それに対応する の列ベク
トルが解 を与える．
なお，ここでの再構成手続きには用いられていないが，行列 は Toeplitz 構造を有している．この性質は，rank $(T)=K$ と併せて標本値に雑音が含まれる場合に重
$T$

$T$

$t_{k}$

$K$

$K$

$V$

$0$

$a$

$T$

要な働きをする．
ひとたび が求まれば，式 (9) は行列$t_{k}$$c_{k}$

に関する線形方程式になる．即ち，Vandermonde$U_{t}=(u_{u_{0}^{P}}^{\frac{u}{0}P+1}0^{P}- u_{u_{1}^{P}}^{\frac{u}{1}P+1}-P1\cdots\cdot \cdot.\cdot u_{K-1}^{-P+1}u_{K-1}^{-.P}u_{K-1}^{P}]$7およびベクトル $c=(c_{0}, \ldots, c_{K-1})^{T}$ を用いて表される方程式(19)$U_{t}c=\hat{d}$をもち，最終的に信号 $s(t)$ を完全再構成できるのである．以上をまとめれば以下の標本化定理を得る:
は， が互いに異なるので単一解$t_{k}$$c$

Theorem 1 [9] 標本化核 $\psi(t)=$ Bsinc $(Bt)$ のが条件 (11) を満たし，標本点数 $N$ について条件 (12) が成立しているとき，式 (10) によって得られる標本値
は，$s(t)$
周期のディラック関数列の十分な特徴付けになっている．
$B$$\{d_{n}\}_{n=0}^{N-1}$$\tau$式 (12), (18) より，$N\geq 2K+1$が成立する．この式は，1 周期あたり $2K+1$ 個以上の標本値から信号 $s(t)$ を完全再構成できることを意味している．この値は未知パラメータ数 $2K$ に非常に近い値になっている．標本数が 1 多い理由は式 (13) にある．仮に$\hat{d}_{p}$を直接観測できれば，$2K$

個のそれらの値からを求める事ができる．しかし，観測できる値は であり，これらの値から式 (13) を使ってを求めなければならない．このとき，実数を表現する場合の DFT の共役対称性から係数は奇数個となる．これらの値を求める
為に，最低で $2K+1$ 個の標本値が必要になるのである．
式 (9) の形式をもつフーリエ係数列から やを求める手順は決して新しいものではない．スペクトル推定や到来方向推定などと数学的に等価な問題になっている$t_{k},$$d_{n}$$c_{k}$$\hat{d}_{p}$$\hat{d}_{p}$$t_{k}$$c_{k}$

ので，これらの問題のための従来手法である MUSIC [17] や ESPRIT [18], Matrix
Pencil [191 などを適用する事も可能である．ただし，これらの手法では冗長な標本値を用いる．その一方で，上記の再構成手続きでは完全再構成に必要な標本値の最小数を解明する為に，アニヒレーティングフィルタ法を用いた．このように従来手法を活用しており新規性が無いように捉えられがちであるが， に同手法を適用できる事を見いだし， を標本値 から求められる事を示した点は，文献 [9] の新規$\hat{d}_{p}$$\hat{d}_{p}$$d_{n}$

かつ重要な知見である．

4 雑音が含まれる場合の標本化と再構成
雑音が含まれる標本値をで表す : $y_{n}=d_{n}+e_{n}$ . 第$N$ 次元ベクトルをそれぞれ $y,$で表せば，
$y_{n}$$n$

番目の要素が$y_{n},$$e_{n}$

である

$e$

(20)

$y=d+e$となる．ここで， の確率密度関数 $p(e)$ が既知であることを仮定する．ベクトル$e$の離散フーリエ変換を$\hat{y}=$$(\hat{y}_{-P,\hat{\mathcal{Y}}-P+\iota}, \cdots,\hat{y}_{P})^{T}=Fy$で表す．フーリエ係数$\hat{d}_{p}$

$y$

の代

8
を用いて構成される式 (17) の左辺を にする は一般に存在しない．こを最小にする をの問題の直ちに考えうる解決策は，残差 $Ta$ の 2 乗ノルム
わりに

$0$

$\hat{y}_{p}$$a$$||Ta||^{2}$$a$

解として用いる事である．この解は，標本値に雑音が含まれない場合と同様に，$T$の特異値分解 $USV^{*}$ の最小特異値に対応するの列ベクトルで与えられる．雑音の影響で最小特異値は一般に 0 にはならない．この点のみが，雑音が含まれない場合と異なっている．この方法を最小 2 乗法 ( $LS$ 法) と呼ぶ事にする．
$LS$ 法の精度をより高める為に Cadzow のアルゴリズム [20] が用いられている [211.
このアルゴリズムは，前節で述べたように式 (16) で定義される行列 が rank $(T)=K$
$V$

$T$

を満たし，かつ Toeplitz 構造を有するという 2 種類の性質を利用したものである．行は一般に矩形であるが，正方である場合に精度が高いことが実験的に示されている [21]. そこで以下では正方行列の場合のアルゴリズムを説明する．

列$T$

1. 標本値$y_{n}$

2. 正方行列の離散フーリエ変換$Y$$\hat{y}_{p}$を計算する．

を以下で構成する :$Y=\{\begin{array}{llll}\hat{\mathcal{Y}}o \hat{\mathcal{Y}}-[ \cdots \hat{\mathcal{Y}}- P\hat{y}_{l} \hat{\mathcal{Y}}o \cdots \hat{\mathcal{Y}}- P+1\vdots \vdots \vdots\hat{y}_{P} \hat{y}_{P- 1} \cdots \hat{\mathcal{Y}}o\end{array}\}$3. 終了条件が満たされるまで以下を繰り返す:$U,$ $S,$(a) 特異値分解 $Y=USV^{*}$ を計算する．ここで，$V$

は $P+1$ 次正方行列

である．
(b) 行列の対角成分に含まれている特異値のうち，大きな値を持つものか
ら順に個を残し，それら以外の $P-K+1$ 個の特異値を にした行列$S’$ を用いて行列 $Y’=US’V^{*}$
を計算する．得られた は Toeplitz 構造を
$S$

$K$

$0$

$Y’$

有していない．
(c) 行列 $Y’$ の対角平行成分毎の平均値で得られる行列でを置き換える事により Toeplitz 構造を回復し (a) に戻る．例えば $P=2$ の場合は以下の操作
$Y$

となる :

$Y’=(\begin{array}{lll}y_{0,0}’ y_{0,1}’ y_{0,2}’y_{1,0}’ y_{1,1}’ y_{1,2}’y_{2,0}’ y_{2,1}’ y_{2,2}\end{array}) \Rightarrow Y:=(\begin{array}{lll}\hat{y}_{0}’ \hat{y}_{-1}’ y_{0,2}’\hat{y}_{1}’ \hat{y}_{0}’ \hat{y}_{-1}’y_{2,0}’ \hat{y}_{1} \hat{y}_{0}’\end{array})$

ここで，$\hat{y}_{1}’,\hat{y}_{0}’,\hat{y}_{-1}’$

はそれぞれ以下の値である :$\hat{y}_{1}’=\frac{y_{1,0}’+y_{2,1}’}{2}, \hat{y}_{0}’=\frac{y_{0,0}’+y_{1,1}’+y_{2,2}’}{3}, \hat{\mathcal{Y}}’-1=\frac{y_{0,1}’+y_{1,2}’}{2}$9
終了条件はステップ (b) に現れる特異値の $K$ 番目と $K+1$ 番目の比が一定値以上になった場合などが考えられるが，通常は $10$ $20$ 回程度の繰り返しで が Toeplitz$Y$〜

構造を有し，rank$(Y)=K$ を満たすようになる．このように構造を持った低ランク行列によってもとの行列を近似する問題は構造制約付き低ランク近似 (Stmctured
Low-Rank Approximation) と呼ばれ，近年盛んに議論されている [22].

ここで注意すべき点は，Toeplitz 構造を有する行列の全体は凸集合であるが，ラン
クが $K$ である行列のそれは凸になっていないということである．この結果，このア

ルゴリズムの収束性は保証されず，得られる解の最適性も保証されない．この問題を解決する為に，筆者等は最尤推定による再構成を提案している [23]. 式 (14), (19)より
$d=F^{-1}U_{t}c$であるので，式 (20) より(21)$e=y-d=y-F^{-1}U_{t}c.$を得る．式 (21) を $p(e)$ に代入して対数をとる事により，対数尤度関数$l(t, c)=\log p^{(}p-F^{-1}U_{t}c)$.

およびを推定できる．
が定義される．この値を最大化する事により，( は未知の正実数， は単位行列) で密度関数 $p(e)$ が平均 , 共分散行列$\{t_{k}\}_{k=0}^{K-1}$$\sigma^{2}I$$0$$\{c_{k}\}_{k=0}^{K-1}$$I$$\sigma$ある場合，対数尤度関数は$l(t, c)=- \frac{||y-F^{-1}U_{t}c||^{2}}{2\sigma^{2}}-N\log(\sqrt{2\pi}\sigma)$となる．この最大化はノルム $||y-F^{-1}U_{t}c||^{2}$ の最小化と等価であり，更に

$F$

のユニ

タリ性より
(22)

$f_{0}(t, c)=||\hat{y}-U_{t}c||^{2}$の最小化に一致する．式 (22) は，最尤推定がフーリエ係数空間における最良近似ベクトルの探索問題に帰着される事を意味している．
でに関して 2 次形式であり，固定した に対する最適解は与えられる．ここで， は の Moore-Penrose 一般逆である．従って，式 (22) の式 (22) は$c$$t$$U_{t}^{\dagger}$$f_{0}(t, c)$を最小にする(23)を最小にする$t$$c=U_{t}^{\dagger}\hat{y}$$U_{t}$

および

$c$

は，

$f(t)=f_{0}(t, U_{t}^{\dagger}\hat{y})=||\hat{y}-U_{t}U_{t}^{\dagger}\hat{y}||^{2}$

および，それを用いて計算されるで求める事ができる．
の順序関れにより，探索パラメータ数を $2K$ から $K$ に削減できている．位置において制限することも可能である．
係を考慮すれば， の探索空間を$t$$c=U_{t}^{\dagger}\hat{y}$こ$\{t_{k}\}_{k=0}^{K-1}$$t$$[0, \tau]^{K}$

10
評価関数 $f(t)$ は凸ではないので，大域的な最適解を探索する事は非常に困難な問題である．本稿では，確率的なヒューリスティック手法である粒子群最適化法 (Particle

Swarm optimization, $PSO$ ) [24] と呼ばれる手法を用いる．各粒子は最適化パラメーを意味している． 個の粒子の位置と速度を探索領域における一様と粒子群全体の最適分布の乱数ベクトルで初期化する．個々の粒子の最適位置
位置をそれぞれ， とそれらの中での最適解で初期化する．終了条件が成立す
タである

$J$

$t$

$t_{j}$$i_{j}$$b_{j}^{(p)}$$b^{(s)}$$t_{1}$

るまで，各粒子の速度を$i_{j}:=wi_{j}+c_{1}r_{1}(b_{j}^{(p)}-t_{j})+c_{2}r_{2}(b^{(s)}-t_{j})$で，位置を$t_{j}:=t_{j}+i_{j}$で更新する．ここで，$c_{1},$$c_{2}$はあらかじめ定められた 1 に近い定数であり，$r_{1},$から 1 までの一様乱数である．もし
であれば
る．また，

$0$

$f(b_{j}^{(p)})<f(b^{(s)})$$f(t_{j})<f(b_{j}^{(p)})$$b^{(s)}$は$b_{j}^{(p)}$であれば$b_{j}^{(p)}$は$t_{j}$$r_{2}$

はで更新されで更新される．最終的に，$b^{(s)}$が探索$PSO$ は大域的かつ乱数的な手法であるので，局所解にとらわれる事が解を与える．
勾配法に比べて少ない．その一方で，計算量の多さが問題である．
計算機シミュレーションの結果を示す．周期 $\tau=1$ で $K=2$ 個のディラック関数からなる信号 $s(t)$ を考える．真値は $t_{0}=0.42,$ $t_{1}=0.52$ および $c_{0}=c_{1}=1.0$ である．
標本数は $N=11$ である．信号対雑音比

SNR$\underline{||d||^{2}}[dB]$$=10\log_{\iota 0_{\sigma^{2}N}}$で定まるガウス分布から雑音ベクトル を 1000 回生成した．得られたベクトル を用いて $t=(t_{0}, t_{1}),$ $c=(c_{0}, c_{1})$ を推定し，その結果の平均 2 乗誤差 (Mean Square Error, MSE)

が $0,5,$$\ldots,$$25$$[dB]$となる$e$$\sigma$$\hat{t},\hat{c}$$y$$MSE(t)=E||\hat{t}-t||^{2}, MSE(c)=E||\hat{c}-c||^{2}$を計算した．結果をそれぞれ図 1(a), (b) に示す．赤線，青線，黒線がそれぞれ，提案
手法，Cadzow 最小 2 乗法，最小 2 乗法による推定結果の平均 2 乗誤差を示す．横軸は SNR $[dB]$ であり，縦軸は MSE をデシベル表示した $10\log_{10}MSE$ である．SNR
が $0dB$ の場合の $MSE(c)$ を除いて，提案手法が最も小さな平均 2 乗誤差を与えていることがわかる．一方，SNR に関係なく 1 回の推定にかかる時間は，Cadzow 最
$+$

$+$

小 2 乗法，最小 2 乗法のどちらも .1 秒以下であるに関わらず，提案手法は約 6.8 秒
かかっている．このように提案手法の計算時間の短縮は重要な課題である．ただし，$PSO$ に比べて精度が若干劣るものの，Cadzow 最
高速化の研究は既に進んでおり，
小 2 乗法と同等の計算量で精度を改善できる手法が報告されている [25].
$O$

$+$

11

(a) MSE $(t)$( ) $MSE(c)$$b$図 1: 平均 2 乗誤差 (a) $MSE(t)$ , および (b) MSE(c).

赤線，青線，黒線それぞれが提案手法，Cadzow 最小 2 乗法，最小 2 乗法による結果 $[dB]$ を示す．
$+$

12

5 結論
本稿では不確定率有限信号の標本化理論を概説した．既知波形を平行移動したものの線形結合で表現される信号に対して，未知パラメータである平行移動量および結合係数の出現頻度を不確定率と定義し，その値が有限である信号を不確定率有限信号と定義した．そして，その信号の不確定率に近い頻度で得られた無雑音の sinc
標本値から，信号を完全再構成できることを示した．また，雑音が含まれる場合には，再構成に用いられるデータ行列の低ランク性と Toeplitz 構造に着目した Cadzowのアルゴリズムを紹介した．しかしこのアルゴリズムでは，低ランク行列集合の非凸性の故に最適性が保証されないことを指摘し，この問題を解決する為の最尤推定を提案した．この場合には，尤度関数が凸にはならないので，大域的確率探索手法
である粒子群最適化法を用いて再構成を行った．その結果，精度では従来手法を改善できるものの，計算時間が劣る事を示した．このため，高速化に関する研究も進められている事を紹介した．なお，本稿では画像直線エッジ抽出などの応用事例は紹介できていないので，関連文献 [26,27] を参照されたい．

謝辞
本研究は一部，JSPS 科研費 23500212 の助成を受けて行った．

参考文献
[1]

C. Shannon, “Communications in the presence of noise,” Proc. IRE, vol.37, pp.1021, 1949.

[2]

染谷勲，波形伝送，修教社，東京，1949.

[3]

小川英光，“標本化定理と染谷勲，” 電子情報通信学会誌，vol.89, no.8, pp.771-773,
2006.

[4] H. Ogawa, “Sampling theory and Isao Someya: $A$ historical note,” Sampling Theory

in Signal and Image Processing, vol.5, no.3, pp.247-256, 2006.
[5] A. Kohlenberg, “Exact interpolation of band-limited functions,” Joumal of Applied

Physics, vol.24, pp.1432-1436, 1953.

“Interferometric profiler for rough surfaces,” Applied Optics, vol.32,
no. 19, pp.3438-3441, 1993.

[6] P. Caber,

13
[7]

A. Hirabayashi, H. Ogawa, and K. Kitagawa, “Fast surface profiler by white-lightinterferometry by use of a new algorithm based on sampling theory,” Applied Optics,
vol.41, no.23, pp.4876-4883, 2002.

[8] “http://www.scn.tv/user/torayins/ $SP$ -500.html.”[9] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innova-

tion,” IEEE Transactions on Signal Processing, vol.50, no.6, pp.1417-1428, 2002.
[10]

C. Knapp and G. Carter, “The generalized correlation method for estimation of timedelay,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol.24, no.4,
pp.320-327, 1976.

[11] 小川英光，“一般標本化定理，” 電子情報通信学会論文誌，vol J71- $A$ , no.2, pp.163-

170, 1988.

[12] Y. Eldar and T. Dvorkind, “$A$ minimum squared-error framework for generalized

sampling,” IEEE Transactions on Signal Processing, vol.54, no.6, pp.2155-2167,

2006.
[13] M. Unser, “Splines:

perfect fit for signal and image processing,” IEEE Signal
Processing Magazine, vol.16, no.6, pp.22-38, 1999.
$A$[14] M. Unser and T. Blu, ”Cardinal exponential splines: Part I–Theory and filtering

algorithms,” IEEE Transactions on Signal Processing, vol.53, no.4, pp.1425-1438,

2005.
[15] P.L. Dragotti, M. Vetterli, and T. Blu, ”Sampling moments and reconstmcting sig-

nals of finite rate of innovation: Shannon meets Strang-Fix,” IEEE Transactions on
Signal Processing, vol.55, no.5, pp.1741-1757, 2007.
[16] R. Tur, Y Eldar, and Z. Friedman, ”Innovation rate sampling of pulse streams withapplication to ultrasound imaging,” IEEE Transactions on Signal Processing, vol.59,
no.4, pp. 1827-1842, 2011.
[17] R. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE

Transactions on Antennas Propagation, vol.34, no.3, pp.276-280, 1986.
[18] R. Roy, A. Paulraj, and T. Kailath, “ESPRIT–a subspace rotation approach to esti-

mation of parameters of cisoids in noise,” IEEE Transactions on Acoustic, Speech,
Signal Processing, vol.34, no.5, pp.1340-1342, 1986.

14[19] Y. Hua and T.K. Sakar, “On SVD for estimating generalized eigenvalues of singu-

lar matrix pencil in noise,” IEEE Transactions on Signal Processing, vol.39, no.4,

pp.892-900, 1991.
[20] J. Cadzow, ”Signal enhancement–a composite property mapping algorithm,” IEEE

Transactions on Acoustic, Speech, and Signal Processing, vol.36, no. 1, pp.49-62,
1988.
[21] T. Blu, P.L. Dragotti, M. Vetterli, P. Marziliano, and L. Coulot, “Sparse samplingof signal innovations,” IEEE Signal Processing Magazine, vol.25, no.2, pp.31 $AO,$

2008.
[22] I. Markovsky, “Structured low-rank approximation and its applications,” Automatica

J. IFAC, vol.44, no.4, pp.891-909, 2008.
[23] A. Hirabayashi, “Sampling and reconstruction of periodic piecewise polynomials

using sinc kemel,” IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, vol.E95- , no. 1, pp.322-329, 2012.
$A$[24] J. Kennedy and R. Eberhart, “Particle swarm optimization,” Proceedings of the 1995

IEEE Intemational Conference on Neural Networks, pp.1942-1948, $nov/dec$ 1995.
[25] L. Condat and A. Hirabayashi, ”Cadzow denoising upgraded:$A$

new projection

method for the recovery of Dirac pulses from noisy linear measurements,” HAL
Archives, http: $\int/hal.archives-$ ouvertes. $fr/docs/OO/75\int 97/16\int PDF\int Condat$ -SLRA.pdf.
[26] L. Baboulaz and P. Dragotti, “Exact feature extraction using finite rate of innovation

principles with an application to image super-resolution,” Image Processing, IEEE
Transactions on, vol.18, no.2, pp.281-298, 2009.
[27] A. Hirabayashi and P.L. Dragotti, “Line-edge extraction based on $E$ -spline acquisi-tion model and a fast optimization algorithm,” Proceedings of Intemational Conference on Image Processing (ICIP 2012), Orlando, pp.89-92, 2012.