ウェーブレットおきもち講座

ウェーブレットおきもち – 導入と実装 – aikiriao July 19, 2021 Available at:
Source: https://github.com/aikiriao/introduction_to_wavelet 1 / 61

もくじ 1. 準備関数はベクトル関数空間の例：フーリエ級数連続ウェーブレット変換直交ウェーブレット変換 2. 多重解像度解析 (MRA)
多重解像度解析の概要多重解像度解析の条件と性質高速ウェーブレット変換 (FWT) 3. 実装 Python による実装応用例 2 / 61

1. 準備 3 / 61

関数はベクトル関数 f : R → C を無限次元のベクトルと見よう：関数のゼロベクトルを ∀x
0(x) := 0 で定める． 4 / 61

関数の内積 2 関数 f, g の内積を次で定める： f(·), g(·) := b
a f(x)g(x)dx (1) g(x) : g(x) の複素共役積分範囲 [a, b] は対象により変わる．省略して f, g とも書く．考え方は，有限次元ベクトルと同じ： v, w := n (v)n (w)n 5 / 61

関数の内積内積に関する性質も有限次元ベクトルと同様関数の直交性：関数 f と g が直交 ⇐⇒ f,
g = 0 (2) 関数のノルム： ||f||2 2 := f, f (3) 6 / 61

関数空間関数の張るベクトル空間を関数空間という．たとえば，ノルムが有限 ||f||2 2 = ∞ −∞ |f(x)|2dx <
∞ な関数 f : R → C の集合は関数空間をなす．これを L2 (R) と書く． 7 / 61

関数空間の例：フーリエ級数周期 T = 2π/ω の関数 f のフーリエ級数：  
       f(t) = ∞ n=−∞ cn exp(jnωt) cn = 1 T T/2 −T/2 f(t) exp(−jnωt)dt (4) f(t) の exp(jnωt) による直交展開 1 √ T exp(jnωt) は周期 T の関数の正規直交基底 1 なので可能 1補足の補題 5 を参照 8 / 61

関数空間の例：フーリエ級数フーリエ係数 cn は exp との内積の定数倍 cn = 1 √
T f(·), 1 √ T exp(jnω·) (5) だから，nω の周波数成分との類似度に相当 9 / 61

連続ウェーブレット変換空間に局在する 2”ざざ波”ψ を使った信号分析手法． ψ のことをウェーブレット (Wavelet) 関数という
ψ を時間軸方向にスケール（伸び縮み）/シフトした関数を基底とする 2注意：無限範囲で値を持つものもある 10 / 61

例：ハールウェーブレットψH
ψH (t) :=    1 (0 ≤ t < 1 2 ) −1 (1 2 ≤ t < 1) 0 otherwise 11 / 61

例：メキシカンハットウェーブレット
ψ(t) := (1 − t2) exp(−t2) 12 / 61

例：シャノンウェーブレット
ψ(t) := 2sinc(2t) − sinc(t) sinc(t) := sin(πt) πt 13 / 61

何故ウェーブレットを使うのか？フーリエ級数は時間的に局在する信号の解析は苦手フーリエ級数は定常な信号（三角関数）による基底分解手法時間のズレは位相をずらすことにより表現フーリエ解析だと厳しいケースの例 1. 単一インパルスは全周波数で係数が等しくなり，どこにインパルスがあるか不明
2. 周波数が時間に応じて変化する場合（例：スイープ信号）も，周波数は解析区間で一定として解析ウェーブレットは，基底の中に時間的な局在性を入れて問題解決を図る 14 / 61

連続ウェーブレット変換 ψ に時間スケールとシフト変換を施した関数を ψa,b (t) := 1 √ a ψ
t − b a (6) と書く．a はスケール，b はシフト位置を設定． 1/ √ a はノルムを保つための定数 3 3（検算）s = (t − b)/a と置換積分すれば， ||ψa,b||2 2 = ∞ −∞ 1 a ψ t − b a ψ t − b a dt = 1 a ∞ −∞ ψ(s)ψ(s)ads = ||ψ||2 2 15 / 61

例：ハールウェーブレット
16 / 61

連続ウェーブレット変換関数 f の連続ウェーブレット変換 F(a, b) は，f と ψa,b の内積によって得られる：
F(a, b) := ∞ −∞ f(t)ψa,b (t)dt (7) = f(·), ψa,b (·) 変換先はスケール a とシフト b の 2 軸（スケログラムという）フーリエ変換の一般化（ψa,b (t) = exp(jat)）に相当 17 / 61

ウェーブレットの離散化 ψa,b （式 (6)）のスケール a とシフト b を a =
2−m, b = n2−m と整数 m, n ∈ Z で離散化して 4， ψm,n (t) := 1 √ 2−m ψ t − n2−m 2−m = 2m 2 ψ(2mt − n) (8) と書く． 4文献によっては m の符号が逆転しているので注意！ 18 / 61

直交ウェーブレット変換もし，ψm,n (t) が正規直交基底をなす ∀i, j, k, l ∈ Z.
ψi,j , ψk,l = δik δjl ならば 5，任意の信号 x ∈ L2 (R) をウェーブレット展開係数 Xm,n = x, ψm,n を用いて x(t) = m,n Xm,n ψm,n (t) (9) と展開できる．これを直交ウェーブレット変換という． 5スケール/シフトの両方で正規直交 19 / 61

2. 多重解像度解析(MRA) 20 / 61

多重解像度解析(MRA)の概要多重解像度解析 (MRA6) は直交ウェーブレット変換の手法のひとつ．スケール M の信号 f(M)(t) を
1 つスケールの落とした信号 f(M−1)(t) と誤差信号 g(M−1)(t) の和で表すことを考える： f(M)(t) = f(M−1)(t) + g(M−1)(t) (10) M はスケールの細かさ．大きくとればどんな信号でも近似できる． 6MultiResolution Analysis 21 / 61

多重解像度解析(MRA)の概要式 (10) を繰り返し適用していくと， f(M)(t) = f(M−1)(t) + g(M−1)(t) =
f(M−2)(t) + g(M−1)(t) + g(M−2)(t) = ... = f(J)(t) + M−1 m=J g(m)(t) (11) 任意のスケール J < M まで信号を分解できる！ 22 / 61

多重解像度解析(MRA)の概要この操作を形式的に考えていく． 23 / 61

多重解像度解析(MRA) MRA は次の部分空間 Vm と関数 φ Vm := n cn
φm,n (t) cn ∈ l2 (Z) (12) φm,n (t) := 2m 2 φ(2mt − n) (13) を基に構成される 78．φ をスケーリング関数 (scaling function) という． 7l2(Z) は二乗総和可能な数列 ( n |cn |2 < ∞) の集合 8文献によっては m の符号が逆転しているので注意！ 24 / 61

MRAの満たすべき条件 I MRA は，集合 Vm に対して以下の条件 (M1)-(M4) を要求する 9： (M1)
V0 は正規直交基底 {φ(t − n)|n ∈ Z}, φ ∈ L2 (R) によって張られる 25 / 61

MRAの満たすべき条件 II M1 お気持ちスケーリング関数 φ は整数シフトに関して正規直交基底となることを要求： φ(· −
n), φ(· − m) = δnm (M1) が成立すれば，任意の m, p, q ∈ Z に対し， φm,p, φm,q = ∞ −∞ 2m 2 φ(2mt − p)2m 2 φ(2mt − q)dt = 2m ∞ −∞ φ(2mt − p)φ(2mt − q)dt = ∞ −∞ φ(s − p)φ(s − q)ds (s = 2mt) = δpq （∵ φ はシフトに関して正規直交）だから，V0 さえ定まれば全ての Vm が張られる． 26 / 61

MRAの満たすべき条件 III (M2) Vm ⊂ Vm+1 M2 お気持ち Vm は入れ子構造をなす．
スケール m を上げると表現可能な信号が増える 27 / 61

MRAの満たすべき条件 IV (M3) m Vm = L2 (R)（A：集合 A の閉包）
M3 お気持ち lim m→∞ Vm = L2(R) であること，つまり，L2(R) のどんな信号であっても任意の精度で近似できる (M4) m Vm = {0} M4 お気持ち lim m→−∞ Vm = {0} であること，つまり，スケールを極限まで落とすと Vm に含まれる信号は定値関数 0 のみになる 9なぜこの要求 (M1)-(M4) で良いかは割愛．[1, 2] を参照のこと 28 / 61

