コードする 3 Quaternion Fourier Transforms ぱらぱらめくる『Quaternion Fourier Transforms for Signal and Image Processing』

• コードするとしたら、四元数をsimplex/perplexという直交する２つの複素数で分けて、それぞれに、普通のフーリエ変換のコードを利用するのが「吉」。なぜなら、普通の高速フーリエ変換のアルゴリズムはとても速いから
• それがうまくいっていることを、ベタに四元数で書いた結果と確認すれば、さらによくわかる
• やってみる
• [q = w + ix + jy + kz = (w+ix) + (y+iz)j = z1 + z2 i]と分解する。ここで、z2をjの左側にくくりだしている場合が、Left-side quaternion fourierに相当するRight-sideではと分解する
• 四元数離散フーリエを、LR(左右)、順逆の４通りで関数化する
  • がLeft-sidedの順
  • がRight-sidedの順
  • がLeft-sidedの逆
  • がLeft-sidedの逆

qft1 <- function(x,m,u,inverse=FALSE,LR="L"){
	n <- length(x)
	s <- -1
	k <- 1
	if(inverse){
		s <- 1
		k <- 1/n
	}
	if(LR=="L"){
		ret <- k*sum(exp(s*m * 2*pi*(0:(n-1))*u/n)*x)
	}else if(LR=="R"){
		ret <- k*sum(x*exp(s*m * 2*pi*(0:(n-1))*u/n))
	}
	return(ret)
}

• 適当に四元数の数列を作り、Left-sided, Right-sidedをやってみる

# データ作成
n.pt <- 2^4
t <- seq(from=0,to=1,length=n.pt)
x <- rquat(n.pt)
Re(x) <- sin(t) + rnorm(n.pt)*0.01
i(x) <- cos(t+2) + t + rnorm(n.pt)*0.01
j(x) <- sin(4*t+2) + rnorm(n.pt)*0.01
k(x) <- 3 + rnorm(n.pt)*0.01
# Left-sided
m <- Hi # ひとまず、この基本純虚四元数について実施する
us <- (0:(length(x)-1))
out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
	out1[i] <- qft1(x,m,us[i],LR="L") # Left-sided
	out1.inv[i] <- qft1(x,m,us[i],inverse=TRUE)

}
# inverseでもとに戻る
for(i in 1:length(out1)){
	out1.inv.inv[i] <- qft1(out1,m,us[i],LR="L",inverse=TRUE)
}
all.equal(x,out1.inv.inv)
# Right-sided
out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
	out1[i] <- qft1(x,m,us[i],LR="R") # Left-sided
	out1.inv[i] <- qft1(x,m,us[i],LR="R",inverse=TRUE)
}
for(i in 1:length(out1)){
	out1.inv.inv[i] <- qft1(out1,m,us[i],LR="R",inverse=TRUE)
}
all.equal(x,out1.inv.inv)

• さて、これは面倒くさい、ということで、普通の複素数フーリエ変換に分解することをやってみる
• と分解する(Left-side)

# Rのfftの逆フーリエの1/Nを取り込んだ関数を作っておく
my.fft <- function(x,inverse=FALSE){
	k <- 1
	if(inverse){
		k <- length(x)
	}
	return(fft(x,inverse=inverse)/k)
}
# 四元数を二つの複素数に分解
x1 <- Re(x) + 1i * i(x)
x2 <- j(x) + 1i * k(x)
# それぞれをフーリエする
tmp.out1 <- my.fft(x1)
tmp.out2 <- my.fft(x2)
tmp.out1.inv <- my.fft(x1,inverse=TRUE)
tmp.out2.inv <- my.fft(x2,inverse=TRUE)
# 四元数に戻す
out2 <- Re(tmp.out1) + Hi * Im(tmp.out1) +  Hj * Re(tmp.out2) + Hk * Im(tmp.out2)
out2.inv <- Re(tmp.out1.inv) + Hi * Im(tmp.out1.inv) + Hj * Re(tmp.out2.inv) + Hk * Im(tmp.out2.inv)
# 四元数フーリエ、Left-sided
m <- Hi
us <- (0:(length(x)-1))

out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
	out1[i] <- qft1(x,m,us[i],LR="L") # Left-sided
	out1.inv[i] <- qft1(x,m,us[i],LR="L",inverse=TRUE)
}
all.equal(out1,out2)
all.equal(out1.inv,out2.inv)

• と分解する(Left-side)

# 分解の仕方が変更されている
x1 <- Re(x) + 1i * i(x)
x2 <- k(x) + 1i * j(x)
tmp.out1 <- my.fft(x1)
tmp.out2 <- my.fft(x2)
tmp.out1.inv <- my.fft(x1,inverse=TRUE)
tmp.out2.inv <- my.fft(x2,inverse=TRUE)
# 四元数への戻し方が変更されている
out2 <- Re(tmp.out1) + Hi * Im(tmp.out1) + Hk * Re(tmp.out2) + Hj * Im(tmp.out2)
out2.inv <- Re(tmp.out1.inv) + Hi * Im(tmp.out1.inv) + Hk * Re(tmp.out2.inv)  + Hj * Im(tmp.out2.inv)
m <- Hi
us <- (0:(length(x)-1))

