Wavelet 変換 伊藤 彰則 aito@fw.ipsj.or.jp 1
参考書 (1) 中村, 山本, 吉田 : ウェーブレットによる信号処理と画像処理, 共立出版 応用の紹介とプログラムリストが中心, 理論的背景はほとんどなし 意味不明の比喩を多用 各時代 各国別に美女を探すのが窓フーリエ変換である 応用テーマ : 不連続信号検出, 相関の検出, ノイズ除去, 画像データ圧縮, 劣化画像復元 芦野, 山本 : ウェーブレット解析 - 誕生 発展 応用, 共立出版 理論編と応用編に分かれている. 比較的わかりやすい Daubechies と Meyer のウェーブレットが中心, その他の紹介はあまりない 応用テーマ : 不連続信号検出, 相関の検出, データ操作 2
参考書 (2) 新井 : ウェーブレット解析の基礎理論 森北出版 前半が理論, 後半が応用だが理論的背景は薄い 理論を飛ばして読む目的には良いが, 理論をこれだけで理解するのは難しい 対象ウェーブレットを幅広く扱っている 応用編はプログラムリストつき 応用テーマ : エッジ抽出, データ圧縮, 積分方程式の数値解法 C.K.Chui, 桜井 新井訳 : ウェーブレット入門, 東京電機大学出版局 理論のみ, 応用事例なし 理論を概観する章がある. 証明抜きで全体像をつかむには良い 扱うのはほとんどスプラインウェーブレット. 3
参考書 (3) 前田, 佐野, 貴家, 原 : ウェーブレット変換とその応用, 朝倉書店 前半が理論, 後半が応用 厳密な理論展開は省略, 最小限の証明 積分 Wavelet 変換についてちょっと詳しい 扱う Wavelet は Haar, MexicanHat,Poisson など 応用テーマ : レーダー, システム同定, 雑音抑圧 B.B.Hubbard, 山田 西野訳 : ウェーブレット入門, 朝倉書店 インタビューにもとづくウェーブレット発見の物語 誰がどういう経緯で何を発見したのかが書いてある 著者は数学者ではない 数学的記述を全く見ないで読むことも可能 4
関数展開 Wavelet の仲間たち 窓関数 Fourier 変換 関数展開 窓 Fourier 変換 離散入力 z 変換 スケーリング関数 窓関数の一般化時間 - 周波数依存性 多重解像度解析 連続 Wavelet 変換 Gabor Wavelet 離散 Wavelet 変換 Fourier 展開 ツースケール関係 離散入力 共役ミラーフィルタ (QMF) 高速 Wavelet 変換 (Mallat アルゴリズム ) Daubechies Wavelet 5
フーリエ変換 無限に長い波形からスペクトルを計算する フーリエ変換 f = f x e i x dx フーリエ逆変換 f x = f x e i x d 6
時系列の解析 時間とともにスペクトルが変化する場合? 上記の例ではだんだん周波数が上がっている 全体を解析したのでは だんだん上がる という分析は不可能 7
窓フーリエ変換 信号に窓関数をかけてフーリエ変換 窓関数をずらすことで時間変化を調べる F t, = = f x w x t e i x dx f x,t x dx 8
時間ー周波数分析 窓関数をずらすことでスペクトルの時間変化を調べる 9
スペクトログラム スペクトルの強度を濃淡で表示 ω t 10
窓関数 信号を時間的に切り取る関数 時間 t が大きいところ 小さいところでは 0 になるのが望ましい ( コンパクト台 ) t 2 2 w x =e w x ={ 1 1 cos t 2 1 t 1 0 otherwise w x ={ 1 1 t 1 0 otherwise w x ={ 0.54 0.46cos t 1 t 1 0 otherwise 11
窓関数 窓関数の条件 x w x dx 窓関数の中心と幅 w 2 = w x 2 dx x = 1 w 2 x w x 2 dx 1 w 2 x x 2lline w x w = 2dx 12
窓フーリエ変換の不確定性 窓が長いほうがスペクトルの解像度が上がる 接近した 2 つの周波数成分の区別が可能に 窓が短いほうが時間の解像度が上がる 周波数の細かい時間変化を捉えることができる 窓関数の時間幅 周波数幅 < 定数 ( デモプログラム ) 13
窓関数で信号を解析する F t, = = f x w x t e i x dx f x,t x dx 窓関数 ψ を使った解析をしているとみなせる ガウス窓 + 指数関数 ガボール (Gabor) 変換 14
連続 Wavelet 変換 窓関数 (Wavelet) を使った変換 W f b, a = 1 a f x x b a dx a: 周波数の逆数に相当 ( ダイレーション ) b: 時間に相当 ( シフト ) p: アナライジングウェーブレット 15
Wavelet の条件 直流信号を解析してもバイアスがないこと x dx=0 逆変換が存在すること x 2 x dx= 1 2 C 16
連続 Wavelet 逆変換 次の式により逆変換が可能 f x = 2 C 0 [ W f b,a 1 a x b a db ] da a 2 ただし f(x) は次の条件を満たす f x L R 2 すなわち f x 2 dx 17
アナライジングウェーブレットの例 Haar wavelet 0 x 0 1 0 x 0.5 x ={ 1 0.5 x 1 0 1 x Mexican hat wavelet x = 1 x 2 e x 2 2 18
演習 cos(kx) を Haar で Wavelet 変換してみよう こんな風になるはず 下の図では a は対数スケールであることに注意 a は周波数の逆数 k が大なら a は小 19
連続 Wavelet 変換の特徴 低い周波数では窓幅が広く 高い周波数では窓幅が狭い 低い周波数に対しては高い周波数分解能 高い周波数に対しては高い時間分解能 ( 対数周波数領域での分解能が一定 ) 20
離散 Wavelet 変換 連続 Wavelet 変換 :1 変数 2 変数 W f b,a = 1 a 離散点の係数の総和で元の関数を表現 離散 Wavelet 変換 ( Wavelet 展開 ) f x = j= W f k k= 2, 1 j 2 j f x x b a dx 2 j/ 2 2 j x k = j = W f k k = 2, 1 j j 2 jk x ここは定数 21
前式の定数の部分 Wavelet 展開係数 c jk = W f k 2 j, 1 2 j f x = j= k= c jk jk x とおくと Wavelet 展開係数 係数は連続 Wavelet 変換の離散点での値 22
Wavelet 展開係数について注意 a は周波数の逆数に比例 a と b でのサンプル点 低周波では時間方向に粗く周波数方向に細かい 周波数 f=1/a と b でのサンプル点 高周波では時間方向に細かく周波数方向に粗い 23
Wavelet 基底関数 (1) Wavelet 展開は 2 次元の関数展開 f x = j= 参考 : Taylor-McLaurin 展開 Fourier 展開 jk x f x = k = f x = k = k= a k x k a k e ikx c jk jk x は Wavelet 展開の基底関数 24
基底関数 Wavelet 基底関数 (2) 基底関数の無限和で( ある条件下の ) 任意の関数が表せる これだけでは係数の一意性は保証されない 正規直交基底 jk x j ' k ' x dx= jj ' kk ' 正規直交基底を使うと関数を一意に展開できる Haar 関数に基づく ψ jk は最も簡単な正規直交基底 2 j / 2 H 2 j x k 2 j '/ 2 H 2 j ' x k dx= jj ' kk ' 演習 : 上のことを証明せよ 25
多重解像度解析 入力信号を周波数帯域ごとに分ける 26
多重解像度解析 各帯域幅で表される関数の集合間の関係 V 2 V 1 V 0 V 1 V 2 V 0 の正規直交基底 : スケーリング関数 f V 0 f x = k = c k x k 最も簡単なスケーリング関数 x ={ H 1 0 x 1 0 otherwise 1 O 1 演習 : 上記のスケーリング関数が正規直交基底であることを示せ 27
スケーリング関数による関数の分解 前述の H x で表現される関数は? 任意の整数 kについて k x k 1 ならば f x = f k 整数点でサンプリングされた関数 28
注意 しばらくは このスケーリング関数 (Haar のスケーリング関数 ) を使って話を進める 簡単だから ほかにもスケーリング関数はある Wavelet との関係は? 多重解像度解析 スケーリング関数 ツースケール関係 Wavelet 関数 離散入力 高速 Wavelet 変換 = 特殊な条件での離散 Wavelet 変換 29
スケーリング関数による関数の分解 この関数をスケーリング関数により分解 c 6 H x 6 c 5 H x 5 c 4 H x 4 c 3 H x 3 各時間の値だけを持つ関数の重ねあわせで関数を表現する もっと細かい関数を表現したい場合は? 30
いろいろな細かさのスケーリング関数 1 O 1 1 O H 0.5x 1 1 O H x 1 H 2x f x V 1 f x V 0 f x V 1 31
スケーリング関数の合成 細かいスケーリング関数から粗いスケーリング関数を作ることができる x = k= p k 2x k Haar のスケーリング関数の場合 H x = H 2x H 2x 1 1 1 1 O 1 = O + 1 O 1 p 0 =1, p 1 =1 p k =0 k 0, k 1 32
スケーリング関数の分解 (1) とを合わせて 2x を表現できる x x x = x V 1 V 0 qk 2x k V 0 V 1 は作れるか? 基底関数 2x k 基底関数 x k W 0 基底関数 x k 33
スケーリング関数の分解 (2) x, x, 2x の関係式 x = x = p k 2x k q k 2x k ツースケール関係 (two-scale relation) ただし q k = 1 k p 1 k これでいい事を証明してみよう Haar のスケーリング関数の場合 x = 2x 2x 1 x = 2x 2x 1 x は Haar ウェーブレット関数 34
スケーリング関数の分解 (3) スケーリング関数 2x をで表す 2x l = 1 2{ k= g k = p k h k =q k Haar のスケーリング関数の場合 x, x g 2k l x k h 2k l x k } H 2x = 1 2 { H x H x } H 2x 1 = 1 2 { H x H x } 35
注意点 W 0 の正規直交基底だけでは Wavelet とは言えない Wavelet になるためには他の条件も必要, 例えば x 2 x dx= 1 2 C スケーリング関数の分解 合成は離散 wavelet 変換と深い関係がある! これからスケーリング関数を使った信号の分解 合成について説明する 36
信号の分解と再構成 (1) 信号 ( 関数 ) f N x V N f N x = k= これを次のように分解する f N 1 x = k= g N 1 x = k= c k N 2 N x k c k N 1 2 N 1 x k d k N 1 2 N 1 x k f N x = f N 1 x g N 1 x 37
信号の分解と再構成 (2) c N,c N 1, d N 1 の関係 分解 c k N 1 = 1 2 l c l N g 2k l d k N 1 = 1 2 l c l N h 2k l 再合成 c k N = l {c l N 1 p k 2l d l N 1 q k 2l } 演習 : 上記の関係を導出せよ. 38
信号の分解と再構成 (3) Haar のスケーリング関数の場合 分解 c k N 1 = 1 2 c 2k N N c 2k 1 d k N 1 = 1 2 c 2k N N c 2k 1 再構成 N c 1 N 2k =c 1 N 1 k d k N c 1 N 2k 1 =c 1 N 1 k d k f N x = k c k N H 2 N x k f N 1 x = k 連続する 2 点の平均 c k N 1 H 2 N 1 x k 連続する 2 点の差の半分 f N 1 x = k d k N 1 H 2 N 1 x k 39
離散 Wavelet 変換との関係 (1) 信号の分解を繰り返すと... f N x = g N 1 x f N 1 x =g N 1 x g N 2 x f N 2 x =g N 1 x g N 2 x g N 3 x f N 3 x... =g N 1 x g N 2 x g N M x f N M x f N がバイアス ( 直流成分 ) を含まなければ N 1 f N x = j= N 1 g j x = j= k d k j 2 j x k 40
離散 Wavelet 変換との関係 (2) 信号の分解 N 1 k f N x = j= 離散 Wavelet 変換 N 1 k f N x = j= d k j 2 j x k W f N k 2 j, 1 2 j 2 j /2 2 j x k 正規直交 Wavelet による展開係数は一意に決まるので d k j = W f N k 2 j, 1 2 j 2 j / 2 つまり f N に対してはスケーリング関数による分解 = 離散 Wavelet 変換 41
高速 Wavelet 変換 f N についての離散 Wavelet 変換...f N って? Haar のスケーリング関数の場合 : 整数の 2 -N 倍の区間で定常な関数 離散点でサンプリングされた信号の場合 : その点を f N とみなすことで分解が可能 高速 Wavelet 変換 (Mallat アルゴリズム ) ( 本当はそれではマズイ場合もあることには注意しておこう ) 42
周波数領域で見た Wavelet ツースケール関係 : 時間領域 周波数領域で見ると? 高次の φ N : 広い周波数帯域 狭い台 低次の φ N : 狭い周波数帯域 広い台 ψ N : φ N を補完する周波数特性 φ N φ N+1 ψ N 43
準備 ツースケール関係 x = p k 2x k x = これらをフーリエ変換すると q k 2x k = k = i k 1 2 p e 2 k 2 =m 0 2 2 = k= i k 1 2 q e 2 k 2 =m 1 2 2 m 0 = k= 1 2 p k e i k m 1 = k= 1 2 q i k k e 44
Wavelet の周波数特性 (1) =m 0 2 2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 H 2 2 = /2 1 e i i /2-15 -10-5 0 5 10 15 2 45
Wavelet の周波数特性 (2) =m 0 2 2 1 m 0 2 2 = 1 e i /2 2 2 0.8 0.6 0.4 0.2 0-15 -10-5 0 5 10 15 46
Wavelet の周波数特性 (3) =m 0 2 2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 H 2 H 2 2-15 -10-5 0 5 10 15 47
Wavelet の周波数特性 (4) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 =m 1 2 2-15 -10-5 0 5 10 15 m 1 2 = 1 e i /2 2 48
Wavelet の周波数特性 (5) =m 1 2 2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 H 2 H 2 2-15 -10-5 0 5 10 15 49
フィルタとしての Wavelet m 0 はローパスフィルタ, m 1 はハイパスフィルタ 1 0.8 0.6 m 0 m 1 共役ミラーフィルタ (QMF) 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 50
Haar 以外の Wavelet Wavelet 解析では Haar 以外にもさまざまな Wavelet が用いられる 連続 Wavelet 変換では Gabor, Maxican hat,... 離散 Wavelet 変換では Spline, Daubechies, Meyer,... 工学的に重要な Daubechies の Wavelet について解説する 51
Daubechies の Wavelet (2 次 ) ツースケール関係 x = 1 3 8 x = 1 3 8 2x 3 3 8 2x 3 3 8 2x 1 3 3 8 2x 1 3 3 8 φ と ϕ は既存の関数では表せない 2x 2 1 3 8 2x 2 1 3 8 2x 3 2x 3 x x 52
Daubechies の Wavelet の周波数特性 周波数領域での表現 m 0 = 1 3 8 3 3 8 e i 3 3 8 e 2i 1 3 8 e 3i 1 0.8 0.6 Haar Daubechies(2) Daubechies は Haar よりも急峻なフィルタ特性 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 53
Wavelet の応用 (1) いくつかの Wavelet の応用について見てみよう 波形中の不連続点の検出 機械の破壊による振動の検知など ノイズ除去 ノイズについての知識を使わずにノイズを除去する 信号圧縮にも有効 画像処理 画像を細かい部分と大まかな部分に分ける ピラミッドアルゴリズム 54
Wavelet の応用 : 不連続点の検出 この波形で不連続なのはどこでしょう? 55
不連続点の検出 現波形に対して Haar と Daubechies(2) Haar Daubechies 56
不連続点の検出 正弦波重畳波形に対して Haar と Daubechies Haar Daubechies Daubechies の方が周波数分離能力が高いことがわかる 57
Wavelet の応用 : ノイズ除去 処理手順 Wavelet 展開係数を求める 閾値を決め, それより小さい係数を 0 にする 元の信号を再構成 先ほどの波形にノイズを載せてみました 58
ノイズ除去結果 Haar Daubechies(2) Daubechies(3) 59
Wavelet の応用 : 画像処理 2 次元 Wavelet の例 (Haar) X,Y 方向にそれぞれ Wavelet 変換を適用 60
画像処理 分析してどうするのか? 高周波部分の量子化ビット数を落とす 圧縮 高周波部分を0に 平滑化 ノイズ除去 再帰的な分解 61