例：ハールウェーブレットψH ψH に対応するスケーリング関数 φH φH (t) := 1 (0 ≤
t < 1) 0 otherwise φHm,n (t) := 2m 2 φH (2mt − n) とかく． 29 / 61

例：ハールウェーブレットψH φH が条件 (M1)-(M4) を満たすかチェック (M1) a, b ∈
Z に対して， φH (· − a), φH (· − b) = ∞ −∞ φH (t − a)φH (t − b)dt = ∞ −∞ φH (s)φH (s + a − b)ds (s = t − a) = δab なので {φH (t − n)|n ∈ Z} は正規直交基底． 30 / 61

例：ハールウェーブレットψH (M2) 任意の m, n ∈ Z に対して， φHm,n (t)
= 2m 2 2−mn ≤ t < 2−m(n + 1) 0 otherwise φHm+1,2n (t) = 2m+1 2 2−(m+1)2n ≤ t < 2−(m+1)(2n + 1) 0 otherwise = 2m+1 2 2−mn ≤ t < 2−m(n + 1 2 ) 0 otherwise だから，任意の φHm,n ∈ Vm は φHm+1,2n , φHm+1,2n+1 ∈ Vm+1 を用いて φHm,n (t) = 2− 1 2 φHm+1,2n (t) + 2− 1 2 φHm+1,2n+1(t) (14) と表せる．よって Vm ⊂ Vm+1 ． 31 / 61

例：ハールウェーブレットψH（説明は省略） (M3) lim m→∞ φHm,0(t) = δ(t)（δ(t)：インパルス関数）の観察から， lim m→∞
Vm はインパルス関数を時間シフトした集合となる．任意の f ∈ L2(R)，s ∈ R に対して δ(t − s) ∈ lim m→∞ Vm が存在して， f(·), δ(· − s) = ∞ −∞ f(t)δ(t − s)dt = f(s) だから，f ∈ lim m→∞ Vm で L2(R) ⊂ lim m→∞ Vm ．一方， φHm,n ∈ L2(R) だから L2(R) ⊃ lim m→∞ Vm ．従って， lim m→∞ Vm = L2(R)． (M4) lim m→−∞ φHm,n (t) = 0 より確かに lim m→−∞ Vm は 0 しか含まない． 32 / 61

Dilation方程式（ツースケール関係） (M2) が成立すれば V−1 ⊂ V0 ．よって φ(t/2) ∈ V−1
は V0 の基底 {φ(t − n)|n ∈ Z} の線形結合で表現でき， φ(t/2) = √ 2 n h[n]φ(t − n) (15) を満たす係数 {h[n]} が存在する 10．これを Dilation 方程式とかツースケール関係と呼ぶ． 10 √ 2 は規格化定数 33 / 61

MRAのウェーブレット関数 (M2) より Vm ⊂ Vm+1 だから，Vm ⊕ Wm =
Vm+1 （⊕：直和）となる Wm が存在する： 34 / 61

MRAのウェーブレット関数 Wm = Vm+1 Vm より，Wm はスケール m のウェーブレットから張られるように定める：
Wm := n cn ψm,n (t) cn ∈ l2 (Z) (16) このとき，Wm は Vm ⊥ Wm 11 を満さなければならない．その成立条件を見ていく． 11Vm, Wm それぞれから選んだ任意の元（関数）が直交している 35 / 61

MRAのウェーブレット関数 W−1 ⊂ V0 だから，ψ(t/2) ∈ W−1 も V0 の基底
{φ(t − n)|n ∈ Z} の線形結合で表せて， ψ(t/2) = √ 2 n g[n]φ(t − n) (17) を満たす係数 {g[n]} が存在する． 36 / 61

MRAのウェーブレット関数このとき， g[n] = (−1)nh[1 − n] (18) とすれば Vm
⊥ Wm が満たされる 12． ⇒ 必要なのは {h[n]} だけ！有名なウェーブレットは {h[n]} の数表が与えられている！ 12証明は長くなる．補足にて述べる（系 1） 37 / 61

例：ハールウェーブレット式 (14) において m = −1, n = 0
とおくと， φH−1,0 (t) = 2−1 2 φH0,0 (t) + 2−1 2 φH0,1 (t) ⇐⇒ 2−1 2 φH (t/2) = 2−1 2 φH (t) + 2−1 2 φH (t − 1) ⇐⇒ φH (t/2) = φH (t) + φH (t − 1) Dilation 方程式（式 (15)）と見比べると h[0] = h[1] = 2−1 2 式 (18) により g[0] = 2−1 2 , g[1] = −2−1 2 が得られ， ψH (t/2) = φH (t) − φH (t − 1) φH から ψH が得られた． 38 / 61

MRAの直交ウェーブレット変換 Vm ⊕ Wm = Vm+1 を繰り返し適用すると， Vm = Vm−1
⊕ Wm−1 = Vm−2 ⊕ Wm−2 ⊕ Wm−1 = ... = VJ ⊕ WJ ⊕ WJ−1 ⊕ · · · ⊕ Wm−1 (J < m) = VJ ⊕ m−1 k=J Wk と書ける．m → ∞ とすると， lim m→∞ Vm = VJ ⊕ ∞ k=J Wk (19) 39 / 61

MRAの直交ウェーブレット変換任意の x(t) ∈ L2 (R) は以下のように直交展開できる 13： x(t)
= n pJ [n]φJ,n (t) ∈VJ + ∞ m=J n qm [n]ψm,n (t) ∈Wm (20) ここで，pm [n] = x, φm,n , qm [n] = x, ψm,n ． 13∵ 式 (19)，(M3)（limm→∞ Vm = L2(R)），Vm ⊥ Wm (m ∈ Z) より 40 / 61

MRAの直交ウェーブレット変換スケール M > J 以上の成分を無視すると以下の近似が成り立つ： x(t) ≈ n
pJ [n]φJ,n (t) + M−1 m=J n qm [n]ψm,n (t) (21) ここで，f(m)(t) := n pm [n]φm,n (t)， g(m)(t) := n qm [n]ψm,n (t) とすると x(t) ≈ f(J)(t) + M−1 m=J g(m)(t) 式 (11)そのもの：多重解像度解析を構成している． 41 / 61

FWTの導出（展開）式 (20) の展開係数 pm [n], qm [n] は h[n],
g[n] を用いて計算できる．補題 1 (展開計算) pm [n] = k h[k − 2n]pm+1 [k] (22) qm [n] = k g[k − 2n]pm+1 [k] (23) 畳み込み演算と間引きになっていることに注目！ 42 / 61

FWTの導出（展開）（証明）準備として，Dilation 方程式（式 (15)）から導かれる関係式に着目する． φm,n (t) = 2 m 2
φ(2mt − n) = 2 m 2 φ 2m+1 t − 2n 2 = 2 m 2 √ 2 k h[k]φ(2m+1 t − 2n − k) （∵ Dilation 方程式） = k h[k]2 m+1 2 φ(2m+1 − 2n − k) = k h[k]φm+1,2n+k (t) (24) ψm,n に関しても同様にして，式 (17) から， ψm,n (t) = k g[k]φm+1,2n+k (t) (25) 任意の x(t) ∈ L2 (R) をとる．式 (24)，式 (25) を用いると， pm [n] = x, φm,n = x, k h[k]φm+1,2n+k = k h[k] x, φm+1,2n+k = k h[k]pm+1 [2n + k] = k h[k − 2n]pm+1 [k] (k ← 2n + k) qm [n] = x, ψm,n = x, k g[k]φm+1,2n+k = k g[k − 2n]pm+1 [k] 43 / 61

FWTの導出（再構成）逆に，pm [n], qm [n] から pm+1 [n] を再構成できる．補題
2 (再構成計算) pm+1 [n] = k h[n − 2k]pm [k] + k g[n − 2k]qm [k] (26) 式 (22)，式 (23)，式 (26) による逐次的な展開・再構成アルゴリズムを高速ウェーブレット変換 (FWT, Fast Wavelet Transform) という． 44 / 61

