Wavelet - mynote

メニュー

recent(15)

2024-02-01

GIMPで縦書き

2024-01-25

TkSQLite

2023-04-18

勉強中のメモとかポインタです。

Haar
Daubechies
カスケードアルゴリズム
多重解像度解析
リンク
ノイズ除去
Tclによるカスケード分割の書き捨てプログラム
TclによるD4 Wavelet変換逆変換の書き捨てプログラム
コメントをどーぞ

Haar†

足して2で割る。または√2で割る。Waveletの紹介でよく紹介される方法は足して２で割っているが、正規化したデータを使う解きは√2の方を用いるようだ？

↑

Daubechies†

DaubechiesのWaveletはHaarを拡張したもので、D2の時はHaarの基底と同じ。D4、D6、D8…と長めのフィルターがある。プログラムするときはすでに係数を計算してあるのを用意して使う。D20くらいまではWebで検索すると見つけられるようです。分解と再構築で同じ係数を使う。
D2の時はともかく端のところで折り返さないとならないので、そのへんが微妙に不安な気がする…。
ドベシーおばさんの若かりし頃→写真
小学校で担任だった先生に激似。

↑

カスケードアルゴリズム†

スケーリング関数を描画する時に使う方法で、N=2のとき

φ(t) = h(3)φ(2t) + h(2)φ(2t-1) + h(1)φ(2t-2) + h(0)φ(2t-3)

を満たすよう最初にφ(0)～φ(3)の4点の値を求め、順次分割して計算していく方法。分割の処理はそんなに難しくないけど分割前の一番最初のtが整数の時のφ(t)を求めるのが難しいです。N=3以上の時の解き方がまだ全然わかりません。N=3の時は6次の代数方程式を解く（解けない）ことになるのでは・・・？N=3以上のプロットとかをWebで見かけるけど近似計算しているのか何か別の計算方法があるのか謎。解きかたがわかりました。

D4の場合
```
+-                -+
|  h0,  0,  0,  0  |
|                  |
|  h2, h1, h0,  0  |
|                  |
|   0, h3, h2, h1  |
|                  |
|   0,  0,  0, h3  |
+-                -+
```
このような行列の固有ベクトルを求めればφ(t)のtが整数の場合の値を求めることができる。固有ベクトルの計算は~~全然わからないので~~めんどくさいのでMuPADを使って事前に計算して、既に書いてあったTclのプログラムに食わせて分割処理。分割はtを2倍した一個上の解像度の値を元に計算できるので、順次解像度を上げながらループさせればＯＫ（？）。
マザーウェーブレットは計算したスケーリング関数の値を使って
```
ψ(t) = -h(0)φ(2t-1) + h(1)φ(2t) - h(2)φ(2t+1) + h(3)φ(2t+2)  
```
で簡単に求めることができる。
D8(N=8)をMuPADで固有ベクトルを計算して、Tclで補間処理して、gnuplotで描画してみた例

↑

多重解像度解析†

Mallat変換とかピラミッドアルゴリズムとか呼ばれる方法は分離された差分には手をつけないで、縮小画像に再びウェーブレット変換していくという方法で計算量がO(NlogN)になる。高速ウェーブレット変換(FWT)などとも呼ばれているような気がしないでもないが、FWTはO(N)で計算できるとか書いているページもあったりしてなんだか用語が曖昧でよくわかりません。
画像へWavelet変換を適用する場合、縦横と折りたたむように変換する。次の画像は自作の変換プログラムでHaarのWavelet変換をレベル-2まで適用してみた例。

レベル-0（何もしない）
レベル-1
レベル-2

↑

リンク†

↑

ノイズ除去†

自作ライブラリによるノイズ除去プログラム->AntiNoise

↑

Tclによるカスケード分割の書き捨てプログラム†

こんなんでいいのかね・・・。まあいいや・・・。gnuplotで描画できるデータを出力します。あとタップ係数と、事前に計算した固有ベクトルのファイルが別に必要です。まあソースを読んでください。

if {$argc <2} {
	error "format : tclsh plot.tcl D4 5 m|s"
}
set name [lindex $argv 0]
set level [lindex $argv 1]
set type [lindex $argv 2]
if {$type != {m}} {
	set type s
}

set f [list]
if {[catch {open $name\_f r} fp]} {
	error "can not find $name\_f"
}
foreach d [read $fp] {lappend f [expr $d * sqrt(2)]}
close $fp

set ev [list]
set fname $name\_ev
if {[catch {open $fname r} fp]} {
	error "can not find $fname"
}
foreach d [read $fp] {lappend ev [expr $d / (sqrt(2) / 2.0)]}
close $fp


if {[llength $f] != [llength $ev]} {
	error "bad data desu..."
}

set n [llength $f]
set len [expr ($n-1) * (1<<($level-1)) +1]

set data [list]
for {set i 0} {$i<$len} {incr i} {
	lappend data 0
}

