Document

池内研究室
Computer Vision Laboratory
Copyright: Daisuke Miyazaki
http://www.cvl.iis.u-tokyo.ac.jp/
池内研究室
Computer Vision Laboratory
Numerical Recipes in C
[日本語版]
初版
宮崎大輔
2004年5月13日（木）
PBVセミナー
池内研究室
Computer Vision Laboratory
第12章フーリエ変換
Fourier transform spectral methods
8th Physics Based Vision Seminar
2004/May/13
Fourier(フーリエ)変換法(Fourier transform
methods), スペクトル法(spectral methods)
• 関数h(t),時間t,時間領域(time domain)
振幅H(f),周波数f(-∞<f<∞),周波数領域(frequency
domain)
h(t)とH(f)は同じ関数の異なる表現
4
8th Physics Based Vision Seminar
2004/May/13
畳込みと相関
• 二つの関数の畳込みまたは合成積
(convolution)g*h
• 二つの関数の相関(correlation)
tは遅れ(lag), H(-f)=H*(f)
• 自己相関(autocorrelation):関数とそれ自身との相
関
5
8th Physics Based Vision Seminar
2004/May/13
パワースペクトル密度
• 全パワー(total power):
Parseval(パーセバル)の定理
(Parseval’s theorem)
• 片側パワースペクトル密度
(one-sided power spectral
density)
• 両側パワースペクトル密度
(two-side power spectral
density)
6
8th Physics Based Vision Seminar
2004/May/13
Nyquist(ナイキスト)の臨界周波数(Nyquist critical frequency),
標本化定理(sampling theorem), エイリアシング(aliasing)
• 標本化周期Δ
• Nyquistの臨界周波数
• 標本化定理
• エイリアシング
– 周波数帯域(-fc,fc)の外側のど
んな周波数成分も，この帯域
内に折り返される
– もし間隔Δで標本化した連続
関数h(t)の周波数帯域がfcよ
り小さい周波数だけなら，関
数h(t)は標本hnから完全に決
定される
7
8th Physics Based Vision Seminar
2004/May/13
離散Fourier変換(discrete Fourier transform)と離散
逆Fourier変換(discrete inverse Fourier transform)
• 離散Fourier変換
離散的な値
だけで推定(-fc～fc)
• 離散逆Fourier変換
• 標本値はN個，標本化間隔はΔ
8
8th Physics Based Vision Seminar
2004/May/13
高速Fourier変換(fast Fourier transform, FFT)
• 離散Fourier変換はO(N2)？
→No! FFTならO(Nlog2N)
• Danielson and Lanczosの補題
Fk  F
even
k
W F
k
odd
k
W  e 2i N
Fkevenは偶数番目の成分
Fkoddは奇数番目の成分
長さNの離散Fourier変換は2個
の長さN/2の離散Fourier変換の
和で書き表せる
• これを帰納的に使える→Nが２の
累乗の場合，どんどん，半分に分
割していける→長さが1のFourier
変換は単なる恒等変換で，入力
をそのまま出力にコピーするだけ
• ビット反転(bit reversal)
a
a
b
e
c
c
d
g
e
b
f
f
g
h
d
h
1 2 3 4 5 6 7 8
a b c d e f g h
奇数
aceg
偶数
bdfh
奇数
ae
偶数
cg
奇数
bf
偶数
dh
a
c
b
d
e
g
f
h
9
8th Physics Based Vision Seminar
2004/May/13
実関数のFFT
• FFTは複素数演算をしているが，虚部を0に設定す
ることにより，N個の実数のデータに対しても，あたり
まえに計算できる
• 出力は複素数の可能性もあるが，出力の自由度はN
のまま変わらない(複素数だったら，実部と虚部合わ
せて自由度2Nになりそうなものである)
• よって，複素数用のFFTルーチンを実数データに使う
のは，実行時間と格納場所の効率が悪い
• 効率的に計算する方法があるが詳細は省略
10
8th Physics Based Vision Seminar
2004/May/13
サイン変換(sine transform),
コサイン変換(cosine transform)
N 1
f j sin(jk N )
• サイン変換： Fk  
• 計算方法はFFTを拡張して
j 1
N 1
求めることもできるが，2倍
• コサイン変換： Fk  
f j cos(jk N )
効率が悪い→もっといい方
j 0
法があるが詳細は省略
• サイン変換：サインだけを使
う
• コサイン変換：コサインだけ
を使う
• 通常のFFT：サイン・コサイ
ン両方を使う
11
8th Physics Based Vision Seminar
2004/May/13
FFTによる畳込みと逆畳込み
•
•
2関数r(t),s(t)の畳込みr*sは数学的
にはs*rと等しい
しかし，多くの応用場面では，以下の
ような意味と性格を持たせる
– sは信号(signal)またはデータの流れ
(stream)，時間について限りなく続く
– rは応答関数(response function)，山
形の関数で両側で0
– 畳込みの役割：s(t)をr(t)でにじませる
•
畳込み
•
rを持続時間(duration)Mの有限イン
パルス応答(finite impulse response,
FIR)フィルタと言う
離散畳込み定理(discrete
convolution theorem)
•
周波数領域で，sjの離散Fourier変換
Snと，rkの離散Fourier変換Rnの積
12
8th Physics Based Vision Seminar
2004/May/13
FFTによる畳込みと逆畳込み
• 畳込みをすると，両端が汚染される
• 両端に0を詰める
• m+とm-のうち大きい値の分だけ0を詰めればよい
13
8th Physics Based Vision Seminar
2004/May/13
FFTによる畳込みと逆畳込み
• 逆畳込み(deconvolution)：測定装置の分解能が完
璧でないためにデータににじみが生じたとき，にじみ
を表す応答関数が既知なら，にじみを取り除くことが
できる
• 周波数領域で単なる割り算
• 入力データの雑音や応答関数rkの精度に非常に敏
感であるという欠点を持つ
14
8th Physics Based Vision Seminar
2004/May/13
FFTによる相関と自己相関
• 2個の標本化された関数gk,hkの離散相関は，これら
が周期Nを持つとき
遅れ(lag)j
• 関数gと，時間の遅れjがあったときの関数hが，よく
似ていると，そのjの値で相関は大きくなる
• 離散相関定理(discrete correlation theorem)
Gk,Hkはgk,hkの離散Fourier変換
*は複素共役を意味する
• 離散自己相関(discrete autocorrelation)：自分自身
との相関
15
8th Physics Based Vision Seminar
2004/May/13
最適フィルタ(optimal filter), Wienerフィルタ
•
汚れていない(unccorrupted)信号u(t)を求めたいが，測定装置から出てくるのは汚
れた(corrupted)信号c(t)．c(t)は，真の信号u(t)を応答関数r(t)でにじまされた信号
s(t)と，雑音n(t)の和
•
真の信号Uの推定値は以下のようになる
•
•
φ(f)は最適フィルタ, c(t)⇔C(f), s(t)⇔S(f), n(t)⇔N(f), r(t)⇔R(f), u(t)⇔U(f)
応答関数r(t)は既知とする
また，にじんだ信号Sと雑音Nを分離できない
といけない
– 単純な方法として，まず，データc(t)を長い時間
にわたって測定し，そのパワースペクトル密度
をプロットする
– 低周波成分には，真のデータが多く載っていて，
高周波成分には雑音が多く含まれるだろう
– グラフを見て，適当にSとNのパワーを計算して
やればいい
16
8th Physics Based Vision Seminar
2004/May/13
FFTによるパワースペクトルの推定
• パワースペクトル(パワースペクトル密度, power spectral density, PSD),
ペリオドグラム(ピリオドグラム, periodogram)推定量P
wjは図にあるような窓関数, 周波数fkの強さがP(fk)となる
• 入力データcjに窓関数wjを掛ける⇔周波数領域での畳込み
17
8th Physics Based Vision Seminar
2004/May/13
ディジタルフィルタ
• 広域フィルタ(high-pass filter):低い周波数の雑音を除去
• 低域フィルタ(low-pass filter):高い周波数の雑音を除去
• 帯域フィルタ(band-pass filter):ある周波数帯だけの信号を取
り出す
• ノッチフィルタ(notch filter):ある周波数帯だけ除去
– Fourier領域でフィルタを掛ける：簡単．データをFFTして，フィルタ関数
H(f)をかけ算して，逆FFTで時間領域に戻す
– 時間領域でフィルタを掛ける：次スライドで説明
18
8th Physics Based Vision Seminar
2004/May/13
時間領域でのディジタルフィルタ
• 一般の線形フィルタ
• 例：特定の周波数w0のまわりの
狭い周波数帯だけを除去するノッ
チフィルタ
入力：データ点の列xk
出力：データ点の列yn
フィルタの応答を決める係数
ck,dj(固定)
新しい出力値を，現在およびM個
の過去の入力値，および自分自
身の過去のN個の出力値から作
る
• 係数ck,djとフィルタ応答関数H(f)
との関係
• H(f)からck,djを求める方法の詳細
は省略
19
8th Physics Based Vision Seminar
2004/May/13
線形予測
•
– 時系列の次の値ynを今までのN個の値yn-j(j=1,…,N)から予測
(predict)する式
– xnは時間ステップnでの予測の残差(外れ, discrepancy)，予測値（和
の部分）と真値ynとの差
– |xn|<<|yn|が成り立つのが良い
– 線形予測(LP)係数(linear prediction (LP) coefficients)d1,…,dN
• 時系列の未来を過去の記録から予測できる．未知の未来の
外れxiを仮に0として，この式をそのまま未来に向かって適用
し続ける．特に滑らかに振動する信号に対して正確に補外で
きる．
• LP係数を求める方法は省略
20
8th Physics Based Vision Seminar
2004/May/13
線形予測符号化
• LP係数を使ってデータを圧縮できる=線形予測符号
化(linear predictive coding, LPC)
–
–
–
–
LP係数の個数N
N個のLP係数
最初のN個のデータ点
それ以降のデータ点は外れ（残差）xiだけ記録する
• この外れは小さい値が多いので，Huffman(ハフマン)符
号化(Huffman coding)などを使って圧縮できる
21
8th Physics Based Vision Seminar
2004/May/13
２次元以上のFFT
• 複素関数h(k1,k2)が２次元の格子で与えられている
とき，その２次元離散Fourier変換H(n1,n2)は，同じ
格子上で与えられる複素関数で，次式のような感じ
になる
• つまり，２次元FFTは，元の関数の添字ごとに１次元
ずつFFTしたもの
• ３次元以上でも同じような感じ
22
8th Physics Based Vision Seminar
2004/May/13
フーリエ変換とラプラス変換とウェーブレット変換
• フーリエ変換
1 
1
iyx
f ( x) 
g ( y )e dy
g ( y) 

