プロシーディングス

2012 RIMS 共同研究
数学・数理科学と諸科学・産業との連携研究ワークショップ
ウェーブレット解析とサンプリング理論
2012 RIMS Joint Research & MEXT Workshop
Wavelet analysis and sampling theory
主催：京都大学数理解析研究所，文部科学省
場所：京都大学数理解析研究所
日程：平成 24 年 9 月 10 日（月）12:50 - 17:00
平成 24 年 9 月 11 日（火） 9:30 - 15:30
2012 RIMS 共同研究
「ウェーブレット解析とサンプリング理論」
研究代表者：芦野 隆一（大阪教育大学）
日時：2012 年 9 月 10 日（月） ∼ 2012 年 9 月 11 日（火）
会場：京都大学数理解析研究所 111 号室
〒 606-8502 京都市左京区北白川追分町
この RIMS 共同研究は，文部科学省と京都大学数理解析研究所が共催する
平成 24 年度 文部科学省
数学・数理科学と諸科学・産業との連携研究ワークショップ
です．連携研究ワークショップに関しては，
http://www.mext.go.jp/a_menu/math/index.htm
をご覧下さい．
ワークショップの開催趣旨
広い意味でウェーブレット解析やサンプリング理論によって解決で
きるかもしれないと期待できるトピックスに関して，講演者の方々
に理論と工学的応用の現状，さらに解決すべき問題を解説していた
だき，その問題提起を端緒として参加者がディスカッションする形
で，ウェーブレット解析やサンプリング理論が実際にどのように応
用されているかをより深く理解することによって，新しい理論と産
業への応用の道が開かれることを目指す．
1
2
プログラム
9 月 10 日（月）
12:50 – 13:00 開会の挨拶
13:00 – 14:30 田中 利幸（京都大学大学院 情報学研究科）Toshiyuki TANAKA
[基調講演] 圧縮センシングの基礎数理
x, y をベクトル，A を行列として，A, y を既知として連立一次
方程式 y = Ax を解いて x を求める問題においては，x に対して
スパースさなどの仮定をおくと，変数よりも方程式のほうが数が
少ない状況でも解が求められる場合がある．圧縮センシングにお
ける基礎数理的に重要な課題のひとつに，x がどのくらいスパー
スであればどのくらい少ない数の方程式から解が求められるかを
調べる問題がある．この問題に関して，最適化，ランダム行列，
積分幾何学，大自由度統計などの多様な分野にまたがる横断的な
視点から，基礎的な結果ならびに最近の進展について概説する．
14:45 – 15:45 平林 晃（山口大学大学院 医学系研究科 応用分子生命科学系専攻）Akira
HIRABAYASHI
連続信号に対するスパースサンプリング
連続信号に対する標本化と言えばこれまで，ナイキスト間隔で行
われることが常識的となってきた．しかし，レーダーやソナー，
エコーなどに現れる信号は広帯域でありながら，特定波形を平
行移動したものの線形結合になっていることに着目すれば，ナイ
キスト間隔より遥かに広いスパースサンプリングが可能になる．
また，区間ごとに多項式や指数関数で定義された信号は帯域制限
を受けてすらいないが，同様のスパースサンプリングが可能であ
る．本講演では，これらの信号を包含する不確定率有限信号とよ
ばれる信号クラスを紹介し，この信号クラスに対する標本化と再
構成の理論を概説する．また，エコー信号や心電図波形など実信
号に対する応用事例を紹介する．
16:00 – 17:00 山田 道夫（京都大学 数理解析研究所）Michio YAMADA
時間周波数解析手法と Wigner 分布
理工学における時間周波数解析手法は，窓付きフーリエ変換，Wigner
分布，ウェーブレット解析，のいずれかを基礎とするものが多い．
ここでは Wigner 分布についての review を試みる．Wigner 分布
と他の時間周波数解析手法の関係，Cohen 分布への拡張，伏見関
数などについて述べる．
懇親会
3
9 月 11 日（火）
9:30 – 10:30 黒石 裕樹 （国土交通省 国土地理院）Yuki KUROISHI
半離散化ウェーブレット変換を応用した測地データの局在化信号分離の
試み
半離散化ウェーブレット変換に基づき，一次元関数及び二次元関
数の局在化信号展開と再構築を行うフレームを構成する．これ
らを測地学分野の時系列（一次元関数）と地理的分布データ（二
次元関数）に適用し，局在化信号に展開する．さらに，異なる二
地点で得られた時系列のウェーブレット展開から，時系列間のコ
ヒーレンシー分析を行う．それらの結果から，半離散化ウェーブ
レット・フレーム表現に基づく入力データのフィルター効果を調
べる．
10:45 – 11:45 池田 和司（奈良先端科学技術大学院大学 情報科学研究科）Kazushi IKEDA
脳情報解析およびバイオインフォマティクスにおけるスパース信号処理
圧縮センシングなどスパース信号処理の応用は様々であるが，本
発表では脳情報解析およびバイオインフォマティクスに応用した
例について紹介する．脳情報解析では脳波を用い，スパースロジ
スティックレグレッションでひらめき判別器を構成することで，
ひらめきに関連する脳活動の時空間パターンを抽出した．バイオ
インフォマティクスでは compound を用いた実験結果を用い，圧
縮センシングにより神経伸長に関連する kinase を特定した．
13:15 – 14:15 吉野 邦生（東京都市大学 知識工学部 自然科学科）Kunio YOSHINO
佐藤超関数による Shannon–染谷の標本化定理の拡張と Ramanujan の
積分公式
佐藤超関数を用いる事により Shannon–染谷の標本化定理を拡張
する．特に，信号（関数）に対する二乗可積分の仮定をはずす事
が出来る．更に複素解析的な手法を用いて Ramanujan の積分公
式との関係を明らかにする．
14:30 – 15:30 岡田 正已（首都大学東京 数理）Masami OKADA
サンプリング値の補間による未知関数の再構成
規則格子でのサンプリング観測値から，もとの信号関数をできる
だけ忠実に再構成するための数学的基礎は古くから知られてい
る．ここでは，比較的に最近の文献に従って，正定型関数を用い
た不規則格子でのサンプリング補間について大まかに解説し，今
後の発展の可能性を探る．
圧縮センシングの基礎数理
田中 利幸
京都大学大学院情報学研究科
概要.
圧縮センシングの数理について，基礎的な事項を整理する．疎性概念の歴史と現
代科学における有効性についても議論する．
Basic mathematical framework for compressed sensing
Toshiyuki Tanaka
Graduate School of Informatics, Kyoto University
Abstract.
We review a basic mathematical framework for compressed sensing. The
notion of sparsity, as well as its historical accounts and its efficiency in modern sciences,
is discussed.
1. 圧縮センシング
圧縮センシングの基本的な問題設定は，未知のベクトル x0 に対して線形観測を行い，
その結果 y にもとづいて未知ベクトル x0 を推定することである．線形観測結果 y を得る
過程が
y = Ax0 + n
(1.1)
と表されるものとする．ここで， A は線形観測を定義する行列であり， n は観測における
ノイズを表す．線形観測結果 y と行列 A とは既知であり，これらから未知ベクトル x0 を
推定することがここでの問題である．
ノイズが x0 や A に対して独立である場合には，x0 を推定する手段として最小二乗推定
(1.2)
x̂(LMS) = arg min y − Ax22
x
がよく用いられる．ノイズベクトル n の各要素が独立に平均が 0, 分散が同一の正規分布
に従うと仮定するならば，最小二乗推定 (1.2) は y を既知としたときの x0 の最尤推定と一
致する．
未知ベクトル x0 の次元を N, 線形観測結果 y の次元を M とする． M < N の場合，より
正確には rank A < N の場合には，最小二乗推定 (1.2) による結果は一意とならない．この
場合には，A の零空間 ker A = {x | Ax = 0} が次元 (N − rank A) の線形空間となり，最小二
乗推定の結果に ker A の任意の元を加えたものもやはり最小二乗推定の解になっているか
らである．このような不定性の問題を回避するためによく用いられてきたのは，正則化に
1
もとづく方法である．古典的なティホノフ正則化のアプローチでは，
x̂(TR) = arg min y − Ax22 + λx22 ,
(1.3)
x
λ≥0
によって x0 の推定結果を得る．このアプローチは，統計学の分野ではリッジ回帰とよば
れている．また， x0 の事前確率として各要素が独立に平均が 0, 分散が同一の正規分布に
従い，さらにノイズベクトルの各要素も独立に平均が 0, 分散が同一の正規分布に従うと
仮定するならば， x0 の最大事後確率 (MAP) 推定はティホノフ正則化 (1.3) においてパラ
メータ λ を適切に定めたものに等価である．ティホノフ正則化 (1.3) の目的関数は x につ
いて 2 次式であるから，その解は代数的に求められ，結果は
x̂(TR) = (AT A + λI)−1 AT y
(1.4)
となる．
正則化のアプローチは，不定性の問題を回避するために未知ベクトル x0 に対してある
種の事前知識を仮定しているものとみなすことができる．ティホノフ正則化では，
「 x0 の
2-ノルムは小さい」という事前知識を仮定しているものと解釈できる．不定性を回避する
ために導入する仮定としては，他にも様々な形のものを考えることができる．なかでも近
年注目されているのが，疎性に関する事前知識である．ベクトル x の非零要素の個数を
x0 によってあらわし，0-「ノルム」とよぶことにする*1 ．ベクトルが疎であることは，
その 0-「ノルム」が小さいことと対応している． x0 は疎である，すなわち，「x0 0 は小
さい」という形の事前知識を仮定するならば，x0 の推定に対する正則化のアプローチとし
て，ティホノフ正則化にならって 0-「ノルム」正則化
(1.5)
x̂(0NR) = arg min y − Ax22 + λx0
x
を考えることができる．
0-「ノルム」正則化は，目的関数が非凸であるという実用上の大きな欠点がある．0-「ノ
ルム」正則化を解くことは x のどの要素が非零であるかに関する組み合わせ最適化問題を
解くことに相当するため，一般には計算複雑度の点から難しい問題であると考えられてい
る．0-「ノルム」の代わりにその凸緩和として 1-ノルムを使う 1-ノルム正則化
(1.6)
x̂(1NR) = arg min y − Ax22 + λx1
x
では，目的関数も全体として凸関数になるため，上記の欠点を回避できる． 1-ノルム正則
化も解が疎性をもつため，疎性に関する事前知識を仮定して未知ベクトル x0 を推定する
問題に対する実用的な取り扱いを可能にする定式化として広く使われている．未知ベクト
ル x0 が疎である場合には，x0 の次元 N よりも線形観測結果の個数 M が小さい場合でも，
1-ノルム正則化によってよい精度の疎な推定結果が実用的な程度の計算複雑度で求められ
*1
0-「ノルム」は数学的な意味のノルムではない．
2
る．圧縮センシングに関する研究の主要な関心は，1-ノルム正則化の性質を議論したり，
効率的な解法を考案したりすることに向けられている．
信号処理や画像処理などの応用事例では，未知ベクトル x0 として信号や画像そのもの
を考えると，それらは一般には疎性をもたないため，上述の枠組みをそのまま適用するこ
とはできない．一方で，我々が処理の対象とするような信号や画像は，ウェーブレット変
換などの変換を施すと疎性を示すことが多い．このような対象を取り扱うには，未知ベク
トル x0 はある行列 Ψ と疎ベクトル z とによって
x0 = Ψz0
(1.7)
と表されるものと仮定すればよい．この場合には，1-ノルム正則化は
(1.8)
ẑ(1NR) = arg min y − AΨz22 + λz1
z
という形で疎ベクトル z0 を推定する問題として定式化できる．行列 Ψ が正則であれば，
Ψ の列ベクトルは RN の基底を与えており， z0 はその基底によって x0 を与える線形結合
の係数である．行列 Ψ として「横長」の長方形行列を使う場合もある．この場合には，Ψ
の列ベクトルの個数は N より大きいからこれらは RN の基底ではない（「過完全基底」あ
るいは「フレーム」とよばれることがある）が， z0 を係数とする Ψ の列ベクトルの線形
結合で x0 があらわされるという点については同様である．いずれの場合も， A の代わり
に AΨ が線形観測を定義しているものとみなせば，変換を考えない場合と同形の問題とし
て取り扱うことができる．
2. 疎性について
最小二乗推定 (1.2) を提案したのはガウスであるが，ノイズや観測誤差が存在する条件
下での推定の問題はガウス以前から議論されており，ガリレオやラプラスは 1-ノルム最
小化
(2.1)
x̂ = arg min y − Ax1
x
に相当する方法を考えていたとされている． 1-ノルム最小化によって推定を行うと，残差
ベクトル y − A x̂ が疎になる傾向があることは早くから認識されていたようで，ガウスは
1-ノルム最小化のこの性質を嫌い，代数的に簡潔に解を与えることができる最小二乗推定
を提案するに至ったようである．
2-ノルムは座標系の回転に関して不変である．また，考える対象によっては物理的なエ
ネルギーとの対応づけが可能である．この場合，残差の 2-ノルムを 2 乗したもの y− Ax22
はノイズのエネルギーに， x の 2-ノルムの 2 乗 x22 は原信号のエネルギーにそれぞれ対
応しているとみなすことができる．このような対応関係のもとでは， 2-ノルム正則化は物
理的には最小エネルギー原理と対応づけられる．
3
一方で，ベクトルの疎性という性質は座標系に依存する概念であり， 2-ノルムの場合と
は異なり，0-「ノルム」や 1-ノルムは座標系の回転に関する不変性をもたない．このこと
は，圧縮センシングの枠組みが物理学の基本的な指導原理のひとつである「物理法則はそ
れを表現する座標系の選び方によらない」という要請と適合していないことを示してい
る．しかし，物理的な領域を離れて広く情報処理の分野に問題を求めるならば，ベクトル
の各要素は往々にして固有の意味を有しており，したがって座標系の回転に意味をもたせ
ることができない場合がほとんどである．さらに，線形回帰分析においてはデータへの適
合度を高めることはもちろん重要であるが，少数の変数によってデータを適切に説明する
ことがもし可能であれば，そのような結果を引き出すことも劣らず重要であり，この目
的を達成するために 1-ノルム正則化などの疎性に注目した定式化を適用できる． LASSO
は，このような定式化の一例である．これらのアプローチはいわば「やや主観的なデータ
解析」ということもできようが，数多くの変数の中からデータを適切に説明する少数の変
数を抽出することは，データの背後にある本質的な構造を我々が理解し活用することにも
つながるであろう．第一原理的な自然法則から種々の結果を演繹できる物理学とは異な
り，実験データからそれに対する簡潔な説明を帰納的に抽出することで現象の理解を深め
ていく科学の諸分野において，疎表現はきわめて有効な概念であると思われるが，このこ
とは方法論として明確に意識され定着しているわけではないようである．本ワークショッ
プでは，諸科学における方法論としての疎表現概念の重要性，有効性について議論したい．
田中 利幸 (京都大学大学院情報学研究科システム科学専攻)
〒606-8501 京都府京都市左京区吉田本町 36-1
E-mail: [email protected]
4
連続信号に対するスパースサンプリング
∗
平林 晃 ∗
山口大学 大学院 医学系研究科 応用分子生命科学系専攻
概要.
連続信号に対する標本化と言えばこれまで，ナイキスト間隔で行われることが常
識的となってきた．しかし，レーダーやソナー，エコーなどに現れる信号は広帯域であり
ながら，特定波形を平行移動したものの線形結合になっていることに着目すれば，ナイキ
スト間隔より遥かに広いスパースサンプリングが可能になる．また，区間ごとに多項式や
指数関数で定義された信号は帯域性制限を受けてすらいないが，同様のスパースサンプリ
ングが可能である．本講演では，これらの信号を包含する不確定率有限信号とよばれる信
号クラスを紹介し，このクラスに対する標本化と再構成の理論を概説する．また，画像特
徴量抽出や心電図波形の圧縮など，本理論の応用事例を紹介する．更に，本理論に基づい
て予想される今後の信号処理におけるパラダイムシフトを展望し，ご意見を伺う．
Sparse sampling of continuous signals with finite rate of
innovation
Akira Hirabayashi∗
∗
Yamaguchi University
Abstract.
We present a survey of sampling theory for signals with a finite rate of
innovations. In contrast to the so-called compressed sensing, this theory focuses on sparse
sampling of continuous signals. Signals in radar, echo, or sonar are expressed by linear
combinations of shifted versions of a know waveform. Those signals can be represented
using a finite number of unknown parameters per unit. Piecewise polynomials or sinusoids
are also expressed in a same way. A signal with finite rate of innovation is a class that
contains signals having such parametric representations. We can uniformly sample and
perfectly reconstruct such signals with much wider interval than the Nyquist, which is
sometimes infinite. We quickly review this theory in noiseless and noisy cases and show
applications on image feature extraction and compression of biomedical signals.
1. 不確定率有限信号
簡単の為に周期的な場合を考える．周期 τ の K(<∞) 個のディラック列 s(t) は，区間
[0, τ) の信号
s0 (t) =
K−1
∑
k=0
1
ck δ(t − tk )
を用いて (0 ≤ t0 < t1 < . . . < tK−1 < τ)，
s(t) =
∞
∑
s0 (t − k0 τ)
k0 =−∞
と表すことができる．デイラック列ではない一般的なパルス列は， s(t) と単一パルスとの
畳込みで表現できるので，サンプリング・再構成も以下と同様に実行できる．
K−1
K−1
この信号の１周期には，位置 {tk }k=0
と対応する係数 {ck }k=0
の合計 2K 個の自由パラ
メータが存在する．そこで不確定率 ρ を 2K/τ と定義する．K が有限であるので，ρ も有
限である．このような信号を不確定率有限信号（Signals with finite rate of innovation）と
呼ぶ [1]．
信号 s(t) は周期的なのでフーリエ級数で表現され，その係数は
(1.1)
1
d̂ p =
τ
∫
τ
1∑
ck ukp
τ k=0
K−1
s(t)e−i2pπt/τ dt =
0
となる．ここで uk = e−i2πtk /τ である．このようにフーリエ係数がべき級数 uk の線形結合
p
になっていることが重要である．
2. 標本化と再構成
N−1
雑音を含まない標本値 {dn }n=0
は，
(2.1)
∫
dn = hs, ψn i =
∞
−∞
f (t)ψ(t − nT )dt
で与えられる．ここで T = τ/N であり，標本化核は ψ(t) = Bsinc(Bt) である．更に
B ≥ ρ = 2K/τ であり
{
sinc(t) =
sin(πt)/(πt) (t , 0)
1
(t = 0)
である．
N−1
K−1
K−1
本論文で論じる課題は，標本値 {dn }n=0
から 2K 個の未知数 {tk }k=0
と {ck }k=0
を求め
ることである．この問題はアニヒレーティングフィルタを用いることによって解決でき
る [1]．概略を述べる． Bτ/2 を超えない最大整数を P で表す： P = bBτ/2c．このとき，
ポアソン総和式
∞
∑
P
1 ∑ −i2pπt/τ
B
sinc{B(t + k τ)} =
e
τ p=−P
k0 =−∞
0
より，
dn =
P
∑
d̂ p ei2pnπ/N
p=−P
2
が成立する．この式は，フーリエ係数 d̂ p が逆離散フーリエ変換によってシンク核による
サンプル値 dn に対応付けられることを意味している．従って，順変換によって dn から d̂ p
を求めることができる：
(2.2)
N−1
1 ∑
d̂ p =
dn e−i2pnπ/N
N n=0
系列 {d̂ p }Pp=−P は z 変換が
A(z) =
K
∑
ak z−k =
k=0
K−1
∏
(1 − uk z−1 )
k=0
であるフィルタ [1, a1 , . . . , aK ] によって零変換される．ここで，a0 = 1 である．零変換関
係は
(2.3)
d̂ p + a1 d̂ p−1 + . . . + aK d̂ p−K = 0
であり，ここで p = K − P, . . . , P である．これらの式を連立して解くことによりフィルタ
係数が求まる．そして，この係数からなる多項式を因数分解することにより uk が求まり，
更に tk が求まる．K 個の式 (2.3) では，2K 個の d̂ p が必要である．ただし，実数 dn を表
現するためのフーリエ級数の対称性のため，係数 d̂ p は奇数個となる．従って，2K + 1 個
の dn を必要とする．tk または uk が求まれば，式 (1.1) を用いて ck を求めることができ
る．これにより，雑音を含まない場合の再構成が完了する．雑音が含まれる場合にはデー
タ行列の特異値分解や，Cadzow アルゴリズムがよく用いられるが [2]，最尤推定を用いれ
ばより高精度の再構成が可能になる [3]．本理論を２次元信号へ拡張すれば，画像におけ
る直線エッジを抽出できる [4, 5]．
謝辞 本研究は JSPS 科研費 23500212 の助成を受けたものである．
参考文献
[1] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,”
IEEE Transactions on Signal Processing, vol.50, no.6, pp.1417–1428, June 2002.
[2] T. Blu, P.L. Dragotti, M. Vetterli, P. Marziliano, and L. Coulot, “Sparse sampling of
signal innovations,” IEEE Signal Processing Magazine, vol.25, no.2, pp.31–40, March
2008.
[3] A. Hirabayashi, T. Iwami, S. Maeda, and Yosuke Hironaga, “Reconstruction of the
sequence of Diracs from noisy samples via maximum likelihood estimation,” Proceedings of the 37th International Conference on Acoustics, Speech, and Signal Processing
(ICASSP 2012), pp. 3805–3808, Kyoto, March 2012.
3
[4] L. Baboulaz and P. Dragotti, “Exact feature extraction using finite rate of innovation
principles with an application to image super-resolution,” IEEE Transactions on Image
Processing, vol.18, no.2, pp.281–298, Feb. 2009.
[5] A. Hirabayashi and P.L. Dragotti, “ Line-edge extraction based on E-spline acquisition
model and a fast optimization algorithm,” Proceedings of International Conference on
Image Processing (ICIP 2012), Orlando, Oct. 2012.
平林 晃 (山口大学 大学院 医学系研究科 応用分子生命科学系専攻)
〒755-8611 山口県宇部市常盤台 2-16-1
E-mail: [email protected]
4
時間周波数解析手法と Wigner-Ville 分布
山田 道夫
京都大学数理解析研究所
概要.
信号の時間周波数解析には便宜的なものも含めさまざまな解析手法が考案されて
きた．これらの時間周波数解析の手法は大きく分けて，局所フーリエ変換，Wigner-Ville
分布，ウェーブレット変換，それぞれの系統があるが，これらは別々のものではなく相互
に関係する概念である．窓付きフーリエ変換とウェーブレット変換は信号の線形変換に基
礎を置くのに対し，Wigner-Ville 分布は信号の二次形式に基礎を置いている．ここでは
Wigner-Ville 分布をめぐる時間周波数解析について若干の review を試み，今後の課題
について参加者方々のご意見を伺いたい．
Time-frequency analysis and Wigner-Ville distribution
Michio YAMADA
RIMS, Kyoto University
Abstract.
Analysis methods for time-frequency structure of a signal are classified into
three groups; short-time Fourier transforms, Wigner-Ville distribution and wavelet transforms. The first and the third ones are based on linear transforms, and the second one is
on a quadratic form of the signal. These methods are not independent of but related to
each other. A brief review is given on the Wigner-Ville distribution and related concepts.
1. 窓付きフーリエ変換
関数（信号） f (x) ∈ L2 (R) の周波数構造を表すものとして，理工学ではフーリエ変換
fˆ(ω) =
∞
e−iωx f (x) dx,
−∞
がよく用いられるが，fˆ(ω) は時間変数を含まないため，信号の周波数分布が時間的に変化
する場合にはフーリエ変換は必ずしも便利ではない．そこで信号の時間局所的な周波数構
造を調べるため，窓関数 w(x) ∈ L2 (R) によって x = b 付近に制限した信号 w(x − b)f (x)
のフーリエ変換
Sf (b, ω) =
∞
−∞
e−iωx w(x − b)f (x) dx,
がしばしば用いられる (上線は複素共役)．これを窓付きフーリエ変換 (windowed Fourier
transform) といい，|Sf (b, ω)|2 はスペクトログラム（spectrogram）とよばれる [1]．
1
w = 1 を満たす窓関数に対し，再生核
∞
K(b0 , b, ω0 , ω) =
w(t − b)w(t − b0 )e−i(ω0 −ω)t dt,
−∞
を
K(b0 , b, ω0 , ω) = e−i(ω0 −ω)(b+b0 )/2 Aw (b0 − b, ω0 − ω)
∞
b̃
b̃
w(τ + )w(τ − )e−iω̃τ dτ,
Aw (b̃, ω̃) =
2
2
−∞
と書くとき，Aw は窓関数 w の曖昧さ関数（ambiguity function）であり窓付きフーリエ
変換の時刻と周波数の精度を表す．
2. Wigner-Ville 分布と Cohen 分布
気体の量子統計力学の研究において 1932 年 Wigner は位置と運動量の結合分布を提
案し，1948 年 Ville はこれを信号解析の分野に導入した．Wigner(-Ville) 分布とは信号
f (x) ∈ L2 (R) に対する次の 2 次形式をいう [2]．
∞
x
x
Wf (b, ω) =
f (b + )f (b + )e−iωx dx.
2
2
−∞
Wigner 分布は，信号が複素数のときにも実数値となり，信号が時間軸あるいは周波数軸
に平行移動すれば同様に平行移動する（推移不変性）．Wigner 分布は時間と周波数につ
∞
∞
いての周辺条件 −∞ Wf (b, ω) db = |fˆ(ω)|2 ,， −∞ Wf (b, ω) dω = |f (b)|2 ,, を満たす．ま
ˆ のサポートがコンパクトな場合は Wf のサポートの時間軸 [周波数軸] への射影は
た f [f]
f [fˆ] のサポートに含まれるなど時間周波数面における分布として自然な性質を備えてい
る．さらに f, g ∈ L2 (R) に対し次の等長性（Moyal の関係）が成り立つ．
2
∞ ∞
1
f (x)g(x) dx =
Wf (b, ω)Wg (b, ω) db dω
2π −∞ −∞
−∞
∞
Wigner 分布はこれらの良い性質と共に，2 次形式であることに起因する欠点も伴っ
ている．その一つは，信号の和 f (x) = f1 (x) + f2 (x) の Wigner 分布には Wf =
Wf1 + Wf2 + I(f1 , f2 ) のように干渉項 I(f1 , f2 ) が生じることである．この項がしばし
ば物理的解釈が困難な正負の振動を伴う．またより深刻なのは Wf (b, ω) の正値性が保証
されないことである．これは原理的な問題で，一般に周辺条件 (2)(2) を満たす Hermite
2 次形式の非負値分布 Wf (b, ω) は存在しない (Wigner) ことが知られている．なお正値
性を回復するため，Wigner 関数を Gauss 関数でならすことも行われている（伏見関数，
Q 関数）．
2
Moyal の関係式 (2) で g(x) = w(x − b)eiωx とおくと，スペクトログラムが信号と窓
（連続ウェーブレット変換についても同様の表現
関数の Wigner 分布によって表される．
が可能である．
）
1
|Sf (b, ω)| =
2π
2
∞
−∞
∞
−∞
Wf (b , ω )Ww (b − b, ω − ω) db dω 一般に，曖昧さ関数と Wigner 変換は互いにフーリエ変換の関係にあり，スペクトログラ
ムのフーリエ変換（特性関数）
1
Ŝf (b̃, ω̃) =
2π
∞
−∞
∞
−∞
eib̃ω−iω̃b Sf (b, ω) db dω
は信号および窓関数の曖昧さ関数の積 Ŝf (b̃, ω̃) = Af (b̃, ω̃)Aw (b̃, ω̃) として表される．そ
こでこの式で窓関数の曖昧さ関数を一般の二変数関数（核関数）φ(b̃, ω̃) に置き換えて得
られる分布
1
C(b, ω) =
2π
∞
−∞
∞
−∞
e−ib̃ω+iω̃b Af (b̃, ω̃)φ(b̃, ω̃) db̃ dω̃
は 1966 年に Cohen によって提案され Cohen 分布とよばれている [2]．Cohen 分布の性
格は核関数によって決定される．特に φ(b̃, ω̃) = 1 の場合は Wigner 分布となるが，他の
適当な核関数を選ぶことによってさまざまの時間周波数分布を得ることができる．
短時間フーリエ変換やウェーブレット変換は微分作用素との相性があまり良くないた
め，データ解析に多用される反面，微分方程式に用いられる機会は少ない．しかし Wigner
分布はもともと量子統計力学的な背景のもとに考案されたもので，微分方程式の解析にも
使用されることがある．その一例は，次の Davey-Stewartson 方程式
(2.1)
(2.2)
i
∂ 2A
∂A
∂ 2A
+λ
+
μ
= ν|A|2 A + ν1 AQ,
2
2
∂τ
∂X
∂Y
∂2Q
∂ 2Q
∂ 2 |A|2
+
μ
=
κ
λ1
1
1
∂X 2
∂Y 2
∂Y 2
の解のアンサンブルに対し，準正規近似の下にエネルギースペクトルの定常性を議論した
ものである．多くの時間周波数解析手法が観測実験データや数値解のデータ解析に留まる
中，基礎方程式の理論解析と直接関わっている点で Wigner 分布はデータ解析手法として
独特の位置を占めている．
参考文献
[1] S.Mallat, a wavelet tour of signal processing, The Sparse Way, Academic Press,
2009.
[2] L.Cohen, Time-frequency analysis, Prentice Hall, 1995 (邦訳：時間-周波数解析，吉
．
川昭，佐藤俊輔訳，朝倉書店，1998）
3
半離散化ウェーブレット変換を応用した測地データ
の局在化信号分離の試み
∗
黒石 裕樹 ∗
国土交通省国土地理院
概要.
半離散化ウェーブレット変換に基づき，一次元関数及び二次元関数の局在化信号
展開と再構築を行うフレームを構成する．これらを測地学分野の時系列（一次元関数）と
地理的分布データ（二次元関数）に適用し，局在化信号に展開する．さらに，異なる二地
点で得られた時系列のウェーブレット展開から，時系列間のコヒーレンシー分析を行う．
それらの結果から，半離散化ウェーブレット・フレーム表現に基づく入力データのフィル
ター効果を調べる．
Application of semi-discrete wavelet transform to segregation
of localized signals in geodetic data
Yuki Kuroishi∗
∗
Geospatial Information Authority of Japan,
Ministry of Land, Infrastructure, Transport and Tourism
Abstract.
Based on semi-discrete wavelet transforms, two kinds of snug frame are constructed to expand a one-dimensional or two-dimensional function into a group of wavelets
and to reconstruct the function from them. Such frame representations are applied to
decompose into localized signals time series data and geographically distributed data in
geodesy. Then a pair of the wavelet transforms of time series data at two separate locations
are used to analyze the coherency of signals localized in time-scale domain between the
two. From the results obtained we investigate the effectiveness of semi-discrete wavelet
transforms in filtering data to isolate localized features.
1. はじめに
測地学分野においては，全地球的あるいは地域的な物理量分布を異種観測から決定する
問題や，多様な地球物理的諸現象を含む物理量の連続的観測について，地理的に近接する
地点での共通性変化を分離，抽出する問題などを扱うことがある．ここでは，前者に属す
るものとして日本列島周辺の海域重力場のモデル化，後者に属するものとして潮位観測に
おけるコヒーレンシー分析を対象として考える．
海域重力場に関しては，人工衛星アルチメトリー（海面高度計）による全球海域の重力
場推定モデル（以下，高度計モデル）と海（船）上重力計による計測が利用可能である．
それぞれの手法には固有の誤差要因があり，地理的位置により異なるスケール分布を持
1
つ，局在的な系統的誤差を含んでいる．そのため，これら二種類の観測による測定値の較
差分布について，局在化信号分解を行うことができれば，誤差の局在的特性を把握し，誤
差低減を実現できると考えられる．
一方，潮位の時間変化については，潮汐，気圧応答といった短周期変動，海水の温度・
塩分濃度の変化によるステリックな変化，海流や風などによる変動や氷河の融解による海
水総量の変化，経年的熱膨張などに伴った長期的変化をしている．潮位観測から海面上昇
などの中長期的な変化を知るためには，観測点毎に推定された潮汐変動・気圧応答を除去
された潮位時系列について，隣接する点間でコヒーレントな変化を示すと思われる，ステ
リックな変化と海流などによる変化を推定，除去する必要がある．
本稿では，Kuroishi and Keller(2005) と黒石 (2012) を再編し，これらの問題について，
半離散化ウェーブレット変換に基づくフレームを用いた信号の局在化展開とそれらをコ
ヒーレンシー分析に応用する試みについて紹介する．
2. 半離散化ウェーブレット変換に基づくフレーム
2.1
Morlet ウェーブレットを用いた一次元半離散化ウェーブレットの
フレーム
時系列データについて，progressive な関数であって周波数領域における関数形（特性）
が複雑でなく，フレームをなす際の冗長性が小さなものとして，Morlet ウェーブレットを
母関数とする一次元半離散化ウェーブレット変換を適用する．Morlet ウェーブレットは，
時間領域と周波数領域において，それぞれ，次式で表される：
Ψ(t) = (1/π1/4 ) exp(−t2 /2){exp(−iω0 t) − exp(−ω20 /2)}
√
Ψ̂(ω) = ( 2π1/4 ) exp{−(ω − ω0 )2 /2}.
ただし，ˆ は周波数領域の関数であることを示し，t と ω はそれぞれ時刻と周波数の変数，
i は虚数単位，ω0 は中心周波数であり，以下の値をとる場合，上式の中括弧内の第二項は
無視できる：
p
ω0 = π 2/ ln 2 5.3364.
この値は，周波数領域における関数値において，最大値と第一サイドローブのピーク値の
振幅比を 2 : 1 とするものである．
離散化ウェーブレット変換においてスケール・パラメータを s = 2 j/4 （ j は 整数；1 オ
クターブにつき 4 個のスケールというマルチ・ボイスをなす），時刻パラメータを b = n
（ n は自然数；データ間隔を単位とする）にとるとき，このウェーブレット関数を用いた
フレームの冗長度の下限 A と上限 B は，それぞれ，6.918 と 6.923 であり，復元における
誤差は，冗長度の上限と下限の比より 0.04 % である（Daubechies，1992）
．Kuroishi and
Keller（2005）は，スケール・パラメータについてこの離散化を導入し，時刻パラメータ
2
については等間隔サンプリングによるデジタル・データの間隔のまま用いるという，半離
散化ウェーブレット変換により，同じ精度で復元が可能であることを示した．ここでは，
その結果に従い，Morlet ウェーブレットによる半離散化変換とその復元を手法として採用
する．つまり，元の時系列を X（t），パラメータ対（ s，n）におけるウェーブレット変換を
WnX (s) ，復元関数を X ∼（t）とすると，これらの関係は次式になる：
WnX (s)
N
√ X
= (1/ s)
X(n′ )Ψ∗ {(n − n′ )/s}
n′ =1
X ∼（t）= {2/(A + B)}
(1/C)
X
N X
X
n=1
j
WnX (2 j/4 )Ψ∗ {((t − n)/2 j/4 )}/2 j/8
WnX (2 j/4 )/2 j/8 .
j
ここで，C は定数， N は時系列の個数， ∗ は複素共役をとることを示し，離散化スケール
を定める j の範囲はナイキスト周波数を下限とし，時系列の長さにおいてウェーブレッ
トのパワーの 99.9 % が含まれるようなスケールを上限とするように選ぶ（Kumar and
Foufoula-Georgiou, 1994）．最後の式は復元において同心の半離散化ウェーブレット変換
だけを考慮すればよいことを示しており，この式による復元は離散フーリエ変換とその逆
変換を用いて効率的に算出できる．詳細は Kuroishi and Keller（2005）を参照されたい．
この Morlet ウェーブレットを用いたウェーブレット変換では，元の時系列において取り
出される最大の波長（フーリエ波長）が ω0 によって変化し，次式で与えられる（Torrence
and Compo，1998）：
4πs/(ω0 +
q
2 + ω20 ).
ここで， s は上記のスケール・パラメータである．したがって，上記の ω0 を用いる場合，
フーリエ波長（周期）は 1.157s になることに注意が必要である．
2.2
Halo ウェーブレットを用いた二次元半離散化ウェーブレットのフ
レーム
二次元データについて， Morlet ウェーブレットを二次元関数に拡張した Halo ウェーブ
レット（Dallard and Spedding, 1993）を母関数とする半離散化ウェーブレット展開による
フレームを適用する．このウェーブレットは，球対称な実関数であるため，二次元データ
において位相情報が失われるが，解析において等方性をもつことになる．Halo ウェーブ
レットは周波数領域において次式で表される：
Ψ̂(~
ω) := exp{−(| ω
~ | −ω0 )2 /2}/κ.
ここで，κ は正規化係数，ω
~ は二次元波数ベクトル，ω0 は基準波数の大きさであり，ここ
では前節の Morlet ウェーブレットの場合と同一の値にとることとする．
3
Halo ウェーブレットを母関数として，2.1 における場合と同じ半離散化ウェーブレッ
ト展開を行う．この場合，エネルギーが有限である任意の二次元関数に対する半離散化
ウェーブレット展開によるノルムは，対応する Morlet ウェーブレットによる半離散化
ウェーブレット・フレームにおけるフレーム定数について回転対称をとったものと一致す
ることを Kuroishi and Keller(2005) は示した．従って，Halo ウェーブレットを用いた半
離散化ウェーブレット展開はフレームをなすことになる．また，半離散化ウェーブレット
展開からの再構築は，各データ点において全てのウェーブレット係数を考慮する代わり
に，その点において同心の半離散化ウェーブレット変換だけを考慮することで，十分な精
度で復元されることも示された．従って，半離散化ウェーブレット展開からの再構築は，
近似として，有界な窓関数を持つ離散化フーリエ逆変換と等価になる．
3. 一次元半離散化ウェーブレットのフレームに基づくコ
ヒーレンシー分析
時系列データに Morlet ウェーブレットによる半離散化ウェーブレット展開を行うと，
（離散的な）スケール・パラメータと時刻パラメータの組みごとにウェーブレット変換が
得られ，これは複素数値をとり，振幅と位相を持つ．それは，パラメータに対応する時
刻と時間スケール（周期）において局在化された信号成分を示すものである．半離散化
ウェーブレット展開は冗長なフレームをなすため，得られたウェーブレット変換は，相互
に独立したものではなく，時間と周波数の両方の領域においてスケールに応じた広がりの
間で相関をもつことになる．
二つの時系列に対する相互ウェーブレット変換は，個々の時系列をウェーブレット変換
したものの複素共役積で定義される．フーリエ・スペクトルにおけるコヒーレンシーを
ウェーブレット変換にそのまま拡張すると，相互ウェーブレット変換のパワーを個々の
ウェーブレット変換のパワー積の平方根で割ったものになるが，その値は原理的に常に
1 となってしまい，コヒーレンシー分析に用いることができない．しかし，Torrence and
Webster（1999）は，ウェーブレット変換について時間・周波数領域のそれぞれにおいて
相関を持つ範囲での平滑化を適用するという考えを導入し，ウェーブレット二乗コヒーレ
ンシー（wavelet squared coherency；以下，WS C と記す）を提案した．二つの時系列 X
と Y について，スケール・パラメータ s ，時刻パラメータ n におけるウェーブレット変換
を，それぞれ，WnX (s)，WnY (s) とすると，WS C である R2n (s) は次のように与えられる：
R2n (s) =| S (WnXY (s)/s) |2 /{S (| WnX (s) |2 /s) · S (| WnY (s) |2 /s)}
ここで，WnXY (s) ≡ WnX (s) · WnY (s)∗ は相互ウェーブレット変換，S は平滑化関数であり，ス
ケール・パラメータに関する平滑化 S scale と時刻パラメータに関する平滑化 S time の連続
4
変換で定義される：
S (Wn (s)) = S scale (S time (Wn (s)))
S time (Wn (s)) | s = (Wn (s) · c1 exp(−t2 /2s2 )) | s
S scale (Wn (s)) | s = (Wn (s) · c2 Π(0.6s)) | s
ここで，Π は箱形関数，c1 と c2 は規格化定数である．この定義によって得られる WS C
は 0 から 1 までの値をとる．
二つの時系列の間における局在化信号の WS C の有意性は，信号やノイズの特性によっ
て異なる．本稿では，その閾値として 0.8 を採ることとする．これは，数多くの地球物理
現象にみられる赤色ノイズが 1 次の自己回帰モデルで十分モデル化できることを考慮し，
Grinsted et al.（2004）がモンテカルロ法によるシミュレーションで調査した有意水準の評
価結果に基づいて選んだものである．彼らが Morlet ウェーブレットについて上記の設定
による離散化ウェーブレット変換を対象に行った調査では，WS C の 95 % 信頼水準は，
いずれのスケール・パラメータ値の範囲においても，ほぼ 0.8 以上であることが明らかに
されている．
4. 二次元半離散化ウェーブレットの応用：日本周辺の海域
重力場モデルにおける局在化較差分離
日本列島は四つのプレートが収斂する境界に位置し，陸海における地形起伏が激しく，
重力場とジオイド形状が広い波長域に亘って大きく変化している．そのため，陸域のジオ
イドを重力手法によって詳細に決定する場合，周辺海域における重力場を正確かつ詳細に
把握する必要がある．
日本列島周辺海域においては，海上重力計による稠密な測定が利用可能であるが，それ
らの観測は主に 1970-80 年代に行われ，当時の航法精度の限界のため，船の速度に対する
エトベス効果（コリオリ力の鉛直成分）の補正などに伴う，観測毎に系統的な中・長波長
誤差が含まれているものがある．一方，高度計モデルも利用可能であるが，空間分解能が
十数 km 程度以上で中波長域の復元に限られ，また，沿岸部では精度が劣る．そこで，高
度計モデルを利用して，海上重力計による重力場モデルの局所的系統誤差を推定，除去す
ることで海域重力場の決定精度を向上させることが可能となる．
日本列島周辺の海域について，海上重力計によるモデルと高度計モデルの較差分布を入
力として，二次元半離散化ウェーブレット・フレームにより，異なるスケールでの局在化
較差信号に分離する．その局在化信号の分布特性に従い，海上重力計によるモデルに対す
る補正モデルの作成を行った． Fig. 1 に重力較差分布とそれから作成された補正モデルの
結果を示す．北海道のオホーツク海沿岸域と三陸北部沖にみられる顕著な補正分布は，特
定の海上重力観測航海の群によく対応している結果となった．
5
-60
-45
-30
-15
0
15
30
45
60
-30
50N
50N
45N
45N
40N
40N
35N
35N
30N
30N
25N
25N
20N
-25
-20
-15
-10
-5
0
5
10
15
20
25
30
20N
120E
125E
130E
135E
140E
145E
150E
120E
125E
130E
135E
140E
145E
Fig. 1. Marine gravity difference field between a model obtained using marine gravity
data and an altimetry-derived model shown in the left panel, and gravity correction model
in the right panel, both from Kuroishi and Keller (2005). Units for the scale above the
image are mGal.
Tidal Height (01:Aburatsubo)
3.0
2.9
Daily mean (m)
2.8
2.7
2.6
2.5
2.4
2.3
01/01/78
01/01/79
01/01/80
12/31/80
12/31/81
01/01/83
01/01/84
12/31/84
12/31/85
01/01/87
01/01/88
12/31/88
12/31/89
01/01/91
01/01/92
12/31/92
12/31/93
01/01/95
01/01/96
12/31/96
12/31/97
01/01/99
01/01/00
12/31/00
12/31/01
01/01/03
01/01/04
12/31/04
12/31/05
01/01/07
01/01/08
12/31/08
2.2
Date
Fig. 2. Time series of tidal heights at Aburatsubo for the period 1978-2008 after removing ocean tides and oceanic response to atmospheric pressure change estimated using
BAYT AP − G (taken from Kuroishi, 2012).
6
150E
5. 一次元半離散化ウェーブレットに基づくコヒーレンシー
分析の応用：東京湾近傍の潮位データの平滑化
国土地理院による油壺験潮場と海上保安庁による横須賀験潮所における潮位観測の１時
間間隔時系列について，二段階の平滑化を考える．まず，それぞれの潮位時系列に対し
て，横浜地方気象台における気圧測定の１時間間隔時系列を併せて用い， ABIC (Akaike’s
Bayesian Information Criterion; Akaike, 1980) に基づき，潮汐変化と lag つきの気圧応答
を推定する BAYT AP − G（Tamura et al., 1991）を適用して滑らかに変化する成分を取り
出す．つぎに，得られた成分について，時間に線形なドリフトを除いた上で一次元半離散
化ウェーブレット・フレーム表現を求め，油壺と横須賀の二点のウェーブレット係数間で
WS C を算出し，有意なコヒーレンシーとなる成分を除去することにより，平滑化された
時系列を得る．
coherency phase diff for 01 -73; start 19780101
180
44
2000
135
40
1000
500
32
200
28
100
24
50
20
16
20
12
10
90
Fourier period(day)
Scale integer
36
45
0
-45
-90
-135
5
8
4
-180
1096
2191
3287
4383
5479
6574
7670
8766
9862
10957
Elapsed day
smoothed squared coherency:smoothed tidal height
1.0
44
2000
40
1000
500
32
200
28
100
24
50
20
16
20
12
10
8
5
0.8
Fourier period(day)
Scale integer
36
0.5
0.2
0.0
4
1096
2191
3287
4383
5479
6574
7670
8766
9862
10957
Elapsed day
Fig. 3. Wavelet squared coherency (bottom panel) and its phase difference (top panel)
between two time series of tidal heights at Aburatsubo and at Yokosuka for the period
1978-2008 (taken from Kuroishi, 2012). BAYT AP − G was applied in advance to the time
series for removing ocean tides and oceanic response to atmospheric pressure change. The
Morlet wavelet is used in the analysis. Horizontal axes are the elapsed day since January 1,
1978 with ticks at a one-year interval and corresponding Fourier periods in days are given
on the right vertical axes. Phase differences are given in degrees.
7
Aburatsubo tidal height corrected 1978-2008
2.65
corrected tidal height
11-day averaged corrected height
Tidal height (m)
2.6
2.55
2.5
2.45
2.4
79
82
85
88
91
94
Year
97
00
03
06
09
Fig. 4. Time series of mean tidal heights at Aburatsubo for the period 1978-2008
smoothed by removing the components coherent to those at Yokosuka, based on WS C
analysis (taken from Kuroishi, 2012). Crosses are daily values and the solid curve is
smoothed by taking 11-days moving average of the daily values.
油壺における 1978 年から 2008 年の潮位について，潮汐と気圧応答を除去した時系列
を Fig. 2 に示す．１年を周期とする三角波状の変動が，一定でない振幅で顕著にみられ，
海水のステリック変動であると考えられる．横須賀についても同様の処理を行い，二つの
時系列の間で，第 3 章に示す WS C を求めた結果をその位相差とともに Fig. 3 に示す．年
周スケールにおいては，位相差がなく，安定して高いコヒーレンシーを示しているほか，
広汎なスケールにおいて有意なコヒーレンシーとなっている．有意なコヒーレンシー成分
を除去することで平滑化された，油壺の潮位時系列を Fig. 4 に示す．一部の期間において
は細かな変動が含まれるが，全体として緩やかな上昇を示し，その緩やかな変化は，横須
賀の平滑化された潮位時系列と共通していることが分かった．
1997 年以降においては，別種の測地学的手法によって二つの験潮場とも地盤の上下変
動が連続的に観測されている．その結果と併せると，Fig. 4 にみられる緩やかな変化は，
1997 年以降においては， 2003 年頃を境として異なる速度での地盤の隆起と，ほぼ一定速
度の海面位上昇からなると推定される．
謝辞 大阪教育大学の芦野隆一先生には，本シンポジウムへ参加する機会をご提供いた
だきました．また，一部の図の作成には，GMT（Wessel and Smith, 1991）を用いました．
ここに記して謝意を表します．
参考文献
[1] Kuroishi, Y. and W. Keller, Wavelet approach to improvement of gravity field-geoid
modeling for Japan, J. Geophys. Res., 110, B03402, DOI:10.1029/2004JB003371,
2005.
8
[2] 黒石 裕樹, 平均海面位変化と上下変動の検出に向けた験潮データの平滑化手法，測地
学会誌，58, (2012), 27–42.
[3] Daubechies, I., Ten lectures on wavelets, Soc. for Ind. and Appl. Math., Philadelphia,
(1992), 53–78.
[4] Kumar, P. and E. Foufoula-Georgiou, Wavelet analysis in geophysics: An introduction,
in Wavelets in Geophysics, edited by E. Foufoula-Georgiou and P. Kumar, Academic
Press, (1994), 1–43.
[5] Torrence, C. and G.P. Compo, A practical guide to wavelet analysis, Bull. American
Meteorological Society, 79, (1998), 61–78.
[6] Dallard, T. and G.R. Spedding, 2-D wavelet transforms: generalization of the Hardy
space and application to experimental studies, Eur. J. Mech., B/Fluids, 12-1, (1993),
107–134.
[7] Torrence, C. and P.J. Webster, Interdecadal changes in the ENSO-monsoon system, J.
Climate, 12, (1999), 2679–2690.
[8] Grinsted, A., J.C. Moore, and S. Jevrejeva, Application of the cross wavelet transform
and wavelet coherence to geophysical time series, Nonlinear Processes in Geophysics,
11, (2004), 561–566.
[9] Akaike, H., Likelihood and Bayes Procedure, in Bayesian Statistics, edited by J.M.
Bernardo, M.H. DeGroot, D.V. Lindley and A.F.M. Smith, University Press, Valencia,
(1980), 143–166.
[10] Tamura, Y., T. Sato, M. Ooe, and M. Ishiguro, A procedure for tidal analysis with a
Bayesian information criterion, Geophys. J. Int., 104, (1991), 507–516.
[11] Wessel, P. and W.H.F. Smith, Free software helps map and display data, EOS Trans.
AGU, 72, (1991), 441,445–446.
黒石 裕樹 (国土交通省国土地理院)
〒305-0811 茨城県つくば市北郷 1
E-mail: [email protected]
9
脳情報解析およびバイオインフォマティクスに
おけるスパース信号処理
∗
池田 和司 ∗ 作村 諭一 †
奈良先端科学技術大学院大学情報科学研究科 † 愛知県立大学情報科学部
概要.
圧縮センシングなどスパース信号処理の応用は様々であるが，本発表では脳情報
解析およびバイオインフォマティクスに応用した例について紹介する．脳情報解析では脳
波を用い，sparse logistic regression でひらめき判別器を構成することで，ひらめきに関
連する脳活動の時空間パターンを抽出した．バイオインフォマティクスでは酵素阻害剤を
用いた実験結果を用い，圧縮センシングにより軸索伸長に関連する酵素を同定した．
Sparse Signal Processing in Brain Informatics
and Bioinformatics
Kazushi Ikeda∗ Yuichi Sakumura†
∗
Nara Institute of Science and Technology † Aichi Prefectural University
Abstract.
Compressed sensing and other sparse signal processing techniques have been
applied to various fields. In this study, they were applied to brain informatics and bioinformatics. In the former case, we made a sparse logistic regressor that discriminates electroencephalogram (EEG) signals during insight from those during non-insight in order to
extract the spatio-temporal pattern for insight. In the latter case, we identified the kinases that regulate neurite growth by solving an inverse problem that models compounds,
kinases and neurite growths using compressed sensing.
1. はじめに
近年注目を浴びている圧縮センシングなどのスパース信号処理は，信号のスパース性す
なわち非ゼロ要素が少数であることを仮定して不良設定問題を解く手法である [1–4]. そ
の応用は様々であるが，本発表では脳情報解析およびバイオインフォマティクスに応用し
た例について紹介する．
脳情報解析への応用においては，ひらめきに関与する脳活動の抽出を試みた [5]. ひらめ
きとアルゴリズム的の二通りの解き方があるタスクについて頭皮脳波 (electroencephalo-
gram, EEG) を測定し，EEG 信号を入力とするひらめき判別器を構成した．この判別器に
sparse logistic regression (SLR) [4] を用いることで，SLR の係数の絶対値が大きなところ
には情報が含まれていると考えることができる．これにより，ひらめきに関与する時空間
パターンを抽出した．
バイオインフォマティクスへの応用においては，神経軸索の伸長に関与する酵素の同定
1
を試みた [6]. これは酵素阻害剤を用いた時の神経伸長データから酵素を特定する問題で
あり，一般に酵素阻害剤は複数の酵素の作用を阻害するため，線形の逆問題として定式化
される．これを圧縮センシングを用いて解き，神経伸長に関連する酵素を同定した．
2. 脳情報解析への応用
ひらめき発生時には前頭状皮質や前頭葉前部の活動に変化が生じると言われている
が [7], 従来の研究は通常状態における脳活動との比較であるため，ひらめきと論理的な思
考過程との切り分けができていなかった．そこで本研究ではひらめきとアルゴリズム的の
二通りの解き方があるアナグラム・テストを用い，両解法における脳活動の差異に着目す
ることで，ひらめきに関与する脳活動の時空間パターンを抽出した [5].
2.1 方法
11 名の被験者に対して 4 文字のアナグラム・テストを実施した (Fig. 1). 被験者には問
題が解けた時にボタンを押させ，その後にどちらの方法で解いたかを口頭で報告させた．
実験中，被験者の EEG を測定した (15 チャンネル，サンプル周波数 200Hz).
Fig. 1. The flow of the anagram test.
上記実験で得られた EEG 信号をひらめきとアルゴリズム的の二つに分け，これらを分
類するひらめき判別器を SLR によって構成した．まず，ボタン押しの時刻を時刻 0 とし
て時間軸を揃え，750ms 前–500ms 前，500ms 前–250ms 前，250ms 前–0ms の 3 つの時間
窓を用意した．各時間窓において，15 チャネルの各々について EEG 信号のバンドパワー
(α 波 8–13Hz, β1 波 14–20Hz, β2 波 21–30Hz, γ 波 36–40Hz) を求め，180 次元のベクト
ルを入力とした．
2
2.2 結果
構成された SLR の係数は，各脳部位/時刻/バンドが持つひらめきに関与する量を表して
いると考えられる．これを図示したものが Fig. 2 であり，各チャネルのバーは左から α
波, β1 波, β2 波, γ 波を表している．
Fig. 2. Spatio-temporal pattern extracted by SLR.
3. バイオインフォマティクスへの応用
神経軸索の伸長には多くの酵素が関与しており，酵素阻害剤の投与により伸長が亢進ま
たは抑制される．しかし一般に酵素阻害剤は複数の酵素の作用を阻害するため，このデー
タからはどの酵素がどれだけ伸長に作用するかは自明ではない．医学的生物学的見地から
は，伸長に関与する酵素を同定し，効率よく軸索を伸長する酵素阻害剤の組み合わせを開
発することが望まれる．
一般に酵素阻害剤数は酵素数よりも少ないため，これは不良設定問題となる．そこで圧
縮センシングを用い，伸長に関与する酵素を同定した [6].
3
3.1 方法
この現象を以下のようにモデル化する．酵素阻害剤 m により酵素 n が Amn だけ阻害さ
れるとし，酵素 n は神経軸索を xn だけ伸長させるとする．この作用が線形であると仮定
すれば，阻害剤 m を与えた時の伸長 ym は
ym =
(3.1)
∑
Amn xn
n
と表され，行列・ベクトルで表記すれば y = Ax という線形方程式となる．一般に酵素阻
害剤数 M は酵素数 N よりも小さいため，これは不良設定問題となる．そこで軸索伸長に
関与する酵素は少数であると仮定し，圧縮センシングの手法でこれを解いた．すなわち
(3.2)
min kxk
s.t. y = Ax
という最適化問題を解くことで x を求め，酵素 n の軸索伸長への影響度を求めた．
なお，本研究で用いた実験データは酵素阻害剤 35 種類，酵素 139 種類であり，共同研
究者の Vance Lemmon 教授のグループから提供を受けたものである．
3.2 結果
Leave-n-out 法を用い (青 n = 1, 赤 n = 2), 各酵素の重要度 (xn ) の分布を調べたところ，
Fig. 3 のようになった (一部抜粋). 分布は n にほとんど依存していないことがわかる．
4. まとめ
本発表ではスパース信号処理を脳情報解析およびバイオインフォマティクスに応用した
例を紹介した．ここでの例ではいずれも正解が存在していないため，この手法がどれだけ
有効であったかを確認できていない．今後，数学的な正当性の確認や人工データによる確
認，生物学的な妥当性の検証などが必要である．
参考文献
[1] D.L. Donoho, M. Elad (2003): Optimally sparse representation in general (nonorthogonal) dictionaries via l1 minimization, Proc. Nat. Aca. Sci., 100, 2197–2202.
[2] D.L. Donoho (2006): Compressed sensing, IEEE Trans. Information Theory, 52, 1289–
1306.
4
Fig. 3. Importance of xn in leave-one-out method and leave-two-method.
[3] E. Candes, T. Tao (2007): The Dantzig selector: Statistical estimation when p is much
larger than n, Ann. Statist., 35, 2313–2351.
[4] O. Yamashita, M. Sato, T. Yoshioka, F. Tong, Y. Kamitani (2008): Sparse estimation
automatically selects voxels relevant for the decoding of fMRI activity patterns, Neuroimage, 42, 1414–1429.
[5] R. Kondo, T. Shibata, N. Oosugi, K. Ikeda (2011): Spatio-temporal transition of EEG
activities during insight occurrence, Proc. JNNS 2011, P3-12.
[6] 丸野, 他 (2011): 圧縮センシングによる神経軸索の伸長を制御する因子の同定に関す
る検討, IBIS2011, D-129.
[7] B. Sheth, et al. (2009): Posterior beta and anterior gamma oscillations predict cognitive
insight, J. Cognitive Neuroscience, 21, 1269–1279.
5
池田 和司 (奈良先端科学技術大学院大学情報科学研究科)
〒630-0192 奈良県生駒市高山町 8916-5
E-mail: [email protected]
作村 諭一 (愛知県立大学情報科学部)
〒480-1998 愛知県長久手市茨ヶ廻間 1522-3
E-mail: [email protected]
6
佐藤超関数による
Shannon ー染谷の標本化定理の拡張と
Ramanujan の積分公式
∗
吉野 邦生 ∗
東京都市大学知識工学部自然科学科
概要.
佐藤超関数の理論を応用し、Shannon ー染谷の標本化定理を拡張する。更に,
Shannon ー染谷の標本化定理と Ramanujan の積分公式の関係を明らかにする
A Generalization of Shannon - Someya’s Sampling Theorem by
Sato’s Hyperfunction Theory and
Ramanujan’s Integral Formula
Kunio Yoshino∗
∗
Tokyo City University
Abstract.
We generalize Shannon-Someya’s sampling theorem by using the theory
of Sato’s Hyperfunctions. Moreover we will clarify the relationship between ShannonSomeya’s sampling theorem and Ramanujan’s integral formula.
1.
Shannon ー染谷の標本化定理
佐藤超関数に対する Paley-Wiener の定理を応用することにより次の定理を得る.
¶
³
定理 整関数 f (z) が次の条件を満たしているとする:
∀ε > 0, ∃Cε > 0 s.t.
| f (z)| ≤ Cε exp(b|y| + ε|z|),
(∀z = x + iy ∈ C).
もし 0 ≤ b < π であると, 次が成立する。
∞
X
f (n) sin π(z − n) −δ|n|
e .
δ→0
π(z
−
n)
n=−∞
f (z) = lim
µ
´
1
2. ラマヌジャンの積分公式
¶
³
ラマヌジャンの積分公式
Z
∞
a−1
u
0
∞





π f (−a)
X
n
f
(n)(−u)
du =
.






sin(πa)
n=0
µ
´
この講演では, 佐藤超関数の概念の簡単な紹介をしてから, 標本化定理とラマヌジャンの積
分公式との関係を明らかにする. その際, 信号データの Z-変換によるラマヌジャンの積分
公式の正確な定式化についても紹介する.
参考文献
[1] G. H. Hardy, Ramanujan, Chelsea, New York (1978)
[2] K. Yoshino, Liouville type theorem for entire functions of exponential type,
Complex Variables, 5(1985), p. 21-51.
[3] 吉野邦生ー荒井隆行, ディジタル信号と超関数, 海文堂 (1995)
[4] 吉野邦生, ディジタル信号と母関数, 数学セミナー 10 月号,
日本評論社 (2001), p. 24-27
[5] K. Yoshino, The relation among Shannon’s sampling theorem, Ramanujan’s integral
formula and Plana’s summation formula, Complex Analysis and Potential Theory,
Proceedings of the Conference Satellite to ICM 2006 (Gebze Institute of Technology,
Turkey) (Tahir Aliyev Azeroǧlu and Promarz M. Tamrazov eds.), World Scientific Publishing Co. Pte. Ltd (2007), p. 191-197.
[6] 吉野邦生, ヘビサイドの超関数を知る, 数理科学 5 月号, (2009), p. 54-55
[7] 吉野邦生, シャノンー染谷の標本化定理とラマヌジャンの積分公式,
東京都市大学共通教育センター紀要, 3 月（2010), p. 1-21
吉野 邦生 (東京都市大学 知識工学部 自然科学科)
〒158-8557 東京都世田谷区玉堤 1-28-1
E-mail: [email protected]
2
サンプリング値の補間による未知関数の
再構成
∗
岡田 正已 ∗
首都大学東京 理工学研究科 数理
概要.
複数次元の不規則な格子点集合で与えられたサンプリング観測値から、元の関数
をできるだけ忠実に再現するために、正定型関数を用いる方法を概説し、今後の問題と応
用への課題を議論したい。
Reconstruction by interpolation of sampled scattered data
Masami Okada∗
∗
Tokyo Metropolitan University
Abstract.
We review a recent approach of using positive definite functions to reconstruct
given multivariate functions from sampled data on scattered infinite points. We would like
to discuss future relevant topics including related applications
1. 動機
整数点全体 Z 上で観測された値 ( fk ) に対して
F(x) =
∑
fk
sin π(x − k)
π(x − k)
とおけば, たとえば、 ( fk ) ∈ l2 のとき、F は直線全体 R で二乗可積分 になり、F(k) = fk で
ある（補間）。
逆に、少し一般にして、関数 f (x) が先に与えられたとき、そのフーリエ変換が遠方で
十分に早く減少すれば、関数 S f,Z (x) =
∑
f (k2−n ) sinc(2n x − k) は、 f の滑らかさに比例し
た良い近似関数となる（近似）。
以上の事柄を規則性のない無限の点集合 X に一般化し、多次元ユークリッド空間で考
えたら、どの程度のことが云えるのか調べたい。
比較的最近、よく研究された RBF (Radial Basis Function) の方法を参考にする。
2. positive definite functions
適当に滑らかな偶関数で positive definite (正定型）なもの ϕ を考える。定義は、（任意
の個数の）任意の有限集合 {ξk } と同じ個数の 任意の実数 αk にたいして、二重和が非負
1
0≤
∑
αi α j ϕ(ξi − ξ j ) ただし、ゼロになるのは、全ての αi = 0 のときに限る、である。こ
のとき、与えられた点集合 X = {x1 , · · · , xN } が有限集合のときには、positive definite ゆえ
∑
に行列 (ϕ(ξk − ξ j )) が正則なので、各 x j で観測値 b j をもつ関数 f を f (x) = ak ϕ(x − xk )
として、未知数ベクトル a = (ak ) が求まる。ただし、有限の個数 N といっても、たとえ
ば N = 104 なら、実際の数値計算は工夫を要す。
では、点集合 X が可算無限集合のときは、どうするか、これを課題とする。当然、X は
（無限遠以外の）集積点は持たないものとする。
positive definite function の例としては、一変数の、sinc(x) とか、適当な一次結合スプラ
インとかの直積 タイプ、ガウス関数、(1+ k x k2 )−1/2 とか、コンパクトサポートを持つよ
うな radial 関数などがある。それらは実質的に殆どのところで、実解析的である。以下、
ϕ(0) = 1 とし、コンパクトサポートを持っているとする。
3. 関数空間
sinc 関数によるサンプリング定理で扱った、帯域制限のついた関数のクラスをかなり一
般にして、次の関数の集合と内積を導入する。
∑
H0 = { f (x) = αk ϕ(x − ξk ) | ξk ∈ Rd , αk ∈ R} . まずは、有限和を考え、関数 g(x) =
∑
∑
βl ϕ(x − ηl ) との内積を ( f, g) = αk βl ϕ(ξk − ηl ) で定義する。これでノルムが定まり、こ
れに関して H0 の閉包を H とすれば、これはヒルベルト空間になる。 f (x) = ( f (·), ϕ(· − x))
であり、 ϕ は H の再生核である。
ここで、
∫
( f, g) =
b
f (ω)b
g(ω)(b
ϕ(ω))−1 dω
d の遠方での減少度に応じて、H は良いソボレフ空間に属することが知られ
であり、ϕ(ω)
ている。
4. Cardinal (Lagrange) function
補間するために、X = {x1 , x2 , · · · } の各点 xk に対応する基底関数 u∗k があればよい。ただ
し、u∗k は H の部分空間 VX = { f ∈ H | f (x) =
∑
β j ϕ(x − x j ), x j ∈ X} に属し、u∗k (x j ) = δk j で
ある。存在証明は、例えば、シュミットの直交化による。近似計算するためには、Greedy
アルゴリズムも適用可能。この基底関数を用いれば、 X の各点 xk で 観測値 f = ( fk ) を持
つ関数で VX に属するノルム最小の関数 S f,X は
∑
S f,X (x) =
fk u∗k (x)
で求まる。
2
5. 近似計算
実際に u∗k を近似計算で求めるには、適当な Greedy アルゴリズムも考えられる。近似
誤差評価であるが、X のメッシュの細かさを hX = sup x inf j k x − x j k とおくとき、次が成
立すると期待される。ϕ が十分滑らかで、関数 f がソボレフ空間 W ps (Rd ) に属するとき (但
し、 dp < s) 、小さな値 hX にたいして
k f − S f,X kL p ≤ C hX s k f kW ps .
6. 課題
1. 一般の領域で考察すること。
2. u∗k のより詳しい性質を調べること。
3. 実用 的 な 計 算 ア ル ゴ リズ ム を与 える こ と。例 え ば 参 考文 献 には sn = sn−1 +
Pn fn−1 , fn = fn−1 − Pn fn−1 で lim sn として u∗k を近似計算するアルゴリズムが
述べられている。但し、勿論、X が無限集合の時には、収束することの証明はみあ
たらない。
4. X が適当な対称性を持った場合や、準均一のとき、逆に非常に不均一のときの考察。
5. ϕ が positive definite でなく、conditionally p.d. のときとの比較。
謝辞 芦野先生には、研究会の企画など、お世話になり感謝します。内容は、カーキャ
チャリアン教授とジャファール教授との共同研究に基づきます。The author thanks Prof.
R. Ashino for organizing the conference. The contents are based on an ongoing joint work
with Profs. G. Kerkyacharian and St. Jaffard.
参考文献
[1] H. Wendland, Scattered Data Approximation, Cambridge Monographs on Applied and
Computational Mathematics, Cambridge University Press, 2005 (Paperback 2010)
岡田 正已 (首都大学東京 理工学研究科 数理)
〒192-0397 東京都八王子市南大沢 1-1
E-mail: [email protected]
3