for {set L 0} {$L<$level} {incr L} {
	if {$L==0} {
		for {set i 0} {$i<$n} {incr i} {
			set didx [expr $i * ($len-1) / ($n-1)]
			lset data $didx [expr [lindex $ev $i]]
		}
	} else {
		set offset [expr int(pow(2, $level-$L))]
		for {set i [expr $offset/2]} {$i<[llength $data]} {incr i $offset} {
			set tmp 0
			for {set t 0} {$t<$n} {incr t} {
				set didx [expr $i * 2 - $t*(1<<($level-1))]
				if {$didx<0} {continue}
				if {$didx>=[llength $data]} {continue}
				set tmp [expr $tmp + [lindex $f $t] * [lindex $data $didx]]
			}
			lset data $i $tmp
		}
	}
}

# ---- scaling ----
if {$type != {m}} {
	set i 0
	set ic [expr 2.0/(1<<$level)]
	foreach d $data {
		puts "$i	$d"
		set i [expr $i + $ic]
	}
	exit
}

# ---- mother ----
set i [expr -1 * ([llength $f]/2-1)]
set ic [expr 2.0/(1<<$level)]
for {set t [expr $i*(1<<($level-1))]} {$t<[llength $data]} {incr t} {
	set d 0.0
	set c [expr 2*$t - (1<<($level-1))]
	set ii -1
	foreach {g} $f {
		if {$c>=0 && $c<[llength $data]} {
			set d [expr $d + $ii*$g * [lindex $data $c]]
		}
		set ii [expr $ii * -1]
		incr c [expr 1<<($level-1)]
	}

	puts "$i	$d"
	set i [expr $i + $ic]
}

↑

TclによるD4 Wavelet変換逆変換の書き捨てプログラム†

練習を兼ねてTclで書いてみた検証用テストプログラム。D4をレベル-1だけ分解と再構築することができます。変換するリストは４の倍数でなければなりません。

set h0 0.4829629131445341
set h1 0.8365163037378079
set h2 0.2241438680420133
set h3 -0.1294095225512603
set g0 [expr $h3]
set g1 [expr -$h2]
set g2 [expr $h1]
set g3 [expr -$h0]


proc d4_h {a b c d} {
	global h0 h1 h2 h3
	return [expr $a*$h0 + $b*$h1 + $c*$h2 + $d*$h3]
}

proc d4_g {a b c d} {
	global g0 g1 g2 g3
	return [expr $a*$g0 + $b*$g1 + $c*$g2 + $d*$g3]
}

proc d4_ih {a b c d} {
	global h0 g0 h2 g2
	return [expr $a*$h2 + $b*$g2 + $c*$h0 + $d*$g0]
}

proc d4_ig {a b c d} {
	global h1 g1 h3 g3
	return [expr $a*$h3 + $b*$g3 + $c*$h1 + $d*$g1]
}


proc dwt {data} {
	if {[llength $data]%4 != 0} {return}
	set len [llength $data]
	
	lappend data [lindex $data 0] [lindex $data 1]
	
	set s [list]
	set d [list]
	for {set i 0} {$i < $len} {incr i 2} {
		set a0 [lindex $data $i]
		set a1 [lindex $data [expr $i+1]]
		set a2 [lindex $data [expr $i+2]]
		set a3 [lindex $data [expr $i+3]]
		lappend s [d4_h $a0 $a1 $a2 $a3]
		lappend d [d4_g $a0 $a1 $a2 $a3]
	}
	
	return [concat $s $d]
}


proc idwt {data} {
	global hlist glist
	if {[llength $data]%4 != 0} {return}
	set len [llength $data]
	set o [list]
	
	set s [lrange $data 0 [expr $len/2-1]]
	set d [lrange $data [expr $len/2] end]
	set s [concat [lindex $s end] $s]
	set d [concat [lindex $d end] $d]
	
	for {set i 0} {$i<[llength $s]-1} {incr i} {
		set a0 [lindex $s $i]
		set a1 [lindex $d $i]
		set a2 [lindex $s [expr $i+1]]
		set a3 [lindex $d [expr $i+1]]
		lappend o [d4_ih $a0 $a1 $a2 $a3]
		lappend o [d4_ig $a0 $a1 $a2 $a3]
	}
	
	return $o
}


set data [list 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.]
set idata [dwt $data]
set ndata [idwt $idata]

for {set i 0} {$i < [llength $data]} {incr i} {
	puts [format "%10.5f %10.5f %10.5f" \
		[lindex $data $i] [lindex $idata $i] [lindex $ndata $i]]
}

出力結果

  1.00000    2.31079    1.00000
  2.00000    5.13922    2.00000
  3.00000    7.96764    3.00000
  4.00000   10.79607    4.00000
  5.00000   13.62450    5.00000
  6.00000   16.45292    6.00000
  7.00000   19.28135    7.00000
  8.00000   20.59403    8.00000
  9.00000   -0.00000    9.00000
 10.00000   -0.00000   10.00000
 11.00000   -0.00000   11.00000
 12.00000   -0.00000   12.00000
 13.00000   -0.00000   13.00000
 14.00000   -0.00000   14.00000
 15.00000   -0.00000   15.00000
 16.00000   -5.65685   16.00000

↑

コメントをどーぞ†

CategoryTclTk