FWTの導出（再構成）（証明）pm+1 [n] の定義と x の直交展開式（式 (20)）より， pm+1 [n] =
x, φm+1 [n] = k pm [k]φm,k + ∞ t=m k qt [k]ψt,k, φm+1,n = k pm [k] φm,k, φm+1,n + ∞ t=m k qt [k] ψt,k, φm+1,n ここで，式 (24) より， φm,k, φm+1,n = l h[l]φm+1,2k+l, φm+1,n = s h[s − 2k]φm+1,s, φm+1,n (s = 2k + l) = s h[s − 2k] φm+1,s, φm+1,n = s h[s − 2k]δsn （∵ φ のシフトに関する直交性） = h[n − 2k] 45 / 61

FWTの導出（再構成）任意の t ≥ m + 1 に対して (M2) より
Vm+1 ⊂ Vt だから， Wt ⊥ Vt =⇒ Wt ⊥ Vm+1 よって ψt,k, φm+1,n = 0 (t ≥ m + 1)．これより， ∞ t=m ψt,k, φm+1,n = ψm,k, φm+1,n = l g[l]φm+1,2k+l, φm+1,n （∵ 式 (25)） = s g[s − 2k]φm+1,s, φm+1,n (s = 2k + l) = s g[s − 2k] φm+1,s, φm+1,n = s g[s − 2k]δsn （∵ φ のシフトに関する直交性） = g[n − 2k] 以上の結果から， pm+1 [n] = k pm [k] φm,k, φm+1,n + ∞ t=m k qt [k] ψt,k, φm+1,n = k h[n − 2k]pm [k] + k g[n − 2k]qm [k] 46 / 61

FWTのフィルタ表現展開処理（式 (22)，(23)）をフィルタ表現すると H G ↓ 2 qm ↓ 2
H G ↓ 2 qm−1 ↓ 2 G ↓ 2 qm−2 pm+1 pm pm−1 Figure 1: 展開処理．↓ 2 は 2 倍の間引き（1 個飛ばし） H はローパスフィルタ，G はハイパスフィルタ pm はスケールを落とした信号成分，qm は差分成分だから 47 / 61

FWTのフィルタ表現再構成処理（式 (26)）をフィルタ表現すると ↑ 2 H G ↑ 2 qm−1
↑ 2 H G ↑ 2 qm pm−1 pm pm+1 Figure 2: 再構成処理．↑ 2 は 2 倍の補間（0 値挿入） 48 / 61

FWTの計算量 FWT の時間的計算量は，入力データ個数 N に対して O(N)．（証明）pM [n] が
N 個のデータからなり，H, G を長さ L の FIR フィルタとすると，各スケールにおける積和演算回数は次のようになる： 1 段目 pM [n] は N 個のデータからなるため，pM−1 [n] の計算で NL 回の積和演算， qM−1 [n] の計算で NL 回の積和演算 ⇒ 計 2NL 回 2 段目 pM−1 [n] は N/2 個のデータからなるため，pM−2 [n] の計算で NL/2 回の積和演算，qM−2 [n] の計算で NL/2 回の積和演算 ⇒ 計 NL 回 3 段目 pM−2 [n] は N/4 個のデータからなるため，pM−3 [n] の計算で NL/4 回の積和演算，qM−3 [n] の計算で NL/4 回の積和演算 ⇒ 計 NL/2 回 ... よってすべての段の積和演算回数の和は， 2NL + NL + NL 2 + ... ≤ 4NL 再構成処理においても逆向きの手順の計算を行うだけなので計算量は変わらない．従って，計算量は O(N)． 49 / 61

展開係数の入力による近似補題 3 (係数の近似 [3]) M は十分に大きいとする．x が連続かつ φ がコン
パクトサポート a ならば， pM [n] ≈ Ax(2−M n) A = ∞ −∞ φ(t)dt (27) a有界な範囲外では常に 0 となっていること ⇒ 入力信号系列をそのまま分解・再構成に突っ込んで（ほぼ）OK！14 14ハールウェーブレットなら A = 1 50 / 61

展開係数の入力による近似（証明）φ は仮定より，ある範囲 [−L, L] の外では常に 0 とする．このとき pM [n]
は， pM [n] = x, φM,n = ∞ −∞ x(s)φM,n (s)ds = ∞ −∞ x(s)2M φ(2M s − n)ds = ∞ −∞ x(2−M t + 2−M n)2M φ(t)2−M dt (t = 2M s − n) = ∞ −∞ x(2−M t + 2−M n)φ(t)dt = L −L x(2−M t + 2−M n)φ(t)dt （∵ φ はコンパクトサポート） M が十分大きいとき，x の連続性により区間 [−2−M L + 2−M n, 2−M L + 2−M n] で x はほぼ一定とみなせるから， x(2−M t + 2−M n) ≈ x(2−M n) この近似により， pM [n] ≈ L −L x(2−M n)φ(t)dt = x(2−M n) L −L φ(t)dt = x(2−M n) ∞ −∞ φ(t)dt （∵ φ はコンパクトサポート） = Ax(2−M n) 51 / 61

FWTの特徴特徴式 (18)（g[n] := (−1)nh[1 − n]）より実行に必要なのは h[n]
のみ！！有名なウェーブレットは h[n] が公開されている計算量は O(N)！！実用上は入力信号をそのまま使用できる！！注意点スケールは 2 の指数に限る（FFT は等間隔スケール）不確定性原理（補足参照）には縛られる 52 / 61

3. 実装 53 / 61

ウェーブレット係数の計算(fwt.py) 28 def calculate_wavelet_coef(scaling_coefficients): 29 """ スケーリング係数からウェーブレット係数を生成 """ 30 wavelet_coefficients
= scaling_coefficients [::-1].copy() # 順序逆の配列 31 for i in range(len(scaling_coefficients)): 32 wavelet_coefficients[i] *= ((-1) ** i) 33 return wavelet_coefficients フィルタ係数条件（式 (18)）において， g[n] = (−1)nh[1 − n] = (−1)nh[L + 1 − n] だから（L：フィルタ係数長．離散フーリエ変換の性質より周期性を持つ），n の範囲を 0, ..., L − 1 に直すと， g[n] = (−1)nh[L − 1 − n] (n = 0, ..., L − 1) 54 / 61

FWTの実装例(fwt.py) 36 def fwt1d(src, scaling_coef): 37 """ 1次元高速ウェーブレット変換 """ 38
# ウェーブレット係数計算 39 wavelet_coef = calculate_wavelet_coef( scaling_coef) 40 # 入力が整数だと丸め込まれるためfloatに変換 41 src = src.astype(float) 42 # correlate1dはフィルタカーネルの半分（中心）だけ出力が後ろにずれるので 43 # 先に入力を前にずらしておく 44 src = np.roll(src, -len(scaling_coef) // 2) 45 # 畳み込み入力の端点は巡回 46 # フィルタのインデックスが正方向に増加するため correlate1dを使用 47 decomp_src = correlate1d(src, scaling_coef, mode= ’wrap’)[::2] 48 decomp_wav = correlate1d(src, wavelet_coef, mode= ’wrap’)[::2] 49 return [decomp_src, decomp_wav] 55 / 61

IFWTの実装例(fwt.py) 52 def ifwt1d(decomp_src, decomp_wav, scaling_coef): 53 """ 1次元高速ウェーブレット逆変換 """
54 # ウェーブレット係数計算 55 wavelet_coef = calculate_wavelet_coef( scaling_coef) 56 src_len = 2 * len(decomp_src) 57 # 0値挿入 58 scaling_interp = np.zeros(src_len) 59 scaling_interp[::2] = decomp_src 60 wavelet_interp = np.zeros(src_len) 61 wavelet_interp[::2] = decomp_wav 62 # 畳み込み入力の端点は巡回 63 src = convolve1d(scaling_interp, scaling_coef, mode=’wrap’) 64 src += convolve1d(wavelet_interp, wavelet_coef, mode=’wrap’) 65 # convolve1dはフィルタカーネルの半分（中心）だけ出力が前にずれるので 66 # 入力を後ろにずらす 67 src = np.roll(src, len(scaling_coef) // 2) 68 return src 56 / 61

