2013-10-20

準周期と離散フーリエ

まずは１次元離散点をディラックの櫛関数としたうえでその離散フーリエスペクトル化

# 地道にフーリエの式通りに指数関数の複素数乗で計算
my.frie <- function(z,x){
	ret <- rep(0,length(x))
	for(j in 1:length(z)){
		ret <- ret + exp(-2*pi*1i*z[j] * x)
	}
	ret
}
# 結局、実数部分だけがほしいことを思い出せば複素数は不要で三角関数だけ
my.frie.2 <- function(z,x){
	ret <- rep(0,length(x))
	for(j in 1:length(z)){
		ret <- ret + cos(-2*pi*z[j] * x)
	}
	ret
}
# 離散ベクトル同士ならループ要らず
my.frie.3 <- function(z,x){
	tmp <- matrix(z,ncol=1) %*% matrix(x,nrow=1)
	tmp <- cos(-2*pi*tmp)
	apply(tmp,2,sum)
}

いくつかの１次元数列でフーリエスペクトルを取ってみる
- 一点集中と等間隔とは相互に双対ってことか
- 素数が周期性と関係するのは、「素数でない数」は「周期に乗っている数」である、ということの裏返し？？

par(mfcol=c(2,2))
x <- seq(from=-10,to=10,length=1000)
z <- rep(0,100)
y <- my.frie.3(z,x)
plot(x,y,main="１点集中")
z <- seq(from=-10,to=10,length=1000)
y <- my.frie.3(z,x)
plot(x,y,main="等間隔")

library(numbers)
N <- 10000
z <- Primes(N)
x <- seq(from=min(z)*0.1,to=max(z)*10,length=1000)
y <- my.frie.3(z,x)
plot(x,y,main="素数列")
# Prime Powerであるというのは、素因数分解したときの因子が１種類ということなので、1:Nについてそうなっているかの確認をする
powerPrime <- rep(0,N)
for(i in 1:N){
	if(length(table(factorize(i)))==1){
		powerPrime[i] <- 1
	}
}
z <- which(powerPrime==1)
x <- seq(from=min(z)*0.1,to=max(z)*10,length=1000)
y <- my.frie.3(z,x)
plot(x,y,main="素数べき")

par(mfcol=c(1,1))

１次元準結晶の１つの定義の仕方(いろいろ議論があって、この定義でよいのか否か、など色々らしいが)は、点集合をフーリエスペクトル化したときに離散点集合になることらしいので、上記の処理はそれに関すること
さて２次元にする
２(多)次元フーリエ変換は、軸ごとに積分して積み上げる、というのが定義なので重くなりがちだが、離散点(２(多)次元ディラックの櫛関数)ならそうでもない

my.frie.multi <- function(Z,X){
	tmp <- Z %*% t(X)
	tmp <- cos(-2*pi*tmp)
	apply(tmp,2,sum)
}

こちらでcut & projection法による２次元準結晶の例を作ったのでそれに適用してみる
- その２次元準結晶の点の座標をランダムに入れ替えるとどうなるかも見てみる
- まずは準結晶点集合

- そのスペクトル

- そのスペクトルの重み分布

- 座標をばらばらにしたときの点集合

- そのスペクトル

- そのスペクトル重み分布

library(geometry)
a1 <- a2 <- a3 <- a4 <- -20:20
aaaa <- expand.grid(a1,a2,a3,a4)
t1 <- (1+sqrt(5))/2
t2 <- (1-sqrt(5))/2
u <- c(cos(0),sin(0))
v <- c(cos(4*pi/5),sin(4*pi/5))
u. <- c(cos(0),sin(0))
v. <- c(cos(2*pi/5),sin(2*pi/5))

x <- cbind((aaaa[,1]+t1*aaaa[,2])*u[1] + (aaaa[,3]+t1*aaaa[,4])*v[1],(aaaa[,1]+t1*aaaa[,2])*u[2] + (aaaa[,3]+t1*aaaa[,4])*v[2])
y <- cbind((aaaa[,1]+t2*aaaa[,2])*u.[1] + (aaaa[,3]+t2*aaaa[,4])*v.[1],(aaaa[,1]+t2*aaaa[,2])*u.[2] + (aaaa[,3]+t2*aaaa[,4])*v.[2])

s <- which(apply(y^2,1,sum)<1)
plot(x[s,],pch=20,cex=1,xlim=c(-15,15),ylim=c(-15,15))
x<- x[s,]
d.x <- delaunayn(x)
plot(x,pch=20,cex=1,xlim=c(-15,15),ylim=c(-15,15))
d.x.2 <- rbind(d.x[,1:2],d.x[,2:3],d.x[,c(3,1)])
segments(x[d.x.2[,1],1],x[d.x.2[,1],2],x[d.x.2[,2],1],x[d.x.2[,2],2])
Z <- x
xx <- seq(from=min(Z)*0.5, to=max(Z)*2,length=100)
X <- expand.grid(rep(list(xx),length(Z[1,])))
#X <- cbind(xx,xx)
y <- my.frie.multi(Z,X)
image(matrix(y,length(xx),length(xx)))
plot(sort(y))
Z2 <- cbind(sample(Z[,1]),sample(Z[,2]))
plot(Z2,pch=20)
y2 <- my.frie.multi(Z2,X)
image(matrix(y2,length(xx),length(xx)))
plot(sort(y2))