---
title: "四元数とそのフーリエ変換"
author: "ryamada"
date: "Sunday, December 21, 2014"
output: html_document
---
四元数は、$1,i,j,k$という4つの基本成分で構成される数で
$q = w\times 1 + x \times i + y \times j + z \times k; w,x,y,z \in \mathcal{R}$のように表されるもののことである。ただし、
$i^2=j^2=k^2=-1$である。
複素数の拡張とみなせる。
$1,i,j,k$の積には次のような関係がある。
行列の要素は、(行の要素) x (列の要素)という積を表している。
たとえば、$i \times 1 = i, i \times i = -1, i \times j = k, i\times k = -j$である。
$$\begin{pmatrix} 1 & i & j & k \\ i & -1 & k & -j \\ j & -k & -1 & i \\ k & j & -i & -1 \end{pmatrix}$$
複素数の場合は
$$ \begin{pmatrix} 1 & i \\ i & -1\end{pmatrix}$$
である。
複素数の積行列と四元数の積行列との大きな違いは、複素数のそれが、対称行列であるのに対して、四元数のそれは非対称であることである。
複素数は積行列が対称であるので、二つの複素数$z,W$について、$zw=wz$である。
他方、積行列が非対称であるから、積は順序によって値が異なる。
二つの四元数$p,q$があったとき、一般に、$pq \ne qp$である。
この違いが、四元数の積を用いる諸処理において、右から掛ける処理と左から掛ける処理との二通りを考慮することを要求する。
積の順序が交換できないことは非可換という。
非可換とともにそれ以外の四元数代数の特徴を列挙する。
> 非可換
> 複数の逆元を持つ
> Normed algebra ($||qp|| = ||q||\times ||p||$)
> Idempotent($q^2=q$)、Nilpotent($q^2=0$)が、1,0
のみである。
R のonionパッケージを使ってみる。
onionパッケージでは、複数の四元数を4xn行列のように表示する。
4行は、$1,i,j,k$の各成分に対応する。
また、$H1,Hi,Hj,Hk$ は4つの基本成分を表す。
```{r}
library(onion)
e.q <- quaternion(4)
e.q[1] <- H1
e.q[2] <- Hi
e.q[3] <- Hj
e.q[4] <- Hk
e.q
```
積を計算してみる。
```{r}
for(i in 1:length(e.q)){
print(e.q[i] * e.q)
}
```
適当な四元数を作ってみる。
```{r}
rquat(3) # 正整数係数四元数をランダムに作る関数
```
$(1,i,j,k)$成分は次のように取り扱う。
```{r}
q <- rquat(1)
q
Re(q)
Im(q)
i(q)
j(q)
k(q)
q. <- q
i(q.) <- j(q)
q.
```
任意の四元数$q = w + i x + j y + k z$は
$q = (w + i x) + (y + i z) j$ と表せる。
それは、$ij = k$だから。
このようにすると、$i$を普通の複素数の$i$とみなして、2系列の複素数があって、それらは、$(1,j)$という「実・虚」で分離されている、とみなすことができる。
これを、ortho-split/symplectic formと言う。
また、非可換なので、
$q = (w + i x) + k (z + i y)$
のように、別の$(z+iy)$の前から、純虚四元数を掛けるときには、$k$を掛けることになる。
念のため、Rで確認する。
```{r}
q <- rquat(1)
q
(Re(q) + Hi * i(q)) + (j(q) + Hi * k(q)) * Hj
(Re(q) + Hi * i(q)) + Hk * (k(q) + Hi * j(q))
```
実数成分が0の四元数は純虚四元数。
Hi,Hj,Hkは純虚四元数であって、そのノルムが1。
```{r}
Norm(H1)
Norm(Hi)
```
ノルムが1の四元数は単位四元数。
ノルムが1の任意の純虚四元数を作ってみる。
```{r}
# 任意次元回転行列を作る関数
library(GPArotation)
r <- Random.Start(3)
r
round(r %*% t(r),3)
# この回転行列の列ベクトルを成分とする純虚四元数
M <- quaternion(3)
for(i in 1:3){
Im(M[i]) <- r[,i]
print(M[i])
print(Norm(M[i]))
}
```
また、四元数の3個の単位純虚数は、$ij=k,jk=i,ki=j$という関係にあることにも注意。これは、以下のような任意の単位純虚直交基底四元数3個の間にも成り立つ関係である。
三次元正規直交基底をなす3ベクトルを虚成分とする3つの四元数の積には、$i,j,k$の間と同じ関係がある。
```{r}
round(M[1] * M[2] - M[3],3)
round(M[2] * M[3] - M[1],3)
round(M[3] * M[1] - M[2],3)
```
任意の四元数が$q = w + i x + j y + k z=(w + i x) + (y + i z) j$と$i,j,k$を使って表せたように、
正規直交基底関係にある3つの単位純虚四元数$\mu_1,\mu_2,\mu_3$を使って
$q = w' + \mu_1 x' + \mu_2 y' + \mu_3 z' = (w' + \mu_1 x') + (y' + \mu_1 z) \mu_2$
と表すことができる。
$i,j,k$の係数が、四元数と要素四元数との内積(のようなもの)で決まる様子を、以下に示す。
```{r}
q
q.1 <- Re(q)
q.i <- i(q)
q.j <- j(q)
q.k <- k(q)
q.1. <- Re(g.even(q,H1))
q.i. <- Re(-g.even(Im(q),Im(Hi)))
q.j. <- Re(-g.even(Im(q),Im(Hj)))
q.k. <- Re(-g.even(Im(q),Im(Hk)))
q.1 - q.1. # 内積様計算で算出した係数が一致している
q.i - q.i.
q.j - q.j.
q.k - q.k.
```
適当に定めた「正規直交基底」的3純虚四元数の係数を求める。
```{r}
w. <- Re(q)
x. <- -Re(g.even(Im(q),Im(M[1])))
y. <- -Re(g.even(Im(q),Im(M[2])))
z. <- -Re(g.even(Im(q),Im(M[3])))
q. <- w. * H1 + x. * M[1] + y. * M[2] + z. * M[3]
round(q - q.,3)
```
さらに$q = (w' + \mu_1 x') + (y' + \mu_1 z) \mu_2$を確かめておく。
```{r}
q.. <- (w. + x. * M[1]) + (y. +z. * M[1]) * M[2]
round(q - q..,3)
```
以上より、任意の正規直交基底様四元数のトリオを用いて、ortho-splitができることがわかる。
関数化しておく。
```{r}
my.orthosplit.1 <- function(q,M){
w. <- Re(q)
x. <- -Re(g.even(Im(q),Im(M[1])))
y. <- -Re(g.even(Im(q),Im(M[2])))
z. <- -Re(g.even(Im(q),Im(M[3])))
return(c(w.,x.,y.,z.))
}
my.orthosplit.2 <- function(q,M){
v <- my.orthosplit.1(q,M)
return(list(v=v, M1 = M[1],M2=M[2]))
}
out1 <- my.orthosplit.1(q,M)
sum(out1[1]*H1, out1[2:4]*M)
round(q - sum(out1[1]*H1, out1[2:4]*M),3)
out2 <- my.orthosplit.2(q,M)
out2$v[1]*H1 + out2$v[2]*out2$M1 + (out2$v[3]+out2$v[4]*out2$M1)*out2$M2
round(q - (out2$v[1]*H1 + out2$v[2]*out2$M1 + (out2$v[3]+out2$v[4]*out2$M1)*out2$M2),3)
```
複素数の場合
$z = |z|(\cos{\theta}+i \sin{\theta}])=|z|e^{i\theta}$という表現がある。
四元数では、複素数における$i$を、単位純虚四元数$\mu$で置き換えることで、同様の表現を考えることとする。
これにより$q = |q| (\cos{\theta} + \mu \sin{\theta}) = |q|e^{\mu \theta}$となる。
四元数を実部と虚部とに分け、$|q|$は実部が、$\theta$は虚部の「長さ」が、$i$には四元数の虚部方向の単位純虚四元数が割り当てられる。
$$
|q| = \sqrt{w^2+x^2+y^2+z^2},
\mu = Im(q)/\sqrt{x^2+y^2+z^2},
\theta = arctan{\sqrt{x^2+y^2+z^2}/x}
$$
Rで実施すれば、以下のとおりである。
```{r}
q <- rquat(1)
q.len <- sqrt(Norm(q))
Sq <- Re(q)
Vq <- Im(q)
Vq.len <- sqrt(Norm(Vq))
V.st <- Vq/Vq.len
theta <- atan(Vq.len/Sq)
q
q.len*exp(V.st*theta)
q.len*(cos(theta)+V.st*sin(theta))
```
複素数の場合、$z = |z| (\cos{\theta} + i \sin{\theta})$に対して、$e^z = e^{\cos{\theta}} e^{i\sin{\theta}}$ となるが、このうち、$e^{i\theta}$は$\cos{\theta} + i \sin{\theta}$であるが、この$i \theta$を、単位虚数ベクトルと-1から1までの実数の積と見れば、これを四元数に拡張するときには、ある純虚四元数$p = |p| p_{st}$を単位純虚四元数を使って表せば、$e^{p_{st} |p|}$と考えらえる。
結局、$e^q = e^{Re(q)} e^{Im(q)}=e^{Re(q)}(\cos{|q|} p_{st} \sin{|q|}$となる。
onionパッケージでは四元数の指数関数が実装されている。
```{r}
# 指数
q <- rquat(1)
Sq <- Re(q)
Vq <- Im(q)
Vq.len <- sqrt(Norm(Vq))
V.st <- Vq/Vq.len
exp(Sq) * (cos(Vq.len) + V.st*sin(Vq.len))
# Rには実装されている
exp(q)
# 対数
log(q)
log(sqrt(Norm(q))) + log(q/sqrt(Norm(q)))
q
mu <- Im(q)/sqrt(Norm(Im(q)))
phi <- atan(sqrt(Norm(Im(q)))/Re(q))
sqrt(Norm(q)) * exp(phi * mu)
sqrt(Norm(q)) * (cos(phi) + mu * sin(phi))
```
> 平行移動 $q +p$
> 拡大縮小 全方向：$a \times q$、p方向：$\frac{a-1}{2}(q -pqp) + q$、pを法線とする平面方向：$\frac{a-1}{2}(q + pqp)+q$
> 反転 $pqp$
> 回転 $pq\bar{p}$
```{r}
Q <- quaternion(1)
Im(Q) <- c(1,2,3)
```
平行移動は純虚四元数の加算。
```{r}
s <- quaternion(1)
Im(s) <- c(1,2,3)
Qt <- Q + s
Q
Qt
```
拡大縮小には中心が必要。
拡大縮小は、3通り。
> 指定方向に実数倍(v=v,type="d"):vは単位ベクトル。$\frac{a-1}{2}(q-vqv) + q$
> 指定方向を法線とする面に関して実数倍(v=v,type="p")。$\frac{a-1}{2}(q+vqv) + q$
> 全体として実数倍(v=Null)。$\frac{a-1}{2}(2q) + q=a q$
```{r}
library(GPArotation)
V <- Random.Start(3)
v <- V[,1]
v.q <- quaternion(1)
Im(v.q) <- v
a <- 2
Qs0 <- a * Q
Qs1 <- (1+a)/2 * Q + (1-a)/2*v.q * Q * v.q
Qs2 <- (1+a)/2 * Q - (1-a)/2*v.q * Q * v.q
my.q.scale <- function(Q,m,ctr=c(0,0,0),v=NULL,type=Null){
if(is.null(v)){
ctr.q <- quaternion(1)
return(m* (Q-ctr.q) + ctr.q)
}else{
if(type=="d"){
s <- -1
}else if(type=="p"){
s <- 1
}
}
m.q <- quaternion(1)
Im(m.q) <- v
ctr.q <- quaternion(1)
Im(ctr.q) <- ctr
tmp.Q <- Q-ctr.q
(m-1)/2 * (tmp.Q + s*m.q*tmp.Q*m.q) + tmp.Q + ctr.q
}
Qs0. <- my.q.scale(Q[1],a)
Qs1. <- my.q.scale(Q[1],a,v=v,type="d")
Qs2. <- my.q.scale(Q[1],a,v=v,type="p")
Qs0
Qs0.
Qs1
Qs1.
Qs2
Qs2.
```
### 反転
単位純虚四元数が表す方向を法線とする原点を通る面についての反転は、$vqv$と簡単。
```{r}
Qr <- v.q * Q * v.q
Q
Qr
i(Qr-Q)/i(v.q)
j(Qr-Q)/j(v.q)
k(Qr-Q)/k(v.q)
Norm(Q)
Norm(Qr)
round(g.even(v.q,Q),3)
round(g.even(v.q,Qr),3)
```
```{r}
n.pt <-5000
X <- matrix(rnorm(n.pt*3),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
X[,1] <- X[,1]+X[,2]+ X[,3]^4
X <- X + 1.4
library(rgl)
library(GPArotation)
V <- Random.Start(3)
v <- V[,3]
v <- v/sqrt(sum(v^2))
Xh <- quaternion(n.pt)
for(i in 1:n.pt){
Im(Xh[i]) <- X[i,]
}
vh <- quaternion(1)
Im(vh) <- v
Xhv <- vh * Xh * vh
Xv <- t(as.matrix(Im(Xhv)))
plot3d(rbind(X,rep(min(X),3),rep(max(X),3)))
points3d(Xv[,2:4],col=3)
lines3d(matrix(c(rep(0,3),v),byrow=TRUE,ncol=3),col=2)
F <- cbind(rnorm(n.pt),rnorm(n.pt),rep(0,n.pt)) %*% t(V)
points3d(F,col=4)
```
回転軸を表す単位ベクトルを純虚成分とする四元数$m$と回転角$\phi$を用いて、$q = e^{m \phi}=\cos{\phi} + m \sin{\phi}$を定めれば、$p$のその回転による変換は$qp\bar{q}$である。ただし、$\bar{q}$は共役四元数。
```{r}
phi <- pi/3
v.q
rot.q <- exp(v.q * phi)
rot.q
# 念のため以下の式でも同じことであることを確認
cos(phi) + v.q * sin(phi)
# 四元数共役
rot.q
Conj(rot.q)
Q.r <- rot.q * Q * Conj(rot.q)
# 回転前後でノルムが保存されることの確認
Norm(Q)
Norm(Q.r)
# 回転軸ベクトルと回転前後ベクトルの作る内積が等しいことの確認
round(g.even(v.q,Q),3)
round(g.even(v.q,Q.r),3)
```
```{r}
Xhv2 <- rot.q * Xh * Conj(rot.q)
Xv2 <- t(as.matrix(Im(Xhv2)))
plot3d(rbind(X,rep(min(c(X,Xv2)),3),rep(max(c(X,Xv2)),3)))
points3d(Xv2[,2:4],col=3)
points3d(F,col=4)
lines3d(matrix(c(rep(0,3),v),byrow=TRUE,ncol=3),col=2,lw=5)
```
2つの単位３次元ベクトル$a,b$があるとき、その2ベクトルに垂直なベクトル$m$と$a,b$のなす角$\phi$とには次の関係がある。
$$ba^{-1} = \cos{\phi} + m \sin{\phi}$$.
```{r}
a <- b <- quaternion(1)
Im(a) <- runif(3)
Im(b) <- runif(3)
a <- a/sqrt(Norm(a))
b <- b/sqrt(Norm(b))
lambda <- b*a^{-1}
m <- Im(lambda)
# cos(phi)がa,bの内積(単位ベクトルの内積はcos(角))
phi <- acos(Re(lambda))
-g.even(a,b)
cos(phi)
# aとm、bとmとが垂直
-g.even(a,m)
-g.even(b,m)
```
これを用いると、単位球面上の2点間の大円距離は、角$\phi$であって、
3点,a,b,dがあると、acの大円距離が$\lambda_{a,b}\lambda_{b,d}$となる。
```{r}
a <- b <- d <- quaternion(1)
Im(a) <- runif(3)
Im(b) <- runif(3)
Im(d) <- runif(3)
a <- a/sqrt(Norm(a))
b <- b/sqrt(Norm(b))
d <- d/sqrt(Norm(d))
lambda.ab <- b*a^{-1}
lambda.bd <- d*b^{-1}
lambda.da <- d*a^{-1}
lambda.ab.bd <- lambda.bd*lambda.ab
lambda.da
lambda.ab.bd
```
複素数を使ったフーリエ変換では
$$
F(w) = \int f(t)e^{-iwt}dt
$$
$$
f(t) = \frac{1}{2\pi}\int F(w) e^{iwt}dw
$$
と表せた。
四元数でのフーリエ変換では、$i$の変わりに、単位純虚四元数$\mu$を用いる。
また、四元数の場合には、$f(t)e^{-\mu w t}$と$e^{-\mu wt}f(t)$とが異なるので(非可換)、以下のように、「左から」「右から」の変換の二通りがある。
また、変換と逆変換とで係数を按分して書く。
いわゆる時間・周波数変換で言えば、tが時間wが周波数である。
$$
F(w) = \frac{1}{\sqrt{2\pi}}\int f(t) e^{-\mu wt}dt
$$
$$
f(t) = \frac{1}{\sqrt{2\pi}}\int F(w) e^{\mu wt}dw
$$
関数の積の順序を入れ替える。
$$
F(w) = \frac{1}{\sqrt{2\pi}}\int e^{-\mu wt} f(t)dt
$$
$$
f(t) = \frac{1}{\sqrt{2\pi}}\int e^{\mu wt} F(w)dw
$$
二次元で考えてみる。
$$
e^{-\mathbf{\mu} wx}
$$を
$$
e^{-\mathbf{\mu} (w_1x + w_2y)} = e^{-mu w_1 x} e^{-\mu w_2 y}
$$
のように$(x)$ から$(x,y)$の二軸にするのは、一つの拡張である。
また、四元数の演算の$\mu$は$i$と異なり、方向を変えられるので
$$
e^{-(\mathbf{\mu_1} w_1x + \mathbf{\mu_2}w_2 y)}
$$
のように、軸ごとに単位純虚四元数を変えることもできる。
さらに、四元数では、
$$e^{-(\mathbf{\mu_1} w_1x + \mathbf{\mu_2}w_2 y))} \ne e^{-\mathbf{\mu_1} w_1 x} e^{-\mathbf{\mu_2} w_2 y}
$$
であるので、それも別途定義できる。
また、
$$
e^{-\mathbf{\mu_1} w_1 x} e^{-\mathbf{\mu_2} w_2 y}
$$
のような形は、
$$
e^{-\mathbf{\mu_1} w_1 x} f(.)e^{-\mathbf{\mu_2} w_2 y}
$$
のように
変換対象関数をサンドイッチすることができるので、結局以下のように、各拡張に対して「右」「左」「サンドイッチ」を考慮することができ、それだけ、変換方法のバリエーションが生じる。
> $e^{-\mu(w_1 x + w_2 y)} f(.)$
> $f(.) e^{-\mu(w_1 x + w_2 y)}$
> $e^{-\mu w_1 x } f(.) e^{-\mu w_2 y}$
> $e^{-(\mu_1 w_1 x + \mu_2 w_2 y)} f(.)$
> $f(.) e^{- (\mu_1 w_1 x + \mu_2 w_2 y)}$
> これはサンドイッチできない
> $e^{-\mu_1 w_1 x}e^{- \mu_2 w_2 y} f(.)$
> $f(.) e^{- \mu_1 w_1 x} e^{- \mu_2 w_2 y}$
> $e^{- \mu_1 w_1 x} f(.)e^{- \mu_2 w_2 y}$
四元数フーリエ変換は、上記の定義だが、実装上は、Ortho-splitを使って複素数のフーリエ変換の組み合わせにすることができる。
$q = w' + \mu_1 x' + \mu_2 y' + \mu_3 z' = (w' + \mu_1 x') + (y' + \mu_1 z) \mu_2$
のように分解し、$(w' + \mu_1 x')$に関して$w' + i x'$という複素数のフーリエ変換を実行し、$(y' + \mu_1 z)$に対しても同様に実行し、結果を合わせるときに$\mu_1,\mu_2$を用いて線形和として四元数化する。
ここで$\mu_1$は四元数フーリエ変換をするときの単位純虚四元数である。
まず、Ortho-splitせずに、onionパッケージの四元数関数を用いて、四元数フーリエ変換そのままで実装する。
左掛け・右掛けをオプション選択しつつ、順・逆フーリエもオプション選択する条件で、Ortho-splitして四元数フーリエ変換する関数(1次元)を作る。
```{r}
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)
}
```
この関数を適当な1次元データに適用してみる。
はじめに$\mu=Hi$を使って試した上で、任意の単位純虚四元数で実行する。
データを作る。
```{r}
# データ作成
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
```
左掛けして、逆変換で戻してみる。
```{r}
# 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)
}
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)
```
右掛けを実行してみる。
```{r}
out1 <- out1.inv <- out1.inv.inv <- rquat(length(x))
for(i in 1:length(out1)){
out1[i] <- qft1(x,m,us[i],LR="R")
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)
```
次にOrtho-splitを使った方法を試す。
そのために、複素数のフーリエ変換関数fft()を使う。
fft()の順・逆を使いやすくするためにfft()関数を少し修正した関数を作成しておく。
```{r}
# Rのfftの逆フーリエの1/Nを取り込んだ関数を作っておく
my.fft <- function(x,inverse=FALSE){
k <- 1
if(inverse){
k <- length(x)
}
return(fft(x,inverse=inverse)/k)
}
```
Ortho-splitして、複素数フーリエ変換の線形結合を実行する。
```{r}
# 四元数を二つの複素数に分解
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))
```
四元数フーリエ変換の結果と一致することを確かめる。
```{r}
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)
```
上記のOrtho-splitは左掛けに対応していた。
右掛けのOrtho-splitを実行する。Ortho-splitが変わり、四元数フーリエ変換とが変わっていることに注意する。
```{r}
# 分解の仕方が変更されている
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")
out1.inv[i] <- qft1(x,m,us[i],LR="R",inverse=TRUE)
}
all.equal(out1,out2)
all.equal(out1.inv,out2.inv)
```
四元数の純虚基本要素にて実施したが、任意の単位純虚四元数にて同じことを確認する。
$$
\lambda = qp^{-1} = \cos{\theta} \mu \sin{\theta}
$$
として、$\mu$が$q,p$に垂直であったことを思い出せば、次のようにして、ある単位純虚四元数$m1$によるOrtho-splitが構成できる。
```{r}
# 適当に３次元正規直交基底を作る
m1 <- m2 <- m3 <- quaternion(1)
Im(m1) <- rnorm(3)
m1 <- m1/sqrt(Norm(m1))
m1
Im(m2) <- rnorm(3)
m2 <- m2/sqrt(Norm(m2))
m2
m3 <- Im(m2*m1^{-1})
m3 <- m3/sqrt(Norm(m3))
m2. <- Im(m1*m3^{-1})
m2. <- m2./sqrt(Norm(m2))
g.even(m1,m2.)
g.even(m2.,m3)
g.even(m3,m1)
mu <- m1
nu <- m2.
munu <- m3
munu - mu * nu
```
このようにして作った正規直交基底によって座標変換しつつ、複素数フーリエ変換とその線形和を取る。
```{r}
# この(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")
out1.inv[i] <- qft1(x,m,us[i],LR="L",inverse=TRUE)
}
all.equal(out1,out2)
all.equal(out1.inv,out2.inv)
```
関数化しておく。
順・逆と左右のオプションがある。
x はデータ、mは単位純虚四元数。
```{r}
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")
outL.inv.[i] <- qft1(x,m,us[i],LR="L",inverse=TRUE)
outR.[i] <- qft1(x,m,us[i],LR="R")
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)
```
フィルタリングは近傍シグナルの重みづけ和によって、座標における値を変換することなので、フィルタ関数をかければよい。
四元数では、右掛け・左掛け・サンドイッチ掛けができる点が通常の複素数フィルタリングとの違い。
複素数の畳み込み積分は、フーリエ変換して、処理を二分することができるが、四元数の場合は、その点の挙動が違うことに注意。
```{r}
# アレイフィルタリングのためのユーティリティ関数
my.array.address <- function(a){
d <- dim(a)
L <- list()
L[[1]] <- 1:d[1]
for(i in 2:length(d)){
L[[i]] <- 1:d[i]
}
as.matrix(expand.grid(L))
}
a <- array(0,c(2,3,4))
my.array.address(a)
my.array.loc <- function(ad,d){
d. <- c(1,cumprod(d)[1:(length(d)-1)])
apply((t(ad)-1) * d.,2,sum)+1
}
my.array.expansion <- function(a,d1,d2=d1){
D <- dim(a)
ad <- my.array.address(a)
ad.new <- t(t(ad) + d1)
D.new <- D + d1 + d2
ret <- array(0,D.new)
tmp <- ad.new > 0
tmp2 <-t(t(ad.new) - D.new) <=0
s <- which(apply(tmp*tmp2,1,prod)==1)
loc.new <- my.array.loc(ad.new[s,],D.new)
ret[loc.new] <- a[s]
ret
}
d <- c(2,3,4)
a <- array(1:prod(d),d)
d1 <- rep(2,length(d))
d2 <- rep(3,length(d))
d1 <- c(1,1,1)
d2 <- c(2,4,2)
my.array.expansion(a,d1,d2)
my.array.filter <- function(a,f,ctr=(dim(f)+1)/2){
d.a <- dim(a)
d.f <- dim(f)
diff.d1 <- ctr-1
diff.d2 <- d.f-ctr
ad.f <- my.array.address(f)
ad.f. <- ad.f - ctr
ad.a <- my.array.address(a)
D.new <- d.a + diff.d1 + diff.d2
a.big <- array(0,D.new)
max.loc <- prod(D.new)
for(i in 1:length(ad.f.[,1])){
tmp.ad <- t(t(ad.a) + (-1)*ad.f.[i,]+diff.d1)
tmp.v <- a * f[i]
tmp.loc <- my.array.loc(tmp.ad,D.new)
s <- which(tmp.loc>0 & tmp.loc<max.loc)
a.big[tmp.loc[s]] <- a.big[tmp.loc[s]] + tmp.v[s]
}
my.array.expansion(a.big,-diff.d1,-diff.d2)
}
d <- c(2,3,4)
a <- array(1:prod(d),d)
f <- array(1,rep(3,length(d)))
a. <- my.array.filter(a,f)
a. - a
my.array.filter.q <- function(a,f,a.v = c(a),f.v = c(f),ctr=(dim(f)+1)/2,LR="L"){
d.a <- dim(a)
d.f <- dim(f)
diff.d1 <- ctr-1
diff.d2 <- d.f-ctr
ad.f <- my.array.address(f)
ad.f. <- ad.f - ctr
ad.a <- my.array.address(a)
D.new <- d.a + diff.d1 + diff.d2
a.big <- array(1:prod(D.new),D.new)
a.big.v <- quaternion(length(a.big))
max.loc <- prod(D.new)
for(i in 1:length(ad.f.[,1])){
tmp.ad <- t(t(ad.a) + (-1)*ad.f.[i,]+diff.d1)
if(LR=="L"){
tmp.v <- f.v[i] * a.v
}else if(LR=="R"){
tmp.v <- a.v * f.v[i]
}
tmp.loc <- my.array.loc(tmp.ad,D.new)
s <- which(tmp.loc>0 & tmp.loc<max.loc)
a.big.v[tmp.loc[s]] <- a.big.v[tmp.loc[s]] + tmp.v[s]
}
tmp.ad <- my.array.expansion(a.big,-diff.d1,-diff.d2)
a.big.v[c(tmp.ad)]
}
```
```{r}
# zは複素数のベクトル
# int0,int1はintensityの上下限、sat0,sat1はSaturation(彩度)の上下限
my.hsv <- function(z,int0=0.6,sat0=0.3,int1=1,sat1=1){
# 複素数の偏角
arg <- Arg(z)
s <- which(arg<0)
arg[s] <- arg[s]+2*pi
# 複素数の絶対値
r <- Mod(z)
# 絶対値が非常に大きくてもそこそこの色になるように対数変換
s <- which(r>1)
r[s] <- log(r[s])
# 絶対値で周期性が出るように4のmod
r. <- 4*(r%%1)
k <- floor(r.)
r. <- r.-k
# 明度が上限、明度が下限、彩度が上限、彩度が下限の４パターンを
# 4のmodに対応づける
# 明度・彩度を動かすときは、複素数の絶対値で１次線形変化
inten <- sat <- rep(0,length(r))
s <- which(k==0)
inten[s] <- int1
sat[s] <- sat1-(sat1-sat0)*r.[s]
s <- which(k==1)
inten[s] <- int1-(int1-int0)*r.[s]
sat[s] <- sat0
s <- which(k==2)
inten[s] <- int0
sat[s] <- sat1-(sat1-sat0)*(1-r.[s])
s <- which(k==3)
inten[s] <- int1-(int1-int0)*(1-r.[s])
sat[s] <- sat1
return(cbind(arg,inten,sat))
}
my.hsv2rgb <- function(h,s,v){
hi <- floor(h/(2*pi)*6)
hi[which(hi==6)] <- 0
f <- (h/(2*pi)*6) %%1
p <- v*(1-s)
q <- v *(1-f*s)
t <- v *(1-(1-f)*s)
r <- g <- b <- rep(0,length(h))
s <- which(hi==0)
r[s] <- v[s];g[s] <- t[s]; b[s] = p[s];
s <- which(hi==1)
r[s] <- q[s];g[s] <- v[s]; b[s] = p[s];
s <- which(hi==2)
r[s] <- p[s];g[s] <- v[s]; b[s] = t[s];
s <- which(hi==3)
r[s] <- p[s];g[s] <- q[s]; b[s] = v[s];
s <- which(hi==4)
r[s] <- t[s];g[s] <- p[s]; b[s] = v[s];
s <- which(hi==5)
r[s] <- v[s];g[s] <- p[s]; b[s] = q[s];
return(cbind(r,g,b))
}
x <- seq(from=-4,to=4,len=2^8)
xx <- expand.grid(x,x)
z <- xx[,1]+1i * xx[,2]
my.f <- function(z){
(z^2-1)*(z-2-1i)^2/(z^2+2+2*1i)
}
w <- my.f(z)
hsv <- my.hsv(w,int0=0.1,sat0=0.1,int1=1,sat1=1)
col <- my.hsv2rgb(hsv[,1],hsv[,3],hsv[,2])
plot(xx,pch=20,col=rgb(col[,1],col[,2],col[,3]))
```
```{r}
col.q <- quaternion(length(col[,1]))
i(col.q) <- col[,1]
j(col.q) <- col[,2]
k(col.q) <- col[,3]
a.q <- col.q
q <- Hi
f.q <- quaternion(9)
Re(f.q) <- rep(c(1,0,Re(q)),3)
i(f.q) <- rep(c(0,0,i(q)),3)
j(f.q) <- rep(c(0,0,j(q)),3)
k(f.q) <- rep(c(0,0,k(q)),3)
f.q. <- Conj(f.q)
a=matrix(0,length(x),length(x))
f=matrix(0,3,3)
out <- my.array.filter.q(a,f,a.v = a.q,f.v = f.q,ctr=(dim(f)+1)/2,LR="L")
out2 <- my.array.filter.q(a,f,a.v = out,f.v = f.q.,ctr=(dim(f)+1)/2,LR="R")
out2 <- out2/6
R <- i(out2)
G <- j(out2)
B <- k(out2)
R <- (R-min(R))/(max(R)-min(R))
G <- (G-min(G))/(max(G)-min(G))
B <- (B-min(B))/(max(B)-min(B))
plot(xx,pch=20,col=rgb(R,G,B))
dev.new()
plot(xx,pch=20,col=rgb(col[,1],col[,2],col[,3]))
```
ウェーブレット変換は、ハイパス・ローパスの２フィルタをかけて、情報を分解する仕組みなので、フィルタが掛かれば、ウェーブレット変換もできる。
自己相関・相互相関は、ずれを入れつつ、関数を掛け合わせるので、その離散版も同様に定義できる。
積を使うので、四元数の場合には、順序違いが生じることは、フーリエ変換・畳み込みと同じ。
通常の複素数の関数相関の場合には、畳み込みと似ているが共役化した関数を畳み込みに使う、という関係にあるが、四元数でも、畳み込みと相関とには共役などが役割を持つ。
また、四元数関数を単位純虚四元数方向と、それに垂直な成分とに分離して考えると、通常の関数相関のフーリエ変換のように、二つのフーリエ変換の積としての表現が可能で、最終結果は、その線形和になる。
また、得られた値は、amplitudeとphaseとでとらえることができ、phaseは通常の複素数相関の場合はスカラーだが四元数の場合は多次元になる。