2次元FWTの考え方縦/横方向に FWT を適用すれば OK． Figure 3: 結果を横/縦の成分で示すと， 1 =低域/低域，
2 =高域/低域， 3 =低域/高域， 4 =高域/高域 57 / 61

2次元FWTの実装例(fwt.py) 71 def fwt2d(src2d, scaling_coef): 72 """ 2次元高速ウェーブレット変換 """ 81
# src2dを低域（左）と高域（右）に分解 82 for j in range(src_len): 83 src_l, src_h = fwt1d(src2d[j, :], scaling_coef) 84 src2d_l[j, :] = src_l 85 src2d_h[j, :] = src_h 86 # src2d_l, src2d_hを更に左上(ll)、左下(hl)、右上(lh )、右下(hh)に分解 87 for j in range(half_src_len): 88 src_l, src_h = fwt1d(src2d_l[:, j], scaling_coef) 89 src2d_ll[:, j] = src_l 90 src2d_hl[:, j] = src_h 91 src_l, src_h = fwt1d(src2d_h[:, j], scaling_coef) 92 src2d_lh[:, j] = src_l 93 src2d_hh[:, j] = src_h 94 return [src2d_ll, src2d_hl, src2d_lh, src2d_hh] 58 / 61

2次元IFWTの実装例(fwt.py) 97 def ifwt2d(src2d_ll, src2d_hl, src2d_lh, src2d_hh, scaling_coef): 98 """
2次元高速ウェーブレット逆変換 """ 99 src_len = src2d_ll.shape[0] 100 twice_src_len = 2 * src_len 101 src2d = np.zeros((twice_src_len, twice_src_len)) 102 src2d_l = np.zeros((twice_src_len, src_len)) 103 src2d_h = np.zeros((twice_src_len, src_len)) 104 # 左上(ll)、左下(hl)、右上(lh)、右下(hh)から左(l)、右(h)に合成 105 for j in range(src_len): 106 src2d_l[:, j] = ifwt1d(src2d_ll[:, j], src2d_hl[:, j], scaling_coef) 107 src2d_h[:, j] = ifwt1d(src2d_lh[:, j], src2d_hh[:, j], scaling_coef) 108 # 左(l)、右(h)から元を合成 109 for j in range(twice_src_len): 110 src2d[j, :] = ifwt1d(src2d_l[j, :], src2d_h[j , :], scaling_coef) 111 return src2d 59 / 61

静止画の分解低域/低域の結果に繰り返し 2 次元 FWT を適用することでピラミッドを構築できる． 60 / 61

他の応用 NoteBook15 参照．以下のコンテンツを予定．静止画の分解・再構成ノイズ（欠損）除去圧縮（エントロピー観察） 15https://github.com/aikiriao/introduction_to_wavelet/ blob/main/implementation/fwt_demos.ipynb 61 /
61

A. 三角関数の完全性 1 / 59

素朴な疑問周期 T の関数 f(x) が以下のようにフーリエ級数展開できるとしよう． f(x) = n
cn exp(jnωx) cn = 1 T T 2 −T 2 f(x) exp(−jnωx)dx いや待てよ：三角関数の和で fって表現できるの？できたとして，f にはどんな条件が必要？ 2 / 59

ディリクレ核一旦フーリエ級数展開ができると受け入れて，フーリエ級数の式から cn を消して整理してみる： f(x) = n cn exp(jnωx)
= n 1 T T 2 − T 2 f(s) exp(−jnωs)ds exp(jnωx) = 1 T T 2 − T 2 f(s) n exp[−jnω(x − s)]ds = 1 T x+ T 2 x− T 2 f(x − u) n exp(jnωu)du (u = x − s) = 1 T T 2 − T 2 f(x − u) n exp(jnωu)du （∵ f は周期 T） 3 / 59

ディリクレ核ここで，以下のディリクレ核 (Dirichlet Kernel)DN を定義する： DN (x) = N n=−N
exp(jnx) (28) f は DN との畳込みで復元できる： f(x) = 1 T T 2 −T 2 f(x − u)D∞ (ωu)du DN の性質を観察してみよう 4 / 59

ディリクレ核 DN は和を使わない式でも表現できる： DN (x) = sin 1 2 +
N x sin 1 2 x (29) （証明）まず， DN (x) = N n=−N exp(jnx) = exp(0) + N n=1 {exp(jnx) + exp(−jnx)} = 1 + N n=1 {cos(nx) + j sin(nx) + cos(nx) − j sin(nx)} = 1 + 2 N n=1 cos(nx) (30) より，DN は cos の和で表現できる． 5 / 59

ディリクレ核両辺 sin 1 2 x を乗じると， sin 1 2
x DN (x) = 2 sin 1 2 x 1 + 2 N n=1 cos(nx) = sin 1 2 x + 2 sin 1 2 x {cos(x) + cos(2x) + cos(3x) + ... + cos(Nx)} = sin 1 2 x + sin 1 2 − 1 x + sin 1 2 + 1 x + sin 1 2 − 2 x + sin 1 2 + 2 x + ... + sin 1 2 − N x + sin 1 2 + N x （∵ 加法定理 sin(nx) cos(mx) = 1 2 {sin(n + m)x + sin(n − m)x}） = sin 1 2 x + sin − 1 2 x + sin 3 2 x + sin − 3 2 x + ... + sin 1 2 + N x = sin 1 2 + N x よって， DN (x) = sin 1 2 + N x sin 1 2 x 6 / 59

ディリクレ核 DN (ωx) の 1 周期分に渡る積分は T に等しい： T 2
−T 2 DN (ωx)dx = T (31) （証明）式 (30) を使うと， T 2 − T 2 DN (ωx) = T 2 − T 2 1 + 2 N n=1 cos(nωx) dx = T 2 − T 2 dx + 2 N n=1 T 2 − T 2 cos(nωx)dx = T + 2 N n=1 1 nω sin n 2π T x T 2 − T 2 ∵ ω = 2π T = T + 2 N n=1 1 nω {sin(nπ) − sin(−nπ)} = T 7 / 59

フーリエ級数による近似ディリクレ核を使うことで，f の N 次までのフーリエ級数 SN (x) を SN
(x) = 1 T T 2 −T 2 f(x − u)DN (ωu)du (32) と表現できる．では， lim N→∞ |SN (x) − f(x)| = 0 (33) となるのはどんな時か？ 8 / 59

フーリエ級数による近似 SN (x) − f(x) を変形していくと， SN (x) − f(x)
= 1 T T 2 − T 2 f(x − u)DN (ωu)du − f(x) = 1 T T 2 − T 2 f(x − u)DN (ωu)du − 1 T T 2 − T 2 DN (ωu)f(x)du （∵ 式 (31)） = 1 T T 2 − T 2 {f(x − u) − f(x)} DN (ωu)du = 1 T T 2 − T 2 f(x − u) − f(x) sin 1 2 ωu sin 1 2 + N ωu du （∵ 式 (29)）(34) 9 / 59

フーリエ級数による近似式 (34) が 0 に収束するか否かは，次の補題により判定できる：補題 4 (リーマン・ルベーグの補題)
範囲 [a, b] で区分的に連続な関数 f に対して， lim N→∞ b a f(x) sin(Nx)dx = 0 (35) 補題の証明は後ろで行う．これより，f(x−u)−f(x) sin(1 2 ωu) が区分的に連続であれば良い． 10 / 59

フーリエ級数による近似定理 1 (収束定理) f が点 x で微分可能ならば，フーリエ級数による近似（式 (32)）は
f(x) に点収束する．（証明） f(x−u)−f(x) sin( 1 2 ωu ) の分母が 0 となる時に収束すれば良い．u → 0 の極限をとると， lim u→0 f(x − u) − f(x) sin 1 2 ωu = lim u→0 f(x − u) − f(x) u u sin 1 2 ωu だから，この極限が収束するためには， lim u→0 f(x − u) − f(x) u , lim u→0 u sin 1 2 ωu が共に収束することが必要である．第二項はロピタルの定理より， lim u→0 u sin 1 2 ωu = lim u→0 1 1 2 ω cos 1 2 ωu = 2 ω だから収束する．よって，第一項の f の点 x における微分が収束すれば良いことが分かる．そしてその時，リーマン・ルベーグの補題の仮定が満たされ，収束定理が成立する． 11 / 59

