勉強中のメモとかポインタです。
#contents
***Haar [#i6752c3c]
足して2で割る。または√2で割る。Waveletの紹介でよく紹介される方法は足して２で割っているが、正規化したデータを使う解きは√2の方を用いるようだ？
***Daubechie [#t8b4d113]
-DaubechieのWaveletはHaarを拡張したもので、D2の時はHaarの基底と同じ。D4、D6、D8…と長めのフィルターがある。プログラムするときはすでに係数を計算してあるのを用意して使う。D20くらいまではWebで検索すると見つけられるようです。分解と再構築で同じ係数を使う。~
***Daubechies [#t8b4d113]
-DaubechiesのWaveletはHaarを拡張したもので、D2の時はHaarの基底と同じ。D4、D6、D8…と長めのフィルターがある。プログラムするときはすでに係数を計算してあるのを用意して使う。D20くらいまではWebで検索すると見つけられるようです。分解と再構築で同じ係数を使う。~
-D2の時はともかく端のところで折り返さないとならないので、そのへんが微妙に不安な気がする…。
-ドベシーおばさんの若かりし頃→[[写真:http://www.math.washington.edu/~sheetz/Milliman/Daubechies.jpg]]~
小学校で担任だった先生に激似。
***カスケードアルゴリズム [#i81557dd]
スケーリング関数を描画する時に使う方法で、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のプログラム>#tclcascade]]に食わせて分割処理。分割はtを2倍した一個上の解像度の値を元に計算できるので、順次解像度を上げながらループさせればＯＫ（？）。~
マザーウェーブレットは計算したスケーリング関数の値を使って
ψ(t) = -h(0)φ(2t-1) + h(1)φ(2t) - h(2)φ(2t+1) + h(3)φ(2t+2)
で簡単に求めることができる。
-D8(N=8)をMuPADで固有ベクトルを計算して、Tclで補間処理して、gnuplotで描画してみた例
#ref(D8-plot.gif,nolink)
***多重解像度解析 [#dabdad05]
Mallat変換とかピラミッドアルゴリズムとか呼ばれる方法は分離された差分には手をつけないで、縮小画像に再びウェーブレット変換していくという方法で計算量がO(NlogN)になる。高速ウェーブレット変換(FWT)などとも呼ばれているような気がしないでもないが、FWTはO(N)で計算できるとか書いているページもあったりしてなんだか用語が曖昧でよくわかりません。~
画像へWavelet変換を適用する場合、縦横と折りたたむように変換する。次の画像は自作の変換プログラムでHaarのWavelet変換をレベル-2まで適用してみた例。
-レベル-0（何もしない）
#ref(level-0.png,nolink)
-レベル-1
#ref(level-1.png,nolink)
-レベル-2
#ref(level-2.png,nolink)
***リンク [#oc34d761]
-http://www.wavelet.org/
-http://www.bath.ac.uk/elec-eng/pages/sipg/resource/warehouse.htm
-http://bell.kuee.kyoto-u.ac.jp/~yoneyone/wavelet/
***ノイズ除去 [#b7cd1b87]
-自作ライブラリによるノイズ除去プログラム->AntiNoise
***Tclによるカスケード分割の書き捨てプログラム[#tclcascade]
こんなんでいいのかね・・・。まあいいや・・・。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変換逆変換の書き捨てプログラム [#a9f86908]
練習を兼ねて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
***コメントをどーぞ [#ncaf6c97]
#comment
----
[[CategoryTclTk]]
HTML convert time: 0.002 sec.