R　によるランダムフラクタル：タイタンの妖女ここではバカはほめ言葉です：SSブログ

	ブログをはじめるログイン

R　によるランダムフラクタル　[データ解析環境R] [編集]

　ランダムな地形データなどを作成するRの関数を作成した。Rでの複素数の扱いがわからず苦労した。関数の返り値は、N*Nの実行列。ｆｆｔｗ()をコールしているので、rimageパッケージをロードしておく必要がある。

f.ff <- function(N,H=0.85){	# フラクタル次元 : 2-H、条件：0 < H < 1、N は偶数
A <- array(0+0i,c(N,N))
# 左上と右下
for(I in 0:(N/2)) {
	for(J in 0:(N/2)) {
		phase <- 2 * 3.141592 * runif(1)
		if (I != 0 | J != 0) {
			rad <- ((I*I+J*J)^( -(H+1)/2 )) * rnorm(1)
		} else {
			rad <- 0
		}
		A[I+1,J+1] <- complex(1,rad * cos(phase), rad * sin(phase))
		if (I == 0) {
			I0 <- 0
		} else {
			I0 <- N - I
		}
		if (J == 0) {
			J0 <- 0
		} else {
			J0 <- N - J
		}
		A[I0+1,J0+1] <- Conj(A[I+1,J+1])
	}
}
# 虚部は0
A[N/2+1,1]     <- complex(1, Re(A[N/2+1,1])    , 0)
A[1,N/2+1]     <- complex(1, Re(A[1,N/2+1])    , 0)
A[N/2+1,N/2+1] <- complex(1, Re(A[N/2+1,N/2+1]), 0)
# 右上と左下
for(I in 1:(N/2-1)){
	for(J in 1:(N/2-1)){
		phase <- 2 * 3.141592 * runif(1)
		rad <- ((I*I+J*J)^( -(H+1)/2 )) * rnorm(1)
		A[I+1,N-J+1] <- complex(1,rad * cos(phase), rad * sin(phase))
		A[N-I+1,J+1] <- Conj(A[I+1,N-J+1])
	}
}
B <- fftw(A,1)	# 離散逆フーリエ変換
Re(B)			# 虚部のごみを除いた実部を返す
}

　次の本を参考にした。