三角関数の完全性補題 5 (三角関数の完全性) 1 √ T exp(jnωt)|n ∈ Z
はフーリエ級数展開可能な関数の正規直交基底系（証明）収束定理より任意のフーリエ級数展開可能な関数を表現可能だから，基底系をなすことは示されている．次に直交性を示す．n, m ∈ Z に対して 1 周期分の内積をとると， 1 √ T exp(jnω·), 1 √ T exp(jmω·) = 1 T T/2 −T/2 exp(jnωt) exp(−jmωt)dt = 1 T exp {j(n − m)2πt/T} j(n − m)2π/T T/2 −T/2 (∵ ω = 2π/T) = 1 π exp {j(n − m)π} − exp {−j(n − m)π} j2(n − m) = 1 π sin {(n − m)π} n − m （∵ オイラーの公式） = δnm （δ：クロネッカーのデルタ）最後の式変形で n = m のときは x = n − m とおいて x → 0 とする： lim x→0 sin(πx) x = lim x→0 π cos(πx) 1 = π 12 / 59

三角関数の完全性最後に，これ以外の基底はないことを背理法により示す．「{exp(jnωt)|n ∈ Z} の全てに直交する連続関数基底 h(t) が存在する」と仮定する．仮定により h(·),
exp(jnω·) = 0 が成立するが，h(t) = 0 は空間を生成しないゼロベクトルのため，h(t) = 0 で考える． h(t) のフーリエ級数展開を考える．n 次のフーリエ係数は， cn = 1 T T 2 − T 2 h(t) exp(−jnωt)dt = 1 T h(·), exp(jnω·) = 0 （∵ h(t) の仮定）だから，h(t) の N 次までのフーリエ級数近似 SN (t) は， SN (t) = N n=−N cn exp(jnωt) = 0 となる．収束定理を用いると， lim N→∞ SN (t) = h(t) = 0 これは，h(t) = 0 で考えていることに矛盾する．従って {exp(jnωt)|n ∈ Z} 以外に基底は存在せず，定理が成立する． 13 / 59

リーマン・ルベーグの補題の証明区間 [a, b] を M 分割し，区間の小さい方から順に a = x1,
x2, ..., xM−1, xM = b となるようにとる．このとき， b a f(x) sin(Nx)dx = M−1 k=1 xk+1 xk f(x) sin(Nx)dx ≤ M−1 k=1 xk+1 xk f(x) sin(Nx)dx 区間 [xk, xk+1 ] の積分に注目すると， xk+1 xk f(x) sin(Nx)dx = xk+1 xk {f(x) − f(xk ) + f(xk )} sin(Nx)dx = xk+1 xk {f(x) − f(xk )} sin(Nx)dx + xk+1 xk f(xk ) sin(Nx)dx ≤ xk+1 xk {f(x) − f(xk )} sin(Nx)dx + xk+1 xk f(xk ) sin(Nx)dx 14 / 59

リーマン・ルベーグの補題の証明前ページ最後の式の 1 項目は， xk+1 xk {f(x) − f(xk )}
sin(Nx)dx ≤ xk+1 xk | {f(x) − f(xk )} sin(Nx)|dx ≤ xk+1 xk |f(x) − f(xk )|| sin(Nx)|dx ≤ xk+1 xk |f(x) − f(xk )|dx (∵ | sin(x)| ≤ 1) 2 項目については，連続な関数は有界だから |f(x)| < L (a ≤ x ≤ b) を満たす定数 L が存在する仮定のもと， xk+1 xk f(xk ) sin(Nx)dx ≤ max x∈[xk,xk+1] |f(x)| xk+1 xk sin(Nx)dx ≤ max x∈[xk,xk+1] |f(x)| − 1 N cos(Nx) xk+1 xk < L − 1 N [cos(Nx)]xk+1 xk ≤ 2L N 15 / 59

リーマン・ルベーグの補題の証明元に戻って結果をまとめると， b a f(x) sin(Nx)dx ≤ M−1 k=1 xk+1
xk f(x) sin(Nx)dx < M−1 k=1 xk+1 xk |f(x) − f(xk )|dx + 2L N = b a |f(x) − f(xk )|dx + 2LM N 任意の正数 ε > 0 に対して，十分に大きな M をとって分割を細かくすれば， |f(x) − f(xk )| < ε 2(b−a) とできる．またこのとき，N も大きく取って 2LM N < ε 2 とすることができる．よって， b a f(x) sin(Nx)dx < b a ε 2(b − a) dx + ε 2 = ε 2 + ε 2 = ε ε は任意だったから， lim N→∞ b a f(x) sin(Nx)dx = 0 16 / 59

まとめ三角関数の和で fって表現できるの？周期関数であれば表現できる．三角関数は周期関数の直交基底．できたとして，f にはどんな条件が必要？各点で微分可能であることが必要．不連続関数でも強引に当てはめられるけど，不連続点で誤差が増大する（ギブス現象）
17 / 59

B. フィルタ係数条件の証明 18 / 59

式(18)の証明補題6 補題 6 (周波数領域での直交条件) f, g ∈ L2 (R)
のフーリエ変換を F, G とする． 1. 集合 {f(t − n)|n ∈ Z} と {g(t − n)|n ∈ Z} が直交する必要十分条件は k F(ω + 2πk)G(ω + 2πk) = 0 (36) 2. 集合 {f(t − n)|n ∈ Z} が正規直交系となる必要十分条件は k |F(ω + 2πk)|2 = 1 (37) 19 / 59

式(18)の証明補題6 （証明）1. は，任意の n, m ∈ Z に対して， f(·
− n), g(· − m) = ∞ −∞ f(t − n)g(t − m)dt = 1 2π ∞ −∞ F [f(t − n)] F [g(t − m)]dω （∵ パーセバルの等式） = 1 2π ∞ −∞ F(ω) exp(−jnω)G(ω) exp(−jmω)dω （∵ 時間シフトのフーリエ変換） = 1 2π ∞ −∞ F(ω)G(ω) exp[−j(n − m)ω]dω = 1 2π k π −π F(ω + 2πk)G(ω + 2πk) exp[−j(n − m)ω]dω = 1 2π π −π k F(ω + 2πk)G(ω + 2πk) exp[−j(n − m)ω]dω = 0 が成立しているため，直交するときの必要十分条件は {·} の中が 0 になることである．よって，1. が示された． 20 / 59

式(18)の証明補題6 2. は，1. の証明において f = g とすれば， f(·
− n), f(· − m) = ∞ −∞ f(t − n)f(t − m)dt = 1 2π ∞ −∞ F(ω)F(ω) exp[−j(n − m)ω]dω = 1 2π π −π k |F(ω + 2πk)|2 exp[−j(n − m)ω]dω = δnm となり，この等式は k |F(ω + 2πk)|2 = 1 のとき，かつそのときに限り成立する．実際， f(· − n), f(· − m) = 1 2π π −π exp[−j(n − m)ω]dω = 1 2π − exp[−j(n − m)ω] j(n − m) π −π = 1 2π − exp[−j(n − m)ω] j(n − m) π −π = 1 π sin {(n − m)π} n − m から，n = m のときと n = m のときで場合分けすれば確認できる（n = m のときは x = n − m として x → 0 を考える）．よって，2. は示された． 21 / 59

式(18)の証明補題7 補題 7 (Dilation 方程式を満たす H(ω) の性質) Dilation 方程式（式
(15)）を満たす係数 h[n] の離散時間フーリエ変換 H(ω)a は， |H(ω)|2 + |H(ω + π)|2 = 2 (38) を満たす． aH(ω) := n h[n] exp(−jnω) 22 / 59

式(18)の証明補題7 （証明）補題 6 より，φ(t) のフーリエ変換を Φ(ω) と書くと，{φ(t − n)|n
∈ Z} が正規直交基底をなすための必要十分条件は， k |Φ(ω + 2πk)|2 = 1 (39) と書ける．これに，Dilation 方程式の両辺をフーリエ変換した F [φ(t/2)] = √ 2F n h[n]φ(t − n) ⇐⇒ 2Φ(2ω) = √ 2 ∞ −∞ n h[n]φ(t − n) exp(−jωt)dt = √ 2 n h[n] ∞ −∞ φ(t − n) exp(−jωt)dt = √ 2 n h[n]Φ(ω) exp(−jnω) = √ 2Φ(ω)H(ω) ⇐⇒ √ 2Φ(2ω) = Φ(ω)H(ω) (40) を導入することを考える． 23 / 59

