2013-12-09

Chapter 3 UDWT step by step　ぱらぱらめくる『Conceptual Wavelets in Digital Signal Procesing』

haarフィルタをかけて細かく見るのと同時に、(1,1,...)というようなただ足し合わせるフィルタをかけて滑らかにするのとを同時にやるのがUDWT
細かく見る方は逆フィルタをかけ直し、滑らかにする方は、さらに追い打ちで同じ(1,1,...)フィルタをかける
出た結果を足し合わせると、元のデータの情報をすべて持っている
確かめてみる

S <- x
# haarフィルタ
cD1 <- conv(make.haar(2),S)
cD1
# 逆haarフィルタ
d1 <- conv(make.haar(2)[2:1],cD1)
d1
# (1,1)フィルタ
cA1 <- conv(rep(1,2),S)
cA1
# 追い打ちでもう一度(1,1)フィルタ
a1 <- conv(rep(1,2),cA1)
a1
# 足し合わせる
round(d1 + a1)

> round(d1 + a1)
 [1]   0 320 320 320 320   0   0   0   0   0
> S <- x
> cD1 <- conv(make.haar(2),S)
> round(cD1)
[1] -80   0   0   0  80   0   0   0   0
> d1 <- conv(make.haar(2)[2:1],cD1)
> round(d1)
 [1] -80  80   0   0  80 -80   0   0   0   0
> cA1 <- conv(rep(1,2),S)
> round(cA1)
[1]  80 160 160 160  80   0   0   0   0
> a1 <- conv(rep(1,2),cA1)
> round(a1)
 [1]  80 240 320 320 240  80   0   0   0   0
> round(d1 + a1)
 [1]   0 320 320 320 320   0   0   0   0   0

最後の足し合わせがまさに、もとのデータになるためには、長さの調節と値の標準化が必要
- 長さの調節にはwkeep関数
- 値の標準化には、前章でやったように出てくる値を割る、という手もあるが、フィルタの値を調整するのもよい
  - (-1,1)を(-0.5,0.5)にする:hid
  - (1,-1)を(0.5,-0.5)にする:hir
  - (1,1)を(0.5,0.5)にする:lod
  - 追い打ちの(1,1)を(0.5,0.5)にする:lor

S <- x
cD1 <- conv(make.haar(2)/2,S)
round(cD1)
d1 <- conv(make.haar(2)[2:1]/2,cD1)
round(d1)
cA1 <- conv(rep(1,2)/2,S)
round(cA1)
a1 <- conv(rep(1,2)/2,cA1)
round(a1)
round(d1 + a1)
S. <- d1 + a1
round(my.wkeep.v(S.,length(S)))

> cD1 <- conv(make.haar(2)/2,S)
> round(cD1)
[1] -40   0   0   0  40   0   0   0   0
> d1 <- conv(make.haar(2)[2:1]/2,cD1)
> round(d1)
 [1] -20  20   0   0  20 -20   0   0   0   0
> cA1 <- conv(rep(1,2)/2,S)
> round(cA1)
[1] 40 80 80 80 40  0  0  0  0
> a1 <- conv(rep(1,2)/2,cA1)
> round(a1)
 [1] 20 60 80 80 60 20  0  0  0  0
> round(d1 + a1)
 [1]  0 80 80 80 80  0  0  0  0  0
> S. <- d1 + a1
> round(my.wkeep.v(S.,length(S)))
[1] 80 80 80 80  0  0  0  0

フィルタを連続してかける(畳みこむ)*1のは、フィルタの畳みこみをしてから、かけて(畳みこんで)も同じこと

round(conv(conv(c(-1,1)/2,c(1,-1)/2),S))
round(d1)
round(conv(conv(c(1,1)/2,c(1,1)/2),S))
round(a1)

> round(conv(conv(c(-1,1)/2,c(1,-1)/2),S))
 [1] -20  20   0   0  20 -20   0   0   0   0
> round(d1)
 [1] -20  20   0   0  20 -20   0   0   0   0
> round(conv(conv(c(1,1)/2,c(1,1)/2),S))
 [1] 20 60 80 80 60 20  0  0  0  0
> round(a1)
 [1] 20 60 80 80 60 20  0  0  0  0

こうしてみてくると、いかに詳細化と大雑把化の２つのトラックに関して、情報を失わないように分離して畳みこむかのルールを決めて、そこの畳みこみルールの部分(フィルタ)を先に畳みこむことが、「変換」の本質であるとわかる
細かくする方を１段階、粗くする方を２段階にするとき、このフィルタは

f.1 <- conv(c(-1,1),c(1,-1))
# 以下、f.2.1,f.2.2にはconv()のかけ方が違うが結果は同じ２つのやり方があるらしい
#f.2.1 <- conv(conv(c(1,1),c(1,0,-1)),conv(c(1,1),c(-1,0,1)))
f.2.1 <- conv(conv(c(1,1),c(1,0,-1)),conv(c(-1,0,1),c(1,1)))
#f.2.2 <- conv(conv(c(1,1),c(1,0,1)),conv(c(1,1),c(1,0,1)))
f.2.2 <- conv(conv(c(1,1),c(1,0,1)),conv(c(1,0,1),c(1,1)))
f.1
f.2.1
f.2.2

> f.1
[1] -1  2 -1
> f.2.1
[1] -1 -2  1  4  1 -2 -1
> f.2.2
[1] 1 2 3 4 3 2 1

Moving averager
- c(1,1,...)で畳みこむと、その幅の範囲での平均値を返していることになる
- これを称して、「畳みこみ関数がMoving」して、「averageを計算する仕事」をしている、と言う

n <- 10
X <- matrix(0,n,length(x))
for(i in 1:6){
	X[i,] <- my.wkeep.v(conv(rep(1,i)/i,x),length(x))
}

matplot(t(X),type="l")

*1:-1,1)->(1,-1