2013-12-09

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

都合8回のテストの点が、80,80,80,80,0,0,0,0という例(途中でやる気を失った)
Haar フィルタというのは[1,-1]なるフィルタ。これを使うと「差分の検出」をしていることに他ならない

x <- c(rep(80,4),rep(0,4))
x
-diff(c(0,x,0))

> x <- c(rep(80,4),rep(0,4))
> x
[1] 80 80 80 80  0  0  0  0
> -diff(c(0,x,0))
[1] -80   0   0   0  80   0   0   0   0

また、畳みこみとみなしてもよい
- ただし、畳みこみで(-1,1)を使うというのは、haarフィルタとして(1,-1)を使う、ということであることに注意

library(pracma)
haar.1 <- c(-1,1)
round(conv(haar.1,x))

> round(conv(c(-1,1),x))
[1] -80   0   0   0  80   0   0   0   0

ただし、WTではベクトルの長さが変わらないのが望ましい。そのために、wkeepという関数がMATLABやWTでは使われるそうだ(ここにあるように)
- Rには対応する関数が無いようなので作っておく

my.wkeep.v <- function(x,L,type="c"){
	n <- length(x)
	m <- n-L
	ret <- c()
	if(type=="c"){
		R <- L <- 0
		tmp <- m%%2
		R <- (m-tmp)/2
		L <- (m+tmp)/2
		ret <- x[(R+1):(n-L)]
	}else if (type=="r"){
		ret <- x[-(1:m)]
	}else if(type=="l"){
		ret <- x[1:L]
	}
	ret
}

conv.out <- round(conv(haar.1,x))
conv.out
my.wkeep.v(conv.out,length(x))

> conv.out <- round(conv(haar.1,x))
> conv.out
[1] -80   0   0   0  80   0   0   0   0
> my.wkeep.v(conv.out,length(x))
[1] -80   0   0   0  80   0   0   0

haarフィルタを格上げしていく。(-1,1)から(-1,0,1),(-1,-1,1,1)へ

haar.1 <- c(-1,1)
round(conv(haar.1,x))
haar.2 <- c(-1,0,1)
conv.out.2 <- round(conv(haar.2,x))
my.wkeep.v(conv.out.2,length(x))
haar.3 <- c(-1,-1,1,1)
conv.out.3 <- round(conv(haar.3,x))
my.wkeep.v(conv.out.3,length(x))

> haar.1 <- c(-1,1)
> round(conv(haar.1,x))
[1] -80   0   0   0  80   0   0   0   0
> haar.2 <- c(-1,0,1)
> conv.out.2 <- round(conv(haar.2,x))
> my.wkeep.v(conv.out.2,length(x))
[1] -80   0   0  80  80   0   0   0
> haar.3 <- c(-1,-1,1,1)
> conv.out.3 <- round(conv(haar.3,x))
> my.wkeep.v(conv.out.3,length(x))
[1] -160  -80    0   80  160   80    0    0

このhaarフィルタを任意の長さで作ってやれば

make.haar <- function(n){
	tmp <- n%%2
	tmp2 <- (n-tmp)/2
	ret <- c(rep(-1,tmp2),rep(0,tmp),rep(1,tmp2))
	ret
}
for(i in 1:10){
	print(make.haar(i))
}

> for(i in 1:10){
+ print(make.haar(i))
+ }
[1] 0
[1] -1  1
[1] -1  0  1
[1] -1 -1  1  1
[1] -1 -1  0  1  1
[1] -1 -1 -1  1  1  1
[1] -1 -1 -1  0  1  1  1
[1] -1 -1 -1 -1  1  1  1  1
[1] -1 -1 -1 -1  0  1  1  1  1
[1] -1 -1 -1 -1 -1  1  1  1  1  1

その上でWTする

for(i in 1:10){
	tmp.haar <- make.haar(i)
	tmp.conv.out <- round(conv(tmp.haar,x))
	print(my.wkeep.v(tmp.conv.out,length(x)))
}

> for(i in 1:10){
+ tmp.haar <- make.haar(i)
+ tmp.conv.out <- round(conv(tmp.haar,x))
+ print(my.wkeep.v(tmp.conv.out,length(x)))
+ }
[1] 0 0 0 0 0 0 0 0
[1] -80   0   0   0  80   0   0   0
[1] -80   0   0  80  80   0   0   0
[1] -160  -80    0   80  160   80    0    0
[1] -160  -80   80  160  160   80    0    0
[1] -240 -160    0  160  240  160   80    0
[1] -240  -80   80  240  240  160   80    0
[1] -320 -160    0  160  320  240  160   80
[1] -240  -80   80  240  320  240  160   80
[1] -320 -160    0  160  320  320  240  160

さらにこれを各段の値を標準化すると

nn <- 10
M <- matrix(0,nn,length(x))

for(i in 1:nn){
	tmp.haar <- make.haar(i)
	tmp.conv.out <- conv(tmp.haar,x)
	M[i,] <- my.wkeep.v(tmp.conv.out,length(x))
	M[i,] <- M[i,]/sqrt(i)
}

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

ここで、最も谷底と山頂とが大きな値をもつのは第８段。そのときの畳みこみは(-1,-1,-1,-1,1,1,1,1)。これに対応するhaarフィルタは(1,1,1,1,-1,-1,-1,-1)で、この形は、(80,80,80,80,0,0,0,0)と一番よく似ている