式(18)の証明補題7 式 (39)，式 (40) より， 2 = 2 k
|Φ(2ω + 2πk)|2 （式 (39) で ω → 2ω とした） = k |Φ(ω + πk)|2 |H(ω + πk)|2 （∵ 式 (40)） = k |Φ(ω + 2πk)|2 |H(ω + 2πk)|2 + k |Φ(ω + 2πk + π)|2 |H(ω + 2πk + π)|2 = |H(ω)|2 k |Φ(ω + 2πk)|2 + |H(ω + π)|2 k |Φ(ω + 2πk + π)|2 （∵ H(ω) は周期 2π16） = |H(ω)|2 + |H(ω + π)|2 （∵ 第 2 項は式 (39) で ω → ω + π とする）よって補題 7 の成立が確かめられた． 16∵ H(ω + 2π) = n h[n] exp[−jn(ω + 2π)] = n h[n] exp(−jnω) = H(ω) 24 / 59

式(18)の証明補題8 補題 8 (H(ω), G(ω) による直交条件) {φ(t − n)|n
∈ Z} は正規直交系とする． 1. 式 (17) を満たす ψ を用いた {ψ(t − n)|n ∈ Z} が正規直交系となる必要十分条件は， |G(ω)|2 + |G(ω + π)|2 = 2 (41) 2. {φ(t − n)|n ∈ Z} と {ψ(t − n)|n ∈ Z} が直交する必要十分条件は， H(ω)G(ω) + H(ω + π)G(ω + π) = 2 (42) 25 / 59

式(18)の証明補題8 （証明）ウェーブレット関数 ψ(t) のフーリエ変換を Ψ(ω) と書く．このとき，式 (17) の両辺をフーリエ変換すると，式
(40) と同様に， √ 2Ψ(2ω) = G(ω)Φ(ω) (43) が成り立つ．式 (41) の証明は補題 7 と全く同様であり，Dilation 方程式の代わりに式 (17) を用いれば良い．式 (42) は，補題 6 の 1. を φ, ψ に適用すると， 0 = k Φ(2ω + 2πk)Ψ(2ω + 2πk) = k 1 2 H(ω + πk)Φ(ω + πk)G(ω + πk)Φ(ω + πk) （∵ 式 (40)，式 (43)） = 1 2 k H(ω + πk)G(ω + πk)|Φ(ω + πk)|2 = 1 2 k H(ω + 2πk)G(ω + 2πk)|Φ(ω + 2πk)|2 + k H(ω + 2πk + π)G(ω + 2πk + π)|Φ(ω + 2πk + π)|2 26 / 59

式(18)の証明補題8 （前ページの続き） 0 = 1 2 H(ω)G(ω) k |Φ(ω
+ 2πk)|2 + H(ω + π)G(ω + π) k |Φ(ω + 2πk + π)|2 （∵ H(ω), G(ω) は周期 2π） = 1 2 H(ω)G(ω) + H(ω + π)G(ω + π) （∵ 式 (39)）により示される． 27 / 59

式(18)の証明補題9 h, g の離散時間フーリエ変換を H, G と書き， H(ω) :=
1 √ 2 H(ω) H(ω + π) G(ω) G(ω + π) (44) とおく．補題 7 と補題 8 を満たすならば，H(ω) はユニタリ行列（H∗(ω) = H−1(ω)）である 17． 17（証明）H(ω) と H∗(ω) の積をとれば， H(ω)H∗(ω) = 1 2 H(ω) H(ω + π) G(ω) G(ω + π) H(ω) G(ω) H(ω + π) G(ω + π) = 1 2 H(ω)H(ω) + H(ω + π)H(ω + π) H(ω)G(ω) + H(ω + π)G(ω + π) G(ω)H(ω) + G(ω + π)H(ω + π) G(ω)G(ω) + G(ω + π)G(ω + π) = 1 2 2 0 0 2 （∵ 補題 7 と補題 8） = I 28 / 59

式(18)の証明補題9 補題 9 (FWT の行列表現) H(ω) がユニタリ行列であれば，次が成り立つ： Pm (ω)
Qm (ω) = 1 √ 2 H ω 2 Pm+1 ω 2 Pm+1 ω 2 + π (45) Pm+1(ω) Pm+1(ω + π) = √ 2H∗(ω) Pm (2ω) Qm (2ω) (46) ここで，Pm , Qm は pm , qm の離散時間フーリエ変換である． 29 / 59

式(18)の証明補題9 （証明）y[n] := k h[k − n]pm+1 [k] の離散時間フーリエ変換
Y (ω) は， Y (ω) = n k h[k − n]pm+1 [k] exp(−jnω) = k pm+1 [k] n h[k − n] exp(−jnω) = k pm+1 [k] exp(−jkω) l h[l] exp(−jlω) (l = k − n) = Pm+1 (ω)H(ω) pm [n] = y[2n] だから，間引いた数列の離散時間フーリエ変換の公式（式 (56)）を使うと， Pm (ω) = 1 2 Y ω 2 + Y ω 2 + π = 1 2 Pm+1 ω 2 H ω 2 + Pm+1 ω 2 + π H ω 2 + π (47) 同様にして， Qm (ω) = 1 2 Pm+1 ω 2 G ω 2 + Pm+1 ω 2 + π G ω 2 + π (48) 30 / 59

式(18)の証明補題9 式 (47)，式 (48) の結果を行列演算で表すと， Pm (ω) Qm (ω)
= 1 2 H ω 2 H ω 2 + π G ω 2 G ω 2 + π Pm+1 ω 2 Pm+1 ω 2 + π ⇐⇒ Pm (ω) Qm (ω) = 1 √ 2 H ω 2 Pm+1 ω 2 Pm+1 ω 2 + π また，両辺左から √ 2H∗ ω 2 を乗じると， √ 2H∗ ω 2 Pm (ω) Qm (ω) = Pm+1 ω 2 Pm+1 ω 2 + π が得られる．式 (46) は ω を 2ω に置き換えて得られる． 31 / 59

式(18)の証明次の目標補題 8 を満たすならば，任意の m, p, q ∈ Z
に対し， ψm,p , ψm,q = 2m ∞ −∞ ψ(2mt − p)ψ(2mt − q)dt = δpq （∵ ψ はシフトに関して正規直交）だから ψm,n は Wm の正規直交系となる．また， ψm,n で張られる Wm と φm,n で張られる Vm が直交するので Vm ⊥ Wm ． ⇒ 次に，Vm ⊕ Wm = Vm+1 を示す． 32 / 59

式(18)の証明補題10 補題 10 (直交直和分解可能な条件) {φ(t − n)|n ∈ Z}
は正規直交系かつ，補題 7 と補題 8 を満たすならば， Vm ⊕ Wm = Vm+1, Vm ⊥ Wm (m ∈ Z) （証明）Vm ⊕ Wm ⊂ Vm+1, Vm ⊥ Wm の前提において，x(t) ∈ Vm+1 (Vm ⊕ Wm ) なる要素（関数）を任意にとる．これが x(t) = 0 しかないことが示せれば， Vm ⊕ Wm ⊃ Vm+1，従って，Vm ⊕ Wm = Vm+1 が結論できる．いま Vm ⊥ Wm かつ x / ∈ Vm ⊕ Wm だから，Vm, Wm への射影成分はそれぞれ pm [n] = x, φm,n = 0, qm [n] = x, ψm,n = 0 となる．これらを両辺離散時間フーリエ変換すると Pm (ω) = Qm (ω) = 0 を得る．これを式 (46) に代入すると Pm+1 (ω) = 0 となり，これは時間領域で pm+1 [n] = x, φm+1,n = 0 となることを示している． {φm+1,n (t)|n ∈ Z} が Vm+1 の正規直交基底であることを踏まえると，x(t) = 0 であることが導かれる．よって，Vm ⊕ Wm = Vm+1 が示された． 33 / 59

