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

Haar

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

Daubechies

カスケードアルゴリズム

スケーリング関数を描画する時に使う方法で、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で見かけるけど近似計算しているのか何か別の計算方法があるのか謎。解きかたがわかりました。

多重解像度解析

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

リンク

ノイズ除去

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だけ分解と再構築することができます。変換するリストは4の倍数でなければなりません。

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


Attach file: fileD8-plot.gif 350 download [Information] filelevel-2.png 344 download [Information] filelevel-0.png 360 download [Information] filelevel-1.png 375 download [Information]

|New|Edit|Freeze|Diff|Backup|Attach|Copy|Rename|
Last-modified: 2005-02-26 (Sat) 00:00:00 (5014d)
HTML convert time: 0.022 sec.