2014-11-27

私のための離散ウェーブレット変換

離散ウェーブレット変換についてのいくつかのこと
- ウェーブレット変換は、１次元波形をうまく分解する方法の一つ
  - 局所ごとに分解する方法であることが特徴で、全体を分解するフーリエ変換と対照的
- うまく分解するというのは、情報を落とさずに分解するということ
- 多次元に拡張するには、２次元なら「縦→横」、３次元なら「縦→横→高さ」と順繰りに１次元用のウェーブレット変換をすることで対応する
- うまく分解するときには、「AとBとの２つ」
- 「AとBとの２つ」にうまく分けるには、「AとBとは相互に『直交』な関係」であるのがよい
  - そんなよい具合の分け方にはいろいろなものがあり、それをウェーブレットと言う
- この小文では、いろいろなウェーブレットに触れず、一番簡単なHaarウェーブレットのみを扱う。ほかのウェーブレットを使うには、Haarウェーブレットを決める『関数』をほかのウェーブレットに取り換えるだけでよいので、ここでは気にしない
- 変換をしたら、変換結果に対して、さらに同じく「２つに分ける」作業を繰り返す、1->2->4->8と細かくなっていく
和と差とに分解
- もっとも素朴な離散ウェーブレット変換では、数値列の隣接数の差をとる作業Aと和をとる作業Bを行う
- 差と和とは「お互い様」の関係
和・差をとったら、『計算された和・差』を１個飛ばしで採用する(そうすれば、和の列・差の列になったけれど、全部の値の数が増えたりしない)
Rでやってみる
- 和をとるというのは、並ぶ２つの値に1をかけて、それ以外の値には0をかけて、その総和をとることだし、差をとるというのは、並ぶ２つの値の片方に1をかけて、もう片方に-1をかけて、それ以外の値には0をかけて、その総和をとること

# 差・和をとる作業に対応する「フィルタ」f1,f2を作る
f1 <- c(-1,1)
f2 <- c(1,1)
# xは数値ベクトル
# fはフィルタ関数
my.dwt <- function(x,f){
# フィルタをかけた後に、xと同じ長さが出るように頭に0をつける
	x. <- c(rep(0,length(f)-1),x)
	y <- rep(0,length(x))
# フィルタをかけて値を出す
	for(i in 1:length(x)){
		y[i] <- sum(x.[i:(i+length(f)-1)] * f)
	}
# 出た値を１個飛ばしで採用する
	matrix(y,byrow=TRUE,ncol=2)[,2]
}

- 使ってみる

x1 <- c(0,0,0,1,1,0,0,0)
my.dwt(x1,f1)
my.dwt(x1,f2)

- 結果

> x1 <- c(0,0,0,1,1,0,0,0)
> my.dwt(x1,f1)
[1]  0  1 -1  0
> my.dwt(x1,f2)
[1] 0 1 1 0

- ほとんど同じだけれど、セル１個分ずれたベクトルで同じことをしてみる

x2 <- c(0,0,0,0,1,1,0,0)
my.dwt(x2,f1)
my.dwt(x2,f2)

- 結果は前と違う

> x2 <- c(0,0,0,0,1,1,0,0)
> my.dwt(x2,f1)
[1] 0 0 0 0
> my.dwt(x2,f2)
[1] 0 0 2 0

- いわゆる離散ウェーブレットパッケージの結果と比べる

library(waveslim)
dwt(x1,"haar",1)
dwt(x2,"haar",1)

- 結果。"haar"なるフィルタ名で「１段」の処理をすると、２つの数値列が帰ってきて、$d1は「差」、$s1は「和」である。ただし、値は、ちょっと違う。定数倍だけ違う。これは、ウェーブレット変換の情報がうまく合うような補正が入ったから。でも、「そういうものだ」と思っておけばよいのであって、処理の仕方は「和・差」と「定数倍」だ、ということ

(x1,"haar",1)
$d1
[1]  0.0000000  0.7071068 -0.7071068  0.0000000
$s1
[1] 0.0000000 0.7071068 0.7071068 0.0000000
attr(,"class")
[1] "dwt"
attr(,"wavelet")
[1] "haar"
attr(,"boundary")
[1] "periodic"
> dwt(x2,"haar",1)
$d1
[1] 0 0 0 0
$s1
[1] 0.000000 0.000000 1.414214 0.000000
attr(,"class")
[1] "dwt"
attr(,"wavelet")
[1] "haar"
attr(,"boundary")
[1] "periodic"

- 念のため、もっとむちゃくちゃな数値列でもやってみる
  - 定数倍なことはプロットからわかる

x <- runif(100)
my.dwt(x,f1)
my.dwt(x,f2)
dwt(x,"haar",1)

par(mfcol=c(1,2))
plot(my.dwt(x,f1),dwt(x,"haar",1)$d1)
plot(my.dwt(x,f2),dwt(x,"haar",1)$s1)
par(mfcol=c(1,1))