式(18)の証明定理2 定理 2 (MRA を構成する G(ω) の条件) H(ω) は補題
7 を満たすとする．このとき，補題 8 を満たす G(ω) は次の形の関数であり，かつそのときに限られる： G(ω) = γ(ω)H(ω + π) (49) γ(ω) は周期 2π の周期関数で，次を満たす： γ(ω) + γ(ω + π) = 0, |γ(ω)|2 = 1 (50) 34 / 59

式(18)の証明定理2 （⇒ の証明）補題 8 の式 (41) より， H(ω)G(ω)+H(ω +
π)G(ω + π) = 0 ⇐⇒ G(ω) = −H(ω + π)G(ω + π)/H(ω) ⇐⇒ G(ω) = −H(ω + π)G(ω + π)/H(ω) = −γ(ω)H(ω + π) (51) ここで γ(ω) := −G(ω + π)/H(ω) と定めた．G, H は周期 2π だから γ(ω) も周期 2π で， γ(ω + π) = −G(ω + 2π)/H(ω + π) = −G(ω)/H(ω + π) = −γ(ω) （∵ 式 (51)）また，式 (51) を補題 8 の式 (42) に代入すると， |G(ω)|2 + |G(ω + π)|2 = 2 ⇐⇒ |γ(ω)|2 |H(ω + π)|2 + |γ(ω + π)|2 |H(ω + 2π)|2 = 2 ⇐⇒ |γ(ω)|2 |H(ω + π)|2 + |γ(ω)|2 |H(ω)|2 = 2 ⇐⇒ |γ(ω)|2(|H(ω + π)|2 + |H(ω)|2) = 2 ⇐⇒ |γ(ω)|2 = 1 （∵ 補題 7 の式 (38)） 35 / 59

式(18)の証明定理2 （⇐ の証明）γ(ω) + γ(ω + π) = 0,
|γ(ω)|2 = 1 を満たす周期 2π の関数 γ(ω) をおく．いま G(ω) := γ(ω)H(ω + π) により G(ω) を定めると， |G(ω)|2 + |G(ω + π)|2 = |γ(ω)|2 |H(ω + π)|2 + |γ(ω + π)|2 |H(ω + 2π)|2 = |γ(ω)|2 |H(ω + π)|2 + | − γ(ω)|2 |H(ω)|2 = |H(ω + π)|2 + |H(ω)|2 = 2 よって式 (41) が成立している．また， H(ω)G(ω) + H(ω + π)G(ω + π) = H(ω)γ(ω)H(ω + π) + H(ω + π)γ(ω + π)H(ω + 2π) = H(ω)H(ω + π) γ(ω) + γ(ω + π) = H(ω)H(ω + π){γ(ω) + γ(ω + π)} = 0 よって式 (42) が成立している．以上により，補題 8 が成立している． 36 / 59

式(18)の証明系 1 (式 (18) の証明) 式 (18) は定理 2
において γ(ω) = − exp(−jω) とおいたときに相当する．（証明）G(ω) = − exp(−jω)H(ω + π) の両辺を逆離散時間フーリエ変換すると， g[n] = − 1 2π 2π 0 H(ω + π) exp(−jω) exp(jωn)dω = − 1 2π 2π 0 H(s) exp[jω(s − π)(n − 1)]ds (s = ω + π) = − 1 2π exp[−jπ(n − 1)] 2π 0 H(s) exp[jωs(n − 1)]ds = (−1)n 1 2π 2π 0 H(s) exp[jωs(1 − n)]ds = (−1)nh[1 − n] 37 / 59

C. 不確定性原理 38 / 59

不確定性原理信号 f は lim t→±∞ |t|f(t) = 0 と
∞ −∞ |f(t)|2dt = 1 を満たすとする．今，|f(t)|2 を確率密度と見よう．パーセバルの等式 18 より， 1 2π ∞ −∞ |F(ω)|2dω = 1 だから， 1 2π |F(ω)|2 も確率密度と見れる． 18フーリエ変換によってパワーは変わらない．冗長なので末尾で示す． 39 / 59

不確定性原理 t, ω をそれぞれ確率変数として，真値 tm , ωm に対する標準偏差 ∆t,
∆ω を以下で定義する： (∆t)2 := ∞ −∞ (t − tm )2|f(t)|2dt (52) (∆ω)2 := 1 2π ∞ −∞ (ω − ωm )2|F(ω)|2dω (53) ∆t, ∆ω は観測の揺らぎ幅とも解釈できる． 40 / 59

不確定性原理定理 3 (不確定性原理) ∆t, ∆ω に対して以下の不等式が成立する： ∆t∆ω ≥ 1
2 (54) 観察する時間の幅 ∆t と周波数の幅 ∆ω を同時に小さくすることはできない． 41 / 59

不確定性原理現象 1 短時間フーリエ変換でより広い時間情報を得ようとして窓幅を大きく取ると (∆t：小) 周波数分解能が小さく（狭く）(∆ω：大) なる現象
2 ウェーブレット変換でスケール a を小さく取ると細かい振動を捉えられる（∆ω：小）が，時間幅は狭く（∆t：大）なる現象 3 フーリエ変換は周波数分解能が最大 (∆ω → 0) だが，時間分解能は最悪 (∆t → ∞) 42 / 59

不確定性原理時間-周波数領域で現象を描くと，この図に現れる矩形領域をハイゼンベルクの箱と言う． 43 / 59

不確定性原理の証明また， f (t) = d dt 1 2π ∞
−∞ F(ω) exp(jωt)dt = 1 2π ∞ −∞ F(ω) d dt exp(jωt)dt （順序交換可能とする） = jω 2π ∞ −∞ F(ω) exp(jωt)dt = jωf(t) より，以下の関係式が得られる． F f (t) = jωF [f(t)] = jωF(ω) =⇒ F f (t) 2 = ω2 |F(ω)|2 これより， 1 2π ∞ −∞ ω2 |F(ω)|2dω = 1 2π ∞ −∞ F f (t) 2 dω = ∞ −∞ |f (t)|2dt （∵ パーセバルの等式） 45 / 59

不確定性原理の証明不等式左辺（= (∆t)2(∆ω)2 ）を変形していくと， ∞ −∞ t2 |f(t)|2dt 1 2π
∞ −∞ ω2 |F(ω)|2du = ∞ −∞ t2 |f(t)|2dt ∞ −∞ |f (t)|2du ≥ ∞ −∞ tf(t)f (t) 2 dt （∵ シュワルツの不等式 ||v||2||w||2 ≥ || v, w ||2） ≥ 1 4 ∞ −∞ tf(t)f (t) + tf(t)f (t) 2 dt (∵ |ab| ≥ 1 2 |ab + ab|) ≥ 1 4 ∞ −∞ tf(t)f (t) + tf(t)f (t) dt 2 = 1 4 ∞ −∞ t d dt f(t)f(t) dt 2 （∵ 積の微分公式） = 1 4 t|f(t)|2 ∞ −∞ − ∞ −∞ |f(t)|2dt 2 （∵ 部分積分と f(t)f(t) = |f(t)|2 ） = 1 4 {0 − 1}2 （∵ f(t) に関する仮定） = 1 4 よって不確定性原理は示された． 46 / 59

D. フーリエ変換における諸性質 47 / 59

パーセバルの等式（フーリエ級数）周期 2π の関数 f(x), g(x) をそれぞれ以下のようにフーリエ級数展開できるとする： f(x) = n
cn exp(jωx), cn = 1 2π π −π f(x) exp(−jnω)dx g(x) = n dn exp(jωx), dn = 1 2π π −π g(x) exp(−jnω)dx このとき， π −π f(x)g(x)dx = π −π n cn exp(jωx) g(x)dx = n cn π −π g(x) exp(jnx)dx = n cn π −π g(x) exp(−jnx)dx = 2π n cndn が成立する．とくに f = g とすれば， π −π |f(x)|2dx = 2π n |cn|2 48 / 59

