2012-04-01

関数の相関・畳み込み・フーリエ変換・画像処理・処理フィルタ

こちらで時系列データの相関についてコメントしている
関数の似ている程度の評価方法
- である
  - ベクトルの内積の無限長版である
拡張すると相互相関関数(Cross correlation(Wiki))というものがある
- $(f\star g)(t)=\int_{-\infty}^{\infty} f(\tau)g(\tau+t) d\tau$
- タイムラグ $t$ が入って、 $t$ の関数になった
- これが相互相関関数
相互相関関数は畳み込み(Convolution)と呼ばれる計算である(の正負に注意)
- 自己相関(Wiki)
  - $g(\tau)$ の代わりに $f(\tau)$ とすると、 $(f\star f)(t)=\int_{-\infty}^{\infty} f(\tau)f(\tau+t) d\tau$ となる
  - $f(\tau)$ が周期的であると、あるタイムラグ $t$ を持って相互相関関数の値は大きくなる
  - Rで自己相関
処理フィルタのこと
- $\int f(t)g(t)dt$ は $t$ の定義域全体に渡る処理であるが、この世の中の出来事は、( $t$ を位置とすると）、同じ場所での $f(t),g(t)$ の間には関係があり、そのごく近傍( $g(t+\delta t)$ も少し関係があり、それより遠いと関係がない、というようなものが多い
- 画像処理で、ぼやかしたり、くっきりさせたりさせる(こちらでやっているような処理)という処理はまさにこれ(関連記事)
- このようにオリジナルの値を少し改変するような処理において、「フィルターをかける」と呼ぶとき、 $g(t)$ 関数は「フィルター」そのものなので、フィルタリングと言うこともある
畳み込み・相互相関関数とフーリエ変換の関係
- 畳み込みのお話解説
- フーリエ変換は
  - これはととの相関関数(こちらを参照)
    - ちなみに $e^{i 2\pi \tau t}=\cos(2\pi \tau t) - i \sin(2\pi \tau t)$ は無限個の要素を持つ正規直交関数列
    - 『フーリエ変換は，元が互いに正規直交の関係である周期関数列に分解したときに個々の元(関数)につく係数を表している』(こちらを参照)
- また、 $f(t),g(t)$ がフーリエ変換で $F(t),G(t)$ となるとき、 $h(t)=\int_{-\infty}^{\infty} f(\tau)g(\tau+t) dt$ のフーリエ変換 $H(t)$ は $H(t)=F(t)G(t)$ となるという
- このことから、 $h(t)=\int_{-\infty}^{\infty} f(\tau)g(\tau+t) dt$ のときに、（h(t)を出すのが面倒くさいとしたら）、 $F(t),G(t)$ を出して(フーリエ変換して）、 $H(t)=F(t)G(t)$ し、そのうえで、 $H(t)$ を逆フーリエ変換して $h(t)$ を得ることができる(こちら)
- これが、畳み込みの話でフーリエ変換の話が（もれなく）出てくる理由
数列のフーリエ変換は「離散フーリエ変換」(こちら)
さて、Rを使ってみることにする
- １次元データにフィルタをかけてみる
  - フィルタはc(1,2,1)/4 である(フィルタの要素の和は1になっている)
  - これは、自身そのものとその前後とに2:1の重みづけをするフィルタで、結果として、「平滑化」している
    - 平滑化は補間とかとも関係がある(補間についてはこちら)

n <- length(x <- -20:24)
y <- (x-10)^2/1000 + rnorm(x)/8
# c(1,2,1)/4でフィルタをかけている
Han <- function(y) # Hanning
       convolve(y, c(1,2,1)/4, type = "filter")

plot(x,y, main="Using  convolve(.) for Hanning filters")
# x[-c(1,n)]はベクトルxの末端要素を除く処理
lines(x[-c(1  , n)      ], Han(y), col="red")
lines(x[-c(1:2, (n-1):n)], Han(Han(y)), lwd=2, col="dark blue")

- ここで使っているconvolve()関数の中身を見ると、フーリエ変換と逆フーリエ変換が使われていることがわかる

> convolve
function (x, y, conj = TRUE, type = c("circular", "open", "filter")) 
{
    type <- match.arg(type)
    n <- length(x)
    ny <- length(y)
    Real <- is.numeric(x) && is.numeric(y)
    if (type == "circular") {
        if (ny != n) 
            stop("length mismatch in convolution")
    }
    else {
        n1 <- ny - 1
        x <- c(rep.int(0, n1), x)
        n <- length(y <- c(y, rep.int(0, n - 1)))
    }
    x <- fft(fft(x) * (if (conj) 
        Conj(fft(y))
    else fft(y)), inverse = TRUE)
    if (type == "filter") 
        (if (Real) 
            Re(x)
        else x)[-c(1L:n1, (n - n1 + 1L):n)]/n
    else (if (Real) 
        Re(x)
    else x)/n
}

- この関数ソースは自分にとって読みにくいので、書き方を変えてみる

my.convolve <- function(x,y,conj=TRUE,type=c("circular", "open", "filter")){
	type <- match.arg(type)
	n <- length(x)
	ny <- length(y)
	Real <- is.numeric(x) && is.numeric(y)
	if(type == "circular"){
		if(ny !=n){
			stop("length mismatch in convolution")
		}
	}else{
		n1 <- ny-1
		x <- c(rep.int(0,n1),x)
		n <- length(y <- c(y, rep.int(0,n-1)))
	}
	if(conj){
		tmp <- Conj(fft(y))
	}else{
		tmp <- fft(y)
	}
	tmp2 <- fft(x) * tmp
	x <- fft(tmp2,inverse = TRUE)
	ret <- NULL
	if(type == "filter"){
		if(Real){
			ret <- Re(x)[-c(1L:n1, (n - n1 + 1L):n)]/n
		}else{
			ret <- x[-c(1L:n1, (n - n1 + 1L):n)]/n
		}
	}else{
		if(Real){
			ret <- Re(x)/n
		}else{
			ret <- x/n
		}
	}
	ret
}

多項式の畳み込み
- 多項式の掛け算の結果の係数列は、元の多項式の係数列の線形畳み込みになる(Wiki)
- conv()関数で、畳み込んだ多項式がconv.polyに入った

library(pracma)
# c(1,1,1) は x^2 + x + 1 を表し、c(-0.5,1,-1)は-0.5 x^2 + x -1を表す
poly1<-c(1,1,1)
poly2<-c(-0.5,1,-1)
conv.poly<-conv(poly1,poly2)
conv.poly
conv.v<-polyval(conv.poly,3)
conv.v