2 
2
• ラプラス変換
1 
 yx
f ( x) 
g
(
y
)
e
dy

2 
• ウェーブレット変換
1
 t b 
(W f )(b, a) 
f (t ) 
dt


a R
 a 
なお  はウェーブレット  の複素共役



f ( x)e iyxdx
23
池内研究室
Computer Vision Laboratory
第14章データのモデル化
Modeling of data
8th Physics Based Vision Seminar
2004/May/13
最小二乗法とカイ２乗当てはめ
• N個のデータ点(xi,yi)i=1,…,NをM
個のパラメータaj,j=1,…,Mを持つ
モデルで当てはめる
• あるいはσiを「重み」と見なしても
いいだろう
• χ2を各パラメータakで微分すると
• 当てはめ値を求めるには
とすれば良い
• データ点(xi,yi)がそれぞれ標準偏
差σiを持つとする（測定誤差が正
規分布であると仮定）
• 一般の場合の以下のカイ２乗を
最小にすることによってモデルパ
ラメータの最尤推定値(maximum
likelihood estimator)を得る
• この，M個の未知数akについての
M個の非線形方程式を解けば良
い
• σが分からない場合でσを計算し
たいときは，まず，σiに適当な定
数を代入（σi=σ）する．例えば
σi=1とする．これでパラメータを求
めたあと，σを以下の式で計算す
れば良い
25
8th Physics Based Vision Seminar
2004/May/13
データの直線への当てはめ
=線形回帰(linear regression)
• N個のデータ点(xi,yi)を直線
に当てはめる
• 以下を最小化する
• 最小値ではχ2のa,bについ
ての微分はゼロ
• 簡単にすると
ただし
• これを解くと
26
8th Physics Based Vision Seminar
2004/May/13
つづき
• a,bそれぞれの分散
• aとbの共分散(covariance)
 N 2 2 