パーセバルの等式（フーリエ変換）関数 f(x), g(x) のフーリエ変換対を以下のように書く： F(ω) = ∞ −∞ f(x)
exp(−jωx)dx, f(x) = 1 2π ∞ −∞ F(ω) exp(jωx)dω G(ω) = ∞ −∞ g(x) exp(−jωx)dx, g(x) = 1 2π ∞ −∞ G(ω) exp(jωx)dω このとき， ∞ −∞ f(x)g(x)dx = ∞ −∞ 1 2π ∞ −∞ F(ω) exp(jωx)dω g(x)dx = 1 2π ∞ −∞ F(ω) ∞ −∞ g(x) exp(jωx)dx dω = 1 2π ∞ −∞ F(ω) ∞ −∞ g(x) exp(−jωx)dx dω = 1 2π ∞ −∞ F(ω)G(ω)dω が成立する．とくに f = g とすれば， ∞ −∞ |f(x)|2dx = 1 2π ∞ −∞ |F(ω)|2dω 49 / 59

パーセバルの等式（高速ウェーブレット変換）補題 11 (パーセバルの等式) 高速ウェーブレット変換の列 {pm+1 [n]}，{pm [n]}， {qm [n]}
に関して次が成立する： ||pm+1||2 2 = ||pm ||2 2 + ||qm ||2 2 (55) 入力信号のエネルギー ||pm+1||2 2 は ||pm ||2 2 と ||qm ||2 2 に分配される． 50 / 59

パーセバルの等式（高速ウェーブレット変換）（証明）式 (45) の両辺のノルムをとると， Pm (ω) Qm (ω) 2 2
= 1 2 H ω 2 Pm+1 ω 2 Pm+1 ω 2 + π 2 2 = 1 2 Pm+1 ω 2 Pm+1 ω 2 + π 2 2 （∵ H(ω) はユニタリ行列） ⇐⇒ |Pm (ω)|2 + |Qm (ω)|2 = 1 2 Pm+1 ω 2 2 + Pm+1 ω 2 + π 2 離散時間フーリエ変換のパーセバルの等式により， ||pm||2 2 + ||qm||2 2 = 1 2π π −π |Pm (ω)|2 + |Qm (ω)|2 dω = 1 4π π −π Pm+1 ω 2 2 + Pm+1 ω 2 + π 2 dω = 1 2π π 2 − π 2 |Pm+1 (s)|2 + |Pm+1 (s + π)|2 ds (ω = 2s) = 1 2π 3 2 π − π 2 |Pm+1 (s)|2ds = 1 2π π −π |Pm+1 (s)|2ds = ||pm+1||2 2 51 / 59

間引いた数列の離散時間フーリエ変換補題 12 (間引いた数列の離散時間フーリエ変換) x[n] の離散時間フーリエ変換を X(ω) と書く． x[n] のサンプリング間隔を
D 倍に間引いた数列の離散時間フーリエ変換は， 1 D D−1 k=0 X ω − 2πk D (56) 52 / 59

間引いた数列の離散時間フーリエ変換（証明）間引き操作は，x[n] に周期 D のインパルス関数列 δD [n] := 1 n
が D の倍数 0 otherwise = 1 D D−1 k=0 exp j 2πnk D （∵ 証明は次のページ）を乗じてから，D 個おきの出力を得る操作 (y[n] := x[n]δD [n] として y[Dn] が結果) と解釈できる．y[Dn] を離散時間フーリエ変換すると， F [y[Dn]] = n y[Dn] exp(−jnω) = n y[n] exp −j nω D = n x[n]δD [n] exp −j nω D = 1 D D−1 k=0 n x[n] exp j 2πnk D − j nω D （∵ 式 (57)） = 1 D D−1 k=0 n x[n] exp −jn ω − 2πk D = 1 D D−1 k=0 X ω − 2πk D 53 / 59

間引いた数列の離散時間フーリエ変換 δD [n] = 1 D D−1 k=0 exp j
2πnk D (57) （証明）n が D の倍数のときは n = mD (m ∈ Z) と書けるから， 1 D D−1 k=0 exp j 2πnk D = 1 D D−1 k=0 exp(j2πmk) = 1 D D−1 k=0 1 = 1 n が D の倍数ではないときは Wk := exp j 2πnk D とおくと， (Wk )D = exp(j2πnk) = 1 より，Wk (k = 0, 1, ..., D − 1) は 1 の D 乗根，すなわち，方程式 xD − 1 = 0 の解となっている．方程式左辺を因数分解すると， xD − 1 = (x − W0 )(x − W1 )...(x − WD−1 ) = xD − (W0 + W1 + ... + WD−1 )xD−1 + ... + (−1)DW0W1...WD−1 両辺の係数を比較する．xD−1 の係数は 0 だから， W0 + W1 + ... + WD−1 = D−1 k=0 Wk = D−1 k=0 exp j 2πnk D = 0 54 / 59

E. 参考文献 55 / 59

参考資料の概観 web で入手できる資料： [4, 5, 6] 大学講義資料．取っ掛かりに． [7, 8] 分かりやすい．C++ソース解説も充実．
[9] Lifting の簡易な説明．ソースもある． [10, 11, 12] Lifting の資料．とくに [11] は簡明．書籍： [13] 分かりやすい．ただしハールウェーブレットのみ． [14] 分かりやすい．C 言語ソースあり． [15] 分かりやすいが絶版...C++ソースあり． [16] 歴史を概観するのに適する．意外に理論も深い． [2] 一番参考にした．信号処理から見た解説も秀逸． [1] ドベシィ様執筆．基礎理論を網羅．読みきれてない． [17] 関数解析必須．（挫折） 56 / 59

参考文献 I [1] 佐々木文夫 Daubechies Ingrid 山田道夫. ウェーブレット
10 講. 丸善出版, 2012. [2] 前田肇 et al. ウェーブレット変換とその応用. 朝倉書店, 2001. [3] Albert Boggess and Francis J Narcowich. A ﬁrst course in wavelets with Fourier analysis. John Wiley & Sons, 2015. [4] 伊藤彰則. Wavelet 変換. url: http://www.spcom.ecei.tohoku.ac.jp/˜aito/wavelet/slide.pdf (visited on 10/08/2020). [5] 羽石秀昭. Wavelet1. url: http://www.cfme.chiba- u.jp/˜haneishi/class/iyogazokougaku/Wavelet1.pdf (visited on 10/08/2020). [6] 羽石秀昭. Wavelet2. url: http://www.cfme.chiba- u.jp/˜haneishi/class/iyogazokougaku/Wavelet2.pdf (visited on 10/08/2020). [7] Fussy. 圧縮アルゴリズム (8) ウェーブレット変換 -1-. url: http://fussy.web.fc2.com/algo/compress8_wavelet.htm (visited on 10/09/2020). 57 / 59

参考文献 II [8] Fussy. 圧縮アルゴリズム (9) ウェーブレット変換 -2-. url: http://fussy.web.fc2.com/algo/compress9_wavelet2.htm
(visited on 10/09/2020). [9] Ian Kaplan. Basic Lifting Scheme Wavelets. url: http: //bearcave.com/misl/misl_tech/wavelets/lifting/basiclift.html (visited on 10/09/2020). [10] 藤ノ木健介. “リフティングウェーブレットと信号解析”. In: 2014. [11] Geert Uytterhoeven, Dirk Roose, and Adhemar Bultheel. “Wavelet transforms using the lifting scheme”. In: (1997). [12] 藤ノ木健介. “リフティングスキームによるウェーブレットの構成法”. In: 日本応用数理学会論文誌 28.2 (2018), pp. 72–133. [13] 金谷健一. これなら分かる応用数学教室–最小二乗法からウェーブレットまで–. 共立出版. [14] 山本鎮男中野宏毅吉田靖夫. ウェーブレットによる信号処理と画像処理. 共立出版. 58 / 59

参考文献 III [15] 川畑洋昭戸田浩章忠. 新ウェーブレット実践講座
: 入門と応用 : 信号処理の基礎から最新理論まで. C magazine. ソフトバンククリエイティブ, 2005. [16] 西野操 B.B. ハバード山田道夫. ウェーブレット入門―数学的道具の物語. 朝倉書店, 2003. [17] 新井勉 C.K.Chui 桜井明. ウェーブレット入門. 東京電機大学出版局. 59 / 59

ウェーブレットおきもち講座

ウェーブレットおきもち講座

Other Decks in Science

Featured

Transcript