out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
	out1[i] <- qft1(x,m,us[i],LR="R") # Right-sided
	out1.inv[i] <- qft1(x,m,us[i],LR="R",inverse=TRUE)
}
all.equal(out1,out2)
all.equal(out1.inv,out2.inv)

• さて、Hiを純虚単位四元数とする場合の四元数フーリエが普通の複素数フーリエ２つの和に分解できることを見た
• 任意の純虚四元数とそれに直交する純虚四元数とを使ったortho-splitにこれを拡張する
  • という基底をという基底に換えてやればOK

# 適当に３次元正規直交基底を作る
library(GPArotation)
R <- Random.Start(3)
# それにより純虚四元数トリオを作る
mu <- nu <- munu <- quaternion(1)
Im(mu) <- R[,1]
Im(nu) <- R[,2]
Im(munu) <- R[,3]
munu - mu * nu
# この(mu,nu,munu)系に座標変換
# この(mu,nu,munu)系に座標変換
W <- Re(x)
X <- -Re(g.even(Im(x),Im(mu)))
Y <- -Re(g.even(Im(x),Im(nu)))
Z <- -Re(g.even(Im(x),Im(munu)))
# 一致の確かめ
XX <- W * H1 + X * mu + Y * nu + Z * munu
round(XX - x,3)
XXX <- (W + X * mu) + (Y +Z * mu) * nu
round(XXX - x,3)
# 新座標 q = z1 + ze nu にて、２複素数に分解
x.1 <- W + 1i * X
x.2 <- Y + 1i * Z

# 複素数２系列でフーリエ変換
tmp.out1 <- my.fft(x.1)
tmp.out2 <- my.fft(x.2)
tmp.out1.inv <- my.fft(x.1,inverse=TRUE)
tmp.out2.inv <- my.fft(x.2,inverse=TRUE)
# 四元数基底に戻す
out2 <- Re(tmp.out1) + mu * Im(tmp.out1) + nu * Re(tmp.out2) + munu * Im(tmp.out2)
out2.inv <- Re(tmp.out1.inv) + mu * Im(tmp.out1.inv) + nu * Re(tmp.out2.inv) + munu * Im(tmp.out2.inv)
# 四元数でじかにフーリエ変換
m <- mu
us <- (0:(length(x)-1))

out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
  out1[i] <- qft1(x,m,us[i],LR="L") # Left-sided
	out1.inv[i] <- qft1(x,m,us[i],LR="L",inverse=TRUE)
}
# ２方法の結果は一致
all.equal(out1,out2)
all.equal(out1.inv,out2.inv)

• 関数化

qft1.comp <- function(x,m,inverse=FALSE,LR="L"){
  # mを含む正規直交基底4元数を作成
  m1 <- m
	m2 <- m3 <- quaternion(1)
	Im(m2) <- rnorm(3)
	m2 <- m2/sqrt(Norm(m2))
	m3 <- Im(m2*m1^{-1})
	m3 <- m3/sqrt(Norm(m3))
	m2. <- Im(m1*m3^{-1})
	m2. <- m2./sqrt(Norm(m2))
	mu <- m1
	nu <- m2.
	munu <- m3
	# この(mu,nu,munu)系に座標変換
	W <- Re(x)
	X <- -Re(g.even(Im(x),Im(mu)))
	Y <- -Re(g.even(Im(x),Im(nu)))
	Z <- -Re(g.even(Im(x),Im(munu)))
	# 新座標 q = z1 + ze nu にて、２複素数に分解

	if(LR=="L"){
		x.1 <- W + 1i * X
		x.2 <- Y + 1i * Z
	}else if(LR=="R"){
		x.1 <- W + 1i * X
		x.2 <- Z + 1i * Y
	}

	# 複素数２系列でフーリエ変換
	tmp.out1 <- my.fft(x.1,inverse=inverse)
	tmp.out2 <- my.fft(x.2,inverse=inverse)
	# 四元数基底に戻す
	if(LR=="L"){
		out2 <- Re(tmp.out1) + mu * Im(tmp.out1) + nu * Re(tmp.out2) + munu * Im(tmp.out2)
	}else if(LR=="R"){
		out2 <- Re(tmp.out1) + mu * Im(tmp.out1) + munu * Re(tmp.out2) + nu * Im(tmp.out2)
	}
	out2
}
```

```{r}
m <- mu
outL <- qft1.comp(x,m,inverse=FALSE,LR="L")
outR <- qft1.comp(x,m,inverse=FALSE,LR="R")

outL.inv <- qft1.comp(x,m,inverse=TRUE,LR="L")
outR.inv <- qft1.comp(x,m,inverse=TRUE,LR="R")
us <- (0:(length(x)-1))

outL. <- outR. <- outL.inv. <- outR.inv. <- quaternion(length(x))
for(i in 1:length(outL.)){
  outL.[i] <- qft1(x,m,us[i],LR="L") # Left-sided
  outL.inv.[i] <- qft1(x,m,us[i],LR="L",inverse=TRUE)
  outR.[i] <- qft1(x,m,us[i],LR="R") # Left-sided
  outR.inv.[i] <- qft1(x,m,us[i],LR="R",inverse=TRUE)
}
# ２方法の結果は一致
round(outL-outL.,3)
round(outL.inv-outL.inv.,3)
round(outR-outR.,3)
round(outR.inv-outR.inv.,3)