合うか合わないかを数値にする

• フーリエ変換は三角関数と「合うか合わないか」を数値にする
• 三角関数と合うか合わないかは、周波数と位相とをずらして色々試してみればよい
• 色々試してみる、というのは、対象の関数と三角関数との「内積」を取ればよい
  • 関数同士の内積を取るところは「畳み込み」みたいなもの
• 周波数は変えてみないとわからない
• 位相も変えてみないとわからない
  • 位相は、「変えてみないとわからない」けれど、数値化するときに「複素数」に数値化していいのなら、わざわざ位相を変えて試す必要はない
  • という「位相こみこみの関数」と合うか合わないかを試すことで、「位相ずれ」の情報つきで「複素数」になって返ってくるから
• ある周波数について、位相は無視して合っているかを見るには、「返ってきた複素数」の絶対値の大小でわかる
• 合っている位相がいくつかは、その「返ってきた複素数(でその絶対値が大きいもの)」の複素数位相を確認すればよい

• また、同じ周波数の２つの波が合わさったとき、それは、位相と振幅が変わるだけで周波数は変わらないので、周波数スペクトル分解の場合は、同じ周波数の複数の波の合算は、その周波数の「統合的振幅と統合的位相」とともに返ってくる

# 対象とする関数のパラメタ
x <- seq(from=0,to=1,length=1000)*100*pi
# 合わせてみる三角関数のパラメタ
t <- x
# 対象とする関数は「ばっちり三角関数」
f <- sin(x)
# 位相ずれ
diff.t <- seq(from=-1,to=1,length=100)*pi
# 周波数は「ばっちり」な条件で、位相だけずらしてみてみよう
ff <- rep(0,length(diff.t))
for(i in 1:length(diff.t)){
	z <- cos(2*pi*(t+diff.t[i]))+1i * sin(2*pi*(t+diff.t[i]))
	ff[i] <- sum(f*z)
}
# 円になる→絶対値は等しく、位相が違く複素数が返ってきている
plot(ff,type="l")

# 周波数を変えよう
# coef.t <- 0.95^((-10):10)
coef.t <- seq(from=0,to=5,by=0.1)
# 位相は(ゼロ点で)合わせて、周波数だけ変えよう
ff.2 <- rep(0,length(coef.t))
for(i in 1:length(coef.t)){
	z <- cos(2*pi*(coef.t[i]*t))+1i * sin(2*pi*(coef.t[i]*t))
	ff.2[i] <- sum(f*z)

}
plot(ff.2,type="l")
# 位相と周波数とを両方変えよう
diff.coef <- expand.grid(diff.t,coef.t)
ff.full <- rep(0,length(diff.coef[,1]))
for(i in 1:length(ff.full)){
	z <- cos(2*pi*diff.coef[i,2]*t+diff.coef[i,1])+1i * sin(2*pi*diff.coef[i,2]*t+diff.coef[i,1])
	ff.full[i] <- sum(f*z)
}
# 絶対値の大小を表示
image(Mod(matrix(ff.full,ncol=length(coef.t))))
# 返り値の実部を表示
persp(Re(matrix(ff.full,ncol=length(coef.t))))

library(rgl)
# 実部・虚部・絶対値・位相ごとにプロット
plot3d(cbind(diff.coef[,1],diff.coef[,2],Re(ff.full)))
open3d()
plot3d(cbind(diff.coef[,1],diff.coef[,2],Im(ff.full)))
open3d()
plot3d(cbind(diff.coef[,1],diff.coef[,2],Mod(ff.full)))
open3d()
plot3d(cbind(diff.coef[,1],diff.coef[,2],Arg(ff.full)))

• 同一周波数の統合的振幅と統合的位相

x <- seq(from=0,to=1,length=1000)*10*pi
t <- x
library(MCMCpack)
k <- 2
thetas.f <- runif(k)*2*pi
thetas.g <- rep(runif(1),k)
weight.f <- rdirichlet(1,rep(1,k))
weight.g <- rdirichlet(1,rep(1,k))

f <- g <- rep(0,length(x))
for(i in 1:length(thetas.f)){
	f <- f + weight.f[i]*sin(x+thetas.f[i])
	g <- g + weight.g[i]*sin(x+thetas.g[i])
}
matplot(cbind(f,g),type="l")