Q  Q
, 
2 
 2
1  t a 1
Q ( a, x ) 
e t dt

( a ) x

( z )   t z 1e t dt
0
• aの不確かさとbの不確かさとの
間の相関係数(-1～1)
• rabが正ならaとbの誤差はほとん
ど同符号，負ならほとんど異符号
• Qが0.1より大きいなら当てはまり
の良さはかなり信頼できる
• 0.001より大きいなら誤差は正規
分布ではない．あるいは，多少の
誤差を気にしない人にとっては十
分妥当な当てはまりである
• 0.001より小さいなら，モデルか評
価方法かその両方が疑わしい
27
8th Physics Based Vision Seminar
2004/May/13
一般の線形最小２乗法
•
M-1次の多項式
•
に当てはめたい．関数はxのべき乗で
なくとも，サインやコサインでも構わな
い
一般形は，
•
•
X1(x),…,XM(x)は基底関数(basis
functions)
導出は省略するがAa=bをaについて
解けば解が得られる
デザイン行列(design matrix)A（右上
の図）はN×Mの行列
•
•
•
•
•
σiはi番目のデータの測定誤差
bは長さNのベクトル
aはM個のパラメータa1,…,aM
•
•
aについて特異値分解などで解く（σi
は既知，もしくは1とおく）
ATAの逆行列をCとおくと
これは，Cの対角要素は，当てはめ
たパラメータaの分散を意味する
Cjkの非対角要素はajとakの共分散
28
8th Physics Based Vision Seminar
2004/May/13
非線形最小２乗法を解くためのLevenbergMarquardt法（レヴェンバーグ-マーカート法）
•
•
当てはめたいモデル
χ2評価関数
•
•
•
•
•
•
M次のベクトルbのk番目βk
M×M行列Aのkl成分αkl
曲率行列(curvature matrix)Aは
Hesse行列を2で割った行列
ベクトルdを現在の解に加えて次
の解を得る
anew=aold+d
aはM個のパラメータ
•
•
•
•
Ad=bをdについて解いて反復計
算→逆Hesse法
d=定数×bで反復計算→最急降
下法
Ad+λId=bをdについて解いて反
復計算→Levenberg-Marquardt法
Iは単位行列
λが大きいと最急降下法に似て
て，小さいと逆Hesse法に似てる
λは反復のたびにうまく変化させ
る（徐々に減らしていく）
29
8th Physics Based Vision Seminar
2004/May/13
モンテカルロシミュレーションによる信頼限界
(confidence limits)
• パラメータを推定する
• 測定装置の測定誤差に
できるだけ似せて乱数を
生成し，それをもとに模
擬データを作る
• 模擬データからパラメー
タを推定する
• なんども繰り返し，模擬
パラメータ推定値を多数
もとめる
• 右図はパラメータが２つ
の場合（２次元確率分
布）
30
8th Physics Based Vision Seminar
2004/May/13
カイ２乗一定の境界による信頼限界
•
•
•
•
求まったパラメータaではχ2は最小
aから離れればχ2は増加
差分Δχ2を使う
自由度νは信頼区間を調べたいパラメータの数
31
8th Physics Based Vision Seminar
2004/May/13
特異値分解による信頼限界
• 特異値分解によってχ2当てはめをしたとき
• 特異値分解で得られる正規直交なベクトルは，Δχ2=1の楕円
体の主軸となる
• このとき，軸の長さは特異値の逆数となる
• 各軸をα倍するとΔχ2はα2倍になる
32
8th Physics Based Vision Seminar
2004/May/13
ロバストネス(頑健性, robustness)
• M推定量(M-estimators):モデ
ルの当てはめに関するロバスト
な統計的推定量
• L推定量(L-estimators):代表値
や中心的傾向の推定に使う
– メディアン
– Tukey(テューキー)の3項平均
(trimean):1/4,1/2,1/4で加重平均
• R推定量(R-estimators):順位
検定に使う
– Kolmogorov-Smirnov統計量
– Spearmanの順位相関係数
33
0
4.96
4.8
4.64
4.48
4.32
4.16
4
3.84
3.68
3.52
3.36
3.2
3.04
2.88
2.72
2.56
2.4
2.24
2.08
1.92
1.76
1.6
1.44
1.28
1.12
0.96
0.8
0.64
0.48
0.32
0.16
-5
1
0
-5
3.4
3.96
4.24
4.52
4.8
3.96
4.24
4.52
4.8
3.12
3.68
2.84
3.12
3.4
2.84
3.68
2.56
0.32
0.04
-0.2
-0.5
-0.8
-1.1
-1.4
-1.6
-1.9
-2.2
-2.5
-2.8
-3
-3.3
-3.6
-3.9
-4.2
-4.4
-4.7
2.56
0
2.28
0.1
2.28
0.2
2
0.3
1.72
0.4
1.44
0.5
2
0.6
1.72
0.7
1.44
4
1.16
0.8
1.16
0.9
0.88
Lorentz分布
0.88
1
0.6
6
0.6
0.32
0.04
-0.2
-0.5
-0.8
-1.1
-1.4
-1.6
-1.9
-2.2
-2.5
-2.8
-3
-3.3
2
重み
-3.6
正規分布
両側指数分布
Lorentz分布
-3.9
3
-4.2
5
-4.4
4.8
4.52
4.24
3.96
3.68
3.4
3.12
2.84
2.56
2.28
2
1.72
1.44
1.16
0.88
0.6
0.32
0.04
-0.2
-0.5
-0.8
-1.1
-1.4
-1.6
-1.9
-2.2
-2.5
-2.8
-3
-3.3
-3.6
-3.9
-4.2
-4.4
-4.7
-5
• 正規分布
• 両側指数分布
• Cauchy(コーシー)分布(Cauchy
distribution), Lorentz(ローレンツ)分布
(Lorentzian distribution)
-4.7
8th Physics Based Vision Seminar
2004/May/13
局所的(local)M推定量によるパラメータ推定
1
正規分布
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
1
両側指数分布
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
誤差分布
34
8th Physics Based Vision Seminar
2004/May/13
絶対偏差最小化による直線の当てはめ
• 初期値は普通に最小２乗法で解いたaとb
•
を解きたい
• 以下を反復計算すれば良い
35
池内研究室
Computer Vision Laboratory
第15章常微分方程式の数値解法
Integration of ordinary differential
equations
8th Physics Based Vision Seminar
2004/May/13
常微分方程式の問題
• ２階の微分方程式
は新しい変数zを補助的に用いて連立の１階微分方
程式に書き直せる
• 一般の常微分方程式の問題は，関数yi,i=1,2,…,N，
N元連立１階微分方程式
関数fiは既知とする
37
8th Physics Based Vision Seminar
2004/May/13
Euler(オイラー)法(Euler's method)
38
8th Physics Based Vision Seminar
2004/May/13
中点法(2次Runge-Kutta法)
39
8th Physics Based Vision Seminar
2004/May/13
Runge-Kutta(ルンゲ・クッタ)法
(Runge-Kutta method)
40
8th Physics Based Vision Seminar
2004/May/13
Runge-Kuttaに対する適応刻み幅制御
(adaptive stepsize control)
• ステップの２重化(step doubling)を十分な精度が得
られるまで繰り返す
41
8th Physics Based Vision Seminar
2004/May/13
修正中点法(modified midpoint method)
• n個のステップによって点xから点x+Hへy(x)を計算し
ていく
42
8th Physics Based Vision Seminar
2004/May/13
Bulirsch-Stoer(ブリルシュ・ストア)法
(Bulirsch-Soer method)
• Richardson(リチャードソン)の極限遅延接近(Richardson's deferred
approach to the limit):刻み幅hを減らしていったときの値を関数に当ては
めてh=0での関数値を補外する
• 有理関数補外(rational function extrapolation)を使う
• 有理関数補外をhでなくh2で行う
• xをx+Hに進める
• 修正中点法を使う
• 区間Hを分割する数n=2,4,6,8,12,16,24,32,48,64,96,…,[nj=2nj-2],…
43
8th Physics Based Vision Seminar
2004/May/13
予測子・修正子法(predictor-corrector method)
Adams-Bashforth-Moulton(アダムス・バッシュフォー
ス・モールトン)法
• 予測子ステップ(predictor step)
• 修正子ステップ(corrector step)
44
8th Physics Based Vision Seminar
2004/May/13
硬い(stiff)連立方程式
• 安定に解くために陰的(implicit)な差分を使う
• 後退(backward)Euler法
45

Download Report