今回は，Rで素直な離散フーリエ変換 (discrete Fourier transform: DFT)とCooley-Tukey型高速フーリエ変換 (fast Fourier Transform: FFT)の処理時間を比較してみました．
Cooley-Tukey型FFT
FFTは，過去の記事を参照してください．
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
DFTとFFTのRスクリプト
離散フーリエ変換 (DFT)も，Cooley-Tukey型FFTも，私が自作したRスクリプトを使いました．Rのfftは使いません．当然，組み込み関数のfftは，すげー速いです．
以下のRスクリプトを実行すれば，離散フーリエ変換 (DFT)の関数DFT(x)と，Cooley-Tukey型FFTの関数FFT(x)が定義されます．
#######################################
# 離散フーリエ変換
DFT <- function(x){
# 時系列の長さ
N <- length(x)
# 回転因子の計算
W <- exp(-2i*(0:(N-1))*pi/N)
#######################
X <- c()
for(i in 1:N){
X[i] <- (x %*% W^(i-1))[1,1]
}
# 結果を返す
return(X)
}
#######################################
# FFT用のサブルーチン
FFT.sub <- function(x){
# 時系列の長さNを取得
N <- length(x)
# フーリエ変換結果の格納用
X <- c()
# 回転因子の計算 (1のN乗根)
k <- 0:(N/2-1)
W <- exp(-2i*k*pi/N)
################
if(N>2){
# 要素の半分をまとめて指定：1, 2, 3, ..., N/2
ii <- 1:(N/2)
# バタフライ演算の右上がり線 (/)に対応する計算
x1 <- x[ii]+x[ii+N/2] # 足すだけ
# バタフライ演算の右下がり線 (\)に対応する計算
x2 <- (x[ii]-x[ii+N/2])*W[ii] # 引いて回転因子をかける
# 長さが半分N/2のバタフライ演算に引き継ぐ (再帰的呼び出し)
X1 <- FFT.sub(x1)
X2 <- FFT.sub(x2)
# 最終的な結果をXに格納して出力
X[ii] <- X1
X[ii+N/2] <- X2
return(X)
}else{
# N=2の場合
X[1] <- x[1] + x[2] # 右上がり線 (/)に対応する計算（足すだけ）
X[2] <- x[1] - x[2] # 右下がり線 (\)に対応する計算（引くだけ）
return(X)
}
}
#######################################
# ビットリバースで順番を指定
# 10進数を2進数に変換してビット反転し，10進数に戻す
bit.rev <- function(n){
i.rev <- c()
i.tmp <- 0
for(i in 1:n){
pow2 <- 2^(i-1)
j <- 1:pow2
i.rev[2*j-1] <- i.tmp[j]
i.rev[2*j] <- i.tmp[j]+pow2
i.tmp <- i.rev
}
return(i.rev)
}
#######################################
# Cooley-TukeyのRadix-2アルゴリズム
FFT <- function(x){
# 時系列の長さ
N <- length(x)
# 2のn乗
n <- log(N)/log(2)
# バタフライ演算
X <- FFT.sub(x)
# ビットリバースで順番の並び替え
X[bit.rev(n)+1] <- X
# 結果を返す
return(X)
}
処理時間の比較
時系列の長さを変えて，処理時間を計測してみました．図にしたものが以下です．両対数目盛でのプロットになってます．
上の図では，両対数で直線に見えるので処理時間はのべき乗に比例しています．べき指数を推定したところ，DFTの処理時間は
に比例し，DFTの処理時間は
に比例していました．FFTの計算速度は
と言われますので，もっと，
がでかくなるとずれが出てくるのだと思います．
今回の数値実験に使ったRスクリプトはこれです．
# 処理時間保存用
T.FFT <- c()
T.DFT <- c()
# 時系列長
N <- 2^(6:10)
n <- length(N)
# 繰り返し回数
n.rep <- 100
# 処理時間計測
for(i in 1:n){
x <- rnorm(N[i])
T.FFT[i] <- system.time(for(j in 1:n.rep) X <- FFT(x))[[1]]/n.rep
T.DFT[i] <- system.time(for(j in 1:n.rep) X <- DFT(x))[[1]]/n.rep
}
# 最大最小
y.min <- min(T.DFT,T.FFT)
y.max <- max(T.DFT,T.FFT)
# plot
plot(N,T.DFT,log="xy",pch=16,col=2,ylim=c(y.min,y.max),xaxt="n",yaxt="n",xlab="length N",ylab="time [s]",cex.lab=1.5)
points(N,T.FFT,pch=2,col=4)
# 軸の描画
axis(1,at=N,cex.axis=1.2)
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.01)
a123<-paste(collapse=",",paste(sep="","expression(10^",-5:5,")"))
a123b<-paste(sep="","axis(2,las=1,at=10^(-5:5),label=c(",a123,"),cex.axis=1.2,tck=-0.02)")
eval(parse(text=a123b))
#
legend("bottomright", legend = c("DFT","FFT"), col = c(2,4), pch = c(16,2))
fit.DFT <- lm(log(T.DFT)~log(N))
fit.FFT <- lm(log(T.FFT)~log(N))
curve(exp(fit.DFT$coefficients[[1]]+fit.DFT$coefficients[[2]]*x),add=TRUE)
curve(exp(fit.FFT$coefficients[[1]]+fit.FFT$coefficients[[2]]*x),add=TRUE)
lines(N,exp(fit.DFT$coefficients[[1]])*N^fit.DFT$coefficients[[2]],col="pink",lty=2)
lines(N,exp(fit.FFT$coefficients[[1]])*N^fit.FFT$coefficients[[2]],col="skyblue",lty=2)