ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【Rで加速度解析】ATR-Promotionsの小型無線多機能センサ TSND151データの読み込み

今回はATR-Promotionsが販売している小型無線多機能センサTSND151
小型無線多機能センサ「TSND121/151」 |ATR-Promotions
のデータファイルを,Rで読み込むときのお話です.ポイントは,計測時刻をどうやってRの時刻変数に変更するかということです.

TSND151のデータ形式

 TSND151で計測されたデータをcsvファイルに書き出すと,

 センサ名-年月日-時分秒ミリ秒.csv

のルールでファイル名が付けられます.このファイルの中身は,

ags,51745030,-9269,-13,3381,108,12,-27
ags,51745040,-9269,45,3403,93,24,-3
ags,51745050,-9306,60,3456,105,3,9
ags,51745060,-9284,23,3361,81,-1,49
ags,51745070,-9291,-1,3312,56,-7,58
ags,51745080,-9321,28,3429,53,12,86
 ... 

のようになってます.各列は以下の変数を表しています.

データ種別 時刻 加速度 X 加速度 Y 加速度 Z 角速度 X 角速度 Y 角速度 Z
ags 51745030 -9269 -13 3381 108 12 -27
ags 51745040 -9269 45 3403 93 24 -3
... ... ... ... ... ... ... ...

 ここで,時刻は2列目ですが,この値を見ただけでは,何時か分かりません.以下では,Rを使って,この値を日時に変換する方法を説明します.

TSND151の計測日時

 TSND151の計測結果のcsvファイルにある「時刻」は,計測日の 00:00:00から経過時間を表しています.単位は,ミリ秒です.

 ファイルの中には「計測日」は書かれていないです.「計測日」はファイル名に書かれています.

 例えば,2022年10月14日14:30:00から計測を開始したのであれば,ファイル名は,

 subject-20221014-143000000.csv

となります.ですので,「20221014」の部分を読めば計測日が20221014
だということが分かります.

 ということで,ファイル名から計測日を読み取り,ファイル中の2列目の値 (ミリ秒単位になってます)を,計測日の00:00:00から足してやれば,計測日時が特定できます.

Rスクリプトの例

 RでTSND151の計測日時を知りたければ,以下のスクリプトを実行してください.timeという列が追加され,それが計測日時になってます.

# TSND151データのファイル名
FN.ATR <- "subject-20221014-143000000.csv"
# ファイル名を分割し,日付を抽出
tmp <- strsplit(FN.ATR,"-")
tmp.Date <- paste(substr(tmp[[1]][2],1,4),"-",substr(tmp[[1]][2],5,6),"-",substr(tmp[[1]][2],7,8),sep="")
# 計測日の00:00:00を設定
TIME0 <- as.POSIXct(tmp.Date,tz="Japan")
# TSND151データの読み込み
TMP <- read.csv(FN.ATR,header=FALSE)
# 日時をtimeとして追加
TMP$time <- TIME0+TMP$V2/1000

【Rで高速フーリエ変換】素直な離散フーリエ変換とFFTの処理時間の比較

今回は,Rで素直な離散フーリエ変換 (discrete Fourier transform: DFT)とCooley-Tukey型高速フーリエ変換 (fast Fourier Transform: FFT)の処理時間を比較してみました.

素直な離散フーリエ変換

 時系列を\{x[0], x[1], x[2], x[3], \cdots, x[N-1]\}とすれば,離散フーリエ変換は,

 \displaystyle X(k) = \sum_{i = 0}^{N-1} x[i] \, e^{-{\bf i} 2 \pi k i /N}

です.

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)
}

処理時間の比較

 時系列の長さNを変えて,処理時間を計測してみました.図にしたものが以下です.両対数目盛でのプロットになってます.

DFTとFFTの処理時間の比較

 
 上の図では,両対数で直線に見えるので処理時間はNのべき乗に比例しています.べき指数を推定したところ,DFTの処理時間はN^2に比例し,DFTの処理時間はN^1に比例していました.FFTの計算速度はO(N \log N)と言われますので,もっと,Nがでかくなるとずれが出てくるのだと思います.

 今回の数値実験に使った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)

まとめ

 DFTとFFTの処理時間を比較した今回の結果では,FFTが妥当な速さを示しました.ということは,私が書いたRスクリプトFFTはきちんとかけていると思います.

 私は,毎日長時間作業をしますが,夜になると,今日も,だらだらしただけで,たいしたことは何もしなかったなーと感じてしまいます.今日も,夕方にセミナーに参加して,他に何もできなかったなーと思いましたが,そういえば,生物物理の講義をして,企業団体向けの講演もしたことを思い出しました.具体的に思い返すと意外と仕事をこなしていますが,できていない仕事の量が多すぎて,やってない感じがします.高速はいいなと思いました.

【Rで高速フーリエ変換】N=16でCooley-Tukey型FFTを考えてみる

FFTを理解してもらいたい,というか,自分がわかりやすく説明できるようになりたいので,しつこくFFTの話をやっていきます.今回は,時系列の長さN=16でやっていきます.

基本事項

回転因子W

 フーリエ変換の式の表現をシンプルにするために,回転因子W_Nを定義します.

 \displaystyle 
 W_N=e^{-{\bf i} 2 \pi /N}

オイラーの式で,

 \begin{eqnarray} e^{-{\bf i} 2 \pi /N} &=& \cos\left(- \frac{2 \pi}{N} \right) + {\bf i} \sin\left(- \frac{2 \pi}{N} \right) \\
 &=& \cos\left(\frac{2 \pi}{N}\right) - {\bf i} \sin\left( \frac{2 \pi}{N} \right)
\end{eqnarray}

となるので,

 \begin{eqnarray} W_N^{k} = e^{-{\bf i} 2 \pi k/N} &=&  \cos\left(\frac{2 \pi k}{N}\right) - {\bf i} \sin\left( \frac{2 \pi k}{N} \right)
\end{eqnarray}

です.ということで,

 \displaystyle 
 W_N^{0}= 1 \displaystyle 
 W_N^{N}= 1 \displaystyle 
 W_N^{N/2}= -1

となります.

離散フーリエ変換を周波数の偶数番と奇数番に分けると

 時系列を\{x[0], x[1], x[2], x[3], \cdots, x[N-1]\}とすれば,離散フーリエ変換は,

 \displaystyle X(k) = \sum_{i = 0}^{N-1} x[i] \, W_N^{ki}

です.時系列の長さが2の倍数のとき,前に説明したように,kについて,偶数番 (k = 0, 2, \cdots, N/2)と奇数番 (k = 1, 3, \cdots, N/2-1)に分けて以下の式が成り立ちます.

\kappa = 0, 1, 2, \cdots, N/2-1として,
偶数番 k=2\kappaでは,

 \begin{eqnarray}  X(2\kappa) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{\kappa i}
 \end{eqnarray}

奇数番 k=2\kappa+1では,

 \begin{eqnarray}
X(2\kappa+1) &=& \sum_{i = 0}^{N/2-1}\left\{ \left(x[i]-x[N/2+i] \right) W_{N}^{i} \right\} W_{N/2}^{\kappa i}
 \end{eqnarray}

この式の導出については,詳しく説明してあるページがいくつもあるので,検索して調べてください.

N=16でFFTの計算式を確認

 N=16=2^4で,時系列を\{x[0], x[1], x[2], x[3], \cdots, x[15]\}として,FFTの計算過程を見ていきます.

1回目の分割

 周波数の偶数番 k=0, 2, \cdots, 14では,

 \begin{eqnarray}  X(0) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{0 i} \\
 X(2) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{i} \\
 X(4) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{2 i} \\
 X(6) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{3 i} \\
 X(8) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{4 i} \\
 X(10) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{5 i} \\
 X(12) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{6 i} \\
 X(14) &=& \sum_{i = 0}^{7} \left(x[i]+x[8+i] \right) W_{8}^{7 i} \\
 \end{eqnarray}

です.ここで,\left(x[i]+x[8+i] \right)の計算が共通なので,

 x_8^{(0)}[i] = x[i]+x[8+i]

として,i=0, 1, \cdots, 7について足し算すれば,全体で足し算する回数が減ります.xに添え字がたくさんつきますが,計算の確認のため記号を付けておきます.上付き添え字の{}^{(0)}は偶数番周波数の置き換え (足し算),上付き添え字の{}^{(1)}は奇数番周波数の置き換え (引き算&回転因子)を表すことにします.

 周波数の奇数番 k=1, 3, \cdots, 15では,

 \begin{eqnarray}  X(1) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{0 i} \\
 X(3) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{1 i} \\
 X(5) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{2 i} \\
&\vdots&\\
 X(15) &=& \sum_{i = 0}^{7}\left\{ \left(x[i]-x[8+i] \right) W_{16}^{i} \right\} W_{8}^{7 i} \\
 \end{eqnarray}

です.ここでは,\left(x[i]-x[8+i] \right) W_{16}^{i}の計算が共通なので,

 x_8^{(1)}[i] = \left(x[i]-x[8+i] \right) W_{16}^{i}

とします.

 上半分に偶数番目の周波数,下半分に奇数番目の周波数をならべて書くと,

 \begin{eqnarray}  X(0) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\, W_{8}^{0 i} \\
 X(2) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\,  W_{8}^{1 i} \\
 X(4) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{2 i} \\
 X(6) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{3 i} \\
&\vdots& \\
 X(14) &=& \sum_{i = 0}^{7} x_8^{(0)}[i] \, W_{8}^{7 i} \\
X(1) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\, W_{8}^{0 i} \\
 X(3) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\,  W_{8}^{1i} \\
 X(5) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{2 i} \\
 X(7) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{3 i} \\
&\vdots& \\
 X(15) &=& \sum_{i = 0}^{7} x_8^{(1)}[i] \, W_{8}^{7 i} \\
 \end{eqnarray}

となります.ここで,右辺について気づくことは,上半分と下半分のそれぞれで,N=8のときのフーリエ変換になっていると言うことです.このフーリエ変換

 \begin{eqnarray} X_8^{(0)} (k) &=& \sum_{i = 0}^{7} x_8^{(0)}[i]\, W_{8}^{k i}\\
X_8^{(1)} (k) &=& \sum_{i = 0}^{7} x_8^{(1)}[i]\, W_{8}^{k i} \end{eqnarray}

と表し,さらに,偶数番,奇数番で分割すると,

 \begin{eqnarray}  X(0) &=&X_8^{(0)}(0) = \sum_{i = 0}^{3} \left(x_8^{(0)}[i]+x_8^{(0)}[4+i] \right)\, W_{4}^{0 i} \\
 X(2) &=& X_8^{(0)}(1) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{0 i}\\
 X(4) &=& X_8^{(0)}(2) =  \sum_{i = 0}^{3} \left(x_8^{(0)}[i]+x_8^{(0)}[4+i] \right)\, W_{4}^{1 i} \\
 X(6) &=& X_8^{(0)}(3) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{1 i}\\
&\vdots&\\
 X(14) &=& X_8^{(0)}(7) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{3 i}\\
X(1) &=&X_8^{(1)}(0) = \sum_{i = 0}^{3} \left(x_8^{(1)}[i]+x_8^{(1)}[4+i] \right)\, W_{4}^{0 i} \\
 X(3) &=& X_8^{(1)}(1) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{0 i}\\
 X(5) &=& X_8^{(1)}(2) =  \sum_{i = 0}^{3} \left(x_8^{(1)}[i]+x_8^{(1)}[4+i] \right)\, W_{4}^{1 i} \\
 X(7) &=& X_8^{(1)}(3) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{1 i}\\
&\vdots&\\
 X(15) &=& X_8^{(1)}(7) = \sum_{i = 0}^{3}\left\{ \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i} \right\} W_{4}^{3 i}
 \end{eqnarray}

となります.

2回目の分割

 上の式にも共通の部分があることに気づくので,

 \begin{eqnarray}  x_4^{(00)}[i] &=& x_8^{(0)}[i]+x_8^{(0)}[4+i] \\
x_4^{(01)}[i] &=& \left(x_8^{(0)}[i]-x_8^{(0)}[4+i] \right) W_{8}^{i}\\
x_4^{(10)}[i] &=& x_8^{(1)}[i]+x_8^{(1)}[4+i] \\
x_4^{(11)}[i] &=& \left(x_8^{(1)}[i]-x_8^{(1)}[4+i] \right) W_{8}^{i}
 \end{eqnarray}

と置きます.N=4フーリエ変換を,

 \begin{eqnarray}  X_4^{(00)} (k) &=& \sum_{i = 0}^{3} x_4^{(00)}[i]\, W_{4}^{k i}\\
X_4^{(01)} (k) &=& \sum_{i = 0}^{3} x_4^{(01)} [i]\, W_{4}^{k i}\\
X_4^{(10)} (k) &=& \sum_{i = 0}^{3} x_4^{(10)}[i]\, W_{4}^{k i}\\
X_4^{(11)} (k) &=& \sum_{i = 0}^{3} x_4^{(11)} [i]\, W_{4}^{k i}
 \end{eqnarray}

と表し,またまた,偶数番,奇数番で分割すると,

 \begin{eqnarray}  X(0) &=&X_8^{(0)}(0) = X_4^{(00)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\, W_{2}^{0 i} \\
 X(2) &=& X_8^{(0)}(1) = X_4^{(01)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(01)}[i]+x_4^{(01)}[2+i] \right) W_{2}^{0 i}\\
 X(4) &=& X_8^{(0)}(2) = X_4^{(00)}(1)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(00)}[i]-x_4^{(00)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(6) &=& X_8^{(0)}(3) = X_4^{(01)}(1)  = \sum_{i = 0}^{1} \left\{ \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(8) &=& X_8^{(0)}(4) = X_4^{(00)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\, W_{2}^{1 i} \\
 X(10) &=&X_8^{(0)}(5) = X_4^{(01)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(01)}[i]+x_4^{(01)}[2+i] \right)\, W_{2}^{1 i} \\
 X(12) &=& X_8^{(0)}(6) = X_4^{(00)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(00)}[i]-x_4^{(00)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 X(14) &=& X_8^{(0)}(7) = X_4^{(01)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
X(1) &=&X_8^{(1)}(0) = X_4^{(10)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\, W_{2}^{0 i} \\
 X(3) &=& X_8^{(1)}(1) = X_4^{(11)}(0)  = \sum_{i = 0}^{1} \left(x_4^{(11)}[i]+x_4^{(11)}[2+i] \right) W_{2}^{0 i}\\
 X(5) &=& X_8^{(1)}(2) = X_4^{(10)}(1)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(10)}[i]-x_4^{(10)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(7) &=& X_8^{(1)}(3) = X_4^{(11)}(1)  = \sum_{i = 0}^{1} \left\{ \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{0 i}\\
 X(9) &=& X_8^{(1)}(4) = X_4^{(10)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\, W_{2}^{1 i} \\
 X(11) &=&X_8^{(1)}(5) = X_4^{(11)}(2)  = \sum_{i = 0}^{1} \left(x_4^{(11)}[i]+x_4^{(11)}[2+i] \right)\, W_{2}^{1 i} \\
 X(13) &=& X_8^{(1)}(6) = X_4^{(10)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(10)}[i]-x_4^{(10)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 X(15) &=& X_8^{(1)}(7) = X_4^{(11)}(3)  = \sum_{i = 0}^{1}\left\{ \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_{4}^{i} \right\} W_{2}^{1 i}\\
 \end{eqnarray}

となります.

3回目と4回目の分割

 ちょっと込み入ってきたので,自信がなくなってきましたが,ここでも,共通部分を,

 \begin{eqnarray}  x_2^{(000)}[i] &=& \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right)\\
x_2^{(001)}[i] &=& \left(x_4^{(00)}[i]+x_4^{(00)}[2+i] \right) W_4^i\\
x_2^{(010)}[i] &=& \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) \\
x_2^{(011)}[i] &=& \left(x_4^{(01)}[i]-x_4^{(01)}[2+i] \right) W_4^i\\
x_2^{(100)}[i] &=& \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right)\\
x_2^{(101)}[i] &=& \left(x_4^{(10)}[i]+x_4^{(10)}[2+i] \right) W_4^i\\
x_2^{(110)}[i] &=& \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) \\
x_2^{(111)}[i] &=& \left(x_4^{(11)}[i]-x_4^{(11)}[2+i] \right) W_4^i
 \end{eqnarray}

と置きます.さらにさらに,N=2フーリエ変換

 \begin{eqnarray}  X_2^{(000)} (k) &=& \sum_{i = 0}^{1} x_2^{(000)}[i]\, W_{2}^{k i}\\
X_2^{(001)} (k) &=& \sum_{i = 0}^{1} x_2^{(001)} [i]\, W_{2}^{k i} \\
X_2^{(010)} (k) &=& \sum_{i = 0}^{1} x_2^{(010)}[i]\, W_{2}^{k i}\\
&\vdots&\\
X_2^{(111)} (k) &=& \sum_{i = 0}^{1} x_2^{(111)} [i]\, W_{2}^{k i}
 \end{eqnarray}

と表せば,

 \begin{eqnarray}  X(0) &=& X_8^{(0)}(0) = X_4^{(00)}(0) = X_2^{(000)} (0)  \\
&=&  x_2^{(000)}[0]\, W_{2}^{0} + x_2^{(000)}[1]\, W_{2}^{0} = x_2^{(000)}[0]+ x_2^{(000)}[1]\\
 X(2) &=& X_8^{(0)}(1) = X_4^{(01)}(0)  =  X_2^{(010)} (0)\\
&=&  x_2^{(010)}[0]\, W_{2}^{0} + x_2^{(010)}[1]\, W_{2}^{0} = x_2^{(010)}[0]+ x_2^{(010)}[1]\\
&\vdots& \\
 \end{eqnarray}

のようになります.この変形をして,さらに,これまでと同様のルールでxに添え字をつけて記述すれば,

 \begin{eqnarray}  X(0) &=& x_2^{(000)}[0]+ x_2^{(000)}[1] = x_1^{(0000)}[0]\\
X(1) &=& x_2^{(100)}[0]+ x_2^{(100)}[1] = x_1^{(1000)}[0]\\
X(2) &=& x_2^{(010)}[0]+ x_2^{(010)}[1] = x_1^{(0100)}[0]\\
 X(3) &=& x_2^{(110)}[0]+ x_2^{(110)}[1] = x_1^{(1100)}[0]\\
X(4) &=& x_2^{(001)}[0]+ x_2^{(001)}[1] = x_1^{(0010)}[0]\\
X(5) &=& x_2^{(101)}[0]+ x_2^{(101)}[1] = x_1^{(1010)}[0]\\
X(6) &=& x_2^{(011)}[0]+ x_2^{(011)}[1] = x_1^{(0110)}[0]\\
X(7) &=& x_2^{(111)}[0]+ x_2^{(111)}[1] = x_1^{(1110)}[0]\\
X(8) &=& x_2^{(000)}[0]- x_2^{(000)}[1] = x_1^{(0001)}[0]\\
X(9) &=& x_2^{(100)}[0]- x_2^{(100)}[1]  = x_1^{(1001)}[0] \\
X(10) &=& x_2^{(010)}[0]- x_2^{(010)}[1] = x_1^{(0101)}[0]  \\
X(11) &=& x_2^{(110)}[0]- x_2^{(110)}[1]   = x_1^{(1101)}[0]\\
X(12) &=& x_2^{(001)}[0]- x_2^{(001)}[1] = x_1^{(0011)}[0]  \\
X(13) &=& x_2^{(101)}[0]- x_2^{(101)}[1]  = x_1^{(1011)}[0] \\
X(14) &=& x_2^{(011)}[0]- x_2^{(011)}[1] = x_1^{(0111)}[0]  \\
X(15) &=& x_2^{(111)}[0]- x_2^{(111)}[1] = x_1^{(1111)}[0]
 \end{eqnarray}

となりました.最後は,順番に並べて書きました.

ここまでの計算で何がやりたかったのか

 Cooley-Tukey型FFTで何をやっているかと言えば,周波数の偶数番と奇数番を分割を繰り返しているだけです.この分割というのは,それぞれ,

  • 足し算を計算する x[i]+x[N/2+i]
  • 引き算して回転因子をかける  \left(x[i]-x[N/2+i] \right) W_{N}^{i}

に対応します.そして,最後まで,この2パターンの計算を繰り返しました.

 上の変形では最終的に,

 \begin{eqnarray}  X(0) &=& x_1^{(0000)}[0]\\
X(1) &=& x_1^{(1000)}[0]\\
 X(2) &=& x_1^{(0100)}[0]\\
 X(3) &=& x_1^{(1100)}[0]\\
X(4) &=& x_1^{(0010)}[0]\\
&\vdots&\\
X(14) &=& x_1^{(0111)}[0]  \\
X(15) &=& x_1^{(1111)}[0]
 \end{eqnarray}

という表現がえられました.

 xの上付き添え字の部分(0000), (0100), \cdotsは,周波数の偶数番と奇数番の分割で,どちら側だったかという順番を表しています.

 例えば,X(0)に登場する(0000)は,1回目偶数(0),2回目偶数(0),3回目偶数(0),4回目偶数(0)だったということです.

 X(5)に登場する(1010)は,1回目奇数(1),2回目偶数(0),3回目奇数(1),4回目偶数(0)だったということです.

 そして,X(5)に登場する(1010)は,5を2進数にする計算と同じことをしています.ただし,0, 1の値の並べ方が右から左ではなく,左から右になります.これが,なぜ,ビット反転をするとうまくいくのかの理由です.

 今回のお話は,前回の記事で紹介したバタフライダイアグラムを,式で表してみたということです.下の図のようなバタフライダイアグラムの構造と今回の結果の関係が見えてくるでしょうか.図の中に書いた0, 1と,x^{(\cdots)}の上付き添え字の部分(\cdots)との対応を考えてみてください.

N=16のバタフライダイアグラム.マイナスの記号や回転因子は省略してあります.その代わり,周波数の偶数番と奇数番に対応する演算の部分に,それぞれ,0と1を書きました.

まとめ

 Cooley-Tukey型FFTについては,いろいろな説明の仕方があります.例えば,行列を使った説明があります.私は,行列を使った説明はわかりやすいと思います.

しかし,私は個性的な説明を試みてみたかったので,今回は,周波数の分割をすべて書いてみるというアプローチをしました.数式が嫌いな人が多いので,わかりにくくなったかもしれません.私としては,数式からイメージがわいてくるので,知識の整理には役立ちました.

【Rで高速フーリエ変換】バタフライダイアグラムを読む

私は肉まんの白い部分が好きです.ほのかに甘い,あの白い部分だけ食べたいです.

 ということで,今回は最近のFFTのお話
chaos-kiyono.hatenablog.com
の続きです.

 FFTの説明で,バタフライ演算ってのが登場します.バタフライ演算は下の図のようなバタフライダイアグラムで表すことができます.

バタフライ演算の模式図 (N=4のとき).

「バタフライ」ってのは,FFTの計算を図で表すと蝶蝶の羽のように見えることに由来していると思います (知らんけど).今回はバタフライダイアグラムが表す計算過程を説明します.

右上がりは足し,右下がりは引く

 時系列の長さNが4の時系列\{x[0], x[1], x[2], x[3]\}を考えます.前回はRでの実装を想定して,時系列の順番について1から番号をふりましたが,今回は0からです.ビット反転を説明するために,一般的な表現を使います.

 N=4のとき,バタフライダイアグラムに補足の記号 (プラスマイナスと回転因子)を加えると下図のようになります.

バタフライダイアグラム (N=4のとき).

 上の図では,左から右へ計算が進んでいきます.赤 (右上がり線)と青 (右下がり線)で描いた線は,それぞれ,「足し算」と「引き算」を表します.

 前回,長さNが2の倍数であれば,時系列\{x[0], x[1], \cdots, x[N-1]\}の離散フーリエ変換

 \displaystyle X(k) = \sum_{i = 0}^{N-1} x[i] \, W_N^{ki}

は,偶数番目の周波数の離散フーリエ変換

 \begin{eqnarray}  X(2k) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{ki}
 \end{eqnarray}

と,奇数番目の周波数の離散フーリエ変換

 \begin{eqnarray}  X(2k+1) &=& \sum_{i = 0}^{N/2-1}\left\{ \left(x[i]-x[N/2+i] \right) W_{N}^{i} \right\} W_{N/2}^{k i}
 \end{eqnarray}

に分けられることを説明しました (詳しくは説明してないけど).ここで,W_{N}=e^{-{\bf i} 2 \pi /N}です.

 バタフライダイアグラムの上半分は偶数番目の周波数の離散フーリエ変換,下半分は奇数番目の周波数の離散フーリエ変換に対応します.

 ダイアグラムの見方を理解するために,以下ではバタフライダイアグラムのパーツに注目してみます.

右上がりの線 (足し算)

 下の図のように,右上がりの線 (赤)は,でつないだ2本の線に対応する値を足すことを表しています.

2つの平行線をつなぐ右上がりの線は足し算

 これは,偶数番目の周波数の離散フーリエ変換

 \begin{eqnarray}  X(2k) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{ki}
 \end{eqnarray}

に登場する

  \left(x[i]+x[N/2+i] \right)

の計算を表しています.

右下がりの線 (引いて回転演算子をかける)

 下の図のように右下がりの線 (青)は,でつないだ2本の線に対応する値を引くことを表しています.さらに,下のの近くにW_4^1と書いてあるので,それをかけます.

右下がりの線は引き算.回転演算子Wが書いてあれば,それをかける.

 これは,奇数番目の周波数の離散フーリエ変換

 \begin{eqnarray}  X(2k+1) &=& \sum_{i = 0}^{N/2-1}\left\{ \left(x[i]-x[N/2+i] \right) W_{N}^{i} \right\} W_{N/2}^{k i}
 \end{eqnarray}

に登場する

  \left(x[i]-x[N/2+i] \right) W_{N}^{i}

の計算を表しています.

回転演算子

 回転演算W_Nをかける必要があるのは,奇数番目の周波数の離散フーリエ変換です.つまり,バタフライダイアグラムでは下半分の部分です.かける必要がある回転演算子 W_{N}^{i}は, i=0 ,1 , \cdots, N/2-1の場合です.したがって,N=4のときは,上から下に,回転演算子を2つ (N/2=2),つまり,

 W_{4}^{0}, W_{4}^{1}

をかけます.図に書き込めば,下のようになります.

バタフライダイアグラムの回転演算子

ビット反転でフーリエ変換の結果を並び替え

 以上の説明で,バタフライダイアグラムが表している計算アルゴリズムを理解できたでしょうか.例えば,下の図のような,N=8のバタフライダイアグラムが表す計算式をイメージできますか.

N=8のバタフライダイアグラム

 とりあえず,バタフライダイアグラムの読み方は分かったとしても,右端に書いた離散フーリエ変換の結果X(k)の順番がめちゃくちゃです.この順番は,どうやって決められるのでしょうか.そこで,登場するのがビット反転 (ビットリバース)です.

ビット反転の計算

 ビット反転では,10進数を2進数に変換して,2進数の並びを反転してから10進数に戻します.ただし,時系列の長さがNのときは,\log_2 N桁の2進数で表し,上位の桁がないときは,0を詰めます (わかった?).

 例えば,N=8のとき,\log_2 8 = 3桁の2進数を考えます.N=8のとき,1をビット反転してみます.10進数の1_{(10)}は,2進数で表しても1_{(2)}ですが,3桁にしたいので,0を詰めて,

 001_{(2)}

とします.数字の順番を逆にならべると,

 100_{(2)}

になります.これを10進数にすると,

 100_{(2)} = 1\times 2^2 + 0 \times 2^1 + 0 \times 2^0 = 4

になります.確かに,上の図では,x[1]の右端に,X(4)が書いてあります.

なぜ,ビット反転でうまくいくの?

 ビット反転で対応する周波数の順番が決まる理由を考えてみます.N=8の例を表した下の図を見てください.この図を使って考えてみます.

なぜ,ビット反転 (リバース)するとうまくいくの

 まず,バタフライダイアグラムの最初の段 (一番左)で,上半分にあれば,それは偶数番目の周波数です.だって,

 \begin{eqnarray}  X(2k) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{ki}
 \end{eqnarray}

の計算をしているからです.N=8の例では,周波数の番号は 0, 2, 4, 6のどれかです.逆に,下半分は奇数で, 1, 3, 5, 7のどれかです.

 偶数ってことは,2進数で表したとき1の位は0です.奇数ってことは,2進数で表したとき,1の位は1です.この桁は,左端のx[i]の順番iを2進数で表したときの,100の位に一致します.

 後は計算式を追ってやればビット反転との対応が見えてきます.

 例えば,X(4)では,4は偶数なので,最初は偶数のフーリエ変換の式

 \begin{eqnarray}  X(2k) &=& \sum_{i = 0}^{N/2-1} \left(x[i]+x[N/2+i] \right) W_{N/2}^{ki}
 \end{eqnarray}

を使います.これで,2進数で表したときの1桁目は,0に確定です.この式では,周波数が2kフーリエ変換が,半分のkでのフーリエ変換に書き換わっています.

 ということは,X(4)4を2で割って,X(2)フーリエ変換の式を次の段階で使うことになります.2は偶数なので,バタフライダイアグラムの真ん中の段は,上半分のパスを通ります.なので,2進数で表したときの10の位は,0になります.

 最後は,X(2)2を2で割って,X(1)フーリエ変換の式を使うことになります.1は奇数なので,バタフライダイアグラムの一番右側の段は,下半分のパスを通ります.なので,2進数で表したときの100の位は,1になります.

 順番を整理すると,0→0→1です.これは,4を2進数に変換した100_{2}を逆に並べたものと同じです.

 エビフライダイアグラムの左側のx[i]の順番iは,上の図に書いた0と1を左から右へ並べた2進数になってます.一方,右側のX(k)の順番k (フーリエ変換結果)では,左が1の位に対応するので,上の図に書いた0と1を右から左へ並べた2進数になってます.

 こんな説明で,理解できますでしょうか.

まとめ

 Rのコマンドとか,パッケージを使って,それっぽい結果出すだけでなく,使う道具にも興味や疑問をもって基礎知識を身につけてください.

【Rで高速フーリエ変換】2のべき乗の長さの時系列用FFTの実装についての解説

前回,2のべき乗の長さの時系列を,高速に離散フーリエ変換 (FFT)できるアルゴリズムのRスクリプトを作りました (高速じゃないけど).
chaos-kiyono.hatenablog.com

 このFFTアルゴリズムは,Cooley-Tukey型と呼ばれるそうです (知らんけど).今回はこのアルゴリズムについて少し解説します.

要は周波数について奇数番と偶数番に分けるだけ

 以下では,長さが2のべき乗の時系列\{x[1], x[2], \cdots, x[N]\}を離散フーリエ変換することを考えます.時系列解析の本では,データの番号を0からN-1で考えることが多いですが,Rでは,配列の番号が1からはじまるので,データの番号は1からNにしました.他の教科書とかインターネットの解説と番号がずれているので注意してください.

 このとき,離散フーリエ変換 X(k) は,

 \displaystyle X(k) = \sum_{i = 1}^{N} x[i]\, e^{-{\bf i} 2 \pi k (i-1)/N }

です.kは,0からN-1までです.サンプリング間隔が1のときは周波数は,f_k = k/Nです.

 Cooley-Tukey型FFTって,要は,kについて,奇数番と偶数番を分けるだけです.

N=4の例

 簡単な例として,N=4でやってみます (速くないけど).kが偶数で,k=2のとき (eの肩の数値に注目),

 \begin{eqnarray} X(k=2) &=& \sum_{i = 1}^{4} x[i] \, e^{-{\bf i} 2 \pi \cdot 2 (i-1)/N }\\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 1/4} + x[3]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 2/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 2 \cdot 3/4} \\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/2} + x[3]\, e^{-{\bf i} 2 \pi \cdot 4/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 6/4} \\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/2} + x[3]+ x[4]\, e^{-{\bf i} 2 \pi \cdot 1/2} \\
 &=& \left(x[1] + x[3] \right) +  \left( x[2]+ x[4]\right) e^{-{\bf i} 2 \pi \cdot 1/2} \\
 &=& \sum_{i = 1}^{2} \left(x[i]+x[i+2] \right) \, e^{-{\bf i} 2 \pi \cdot (i-1)/2}
\end{eqnarray}

と変形できます.この変形の注目点は

  • 最初は\displaystyle \sum_{i = 1}^{4}で4項の和だった部分が,変形すると\displaystyle \sum_{i = 1}^{2}になり,半分の2項の和になった.
  • e^{-{\bf i} 2 \pi k (i-1)/N}の部分のk=2が半分の1になった (周波数が半分になった).
  • e^{-{\bf i} 2 \pi k (i-1)/N}の部分のN=4が半分の2になった.

ということです.

 また,kが奇数で,k=1のとき,

 \begin{eqnarray} X (k=1) &=& \sum_{i = 1}^{4} x[i] \, e^{-{\bf i} 2 \pi \cdot 1 \cdot (i-1)/N }\\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1 \cdot 1/4} + x[3]\, e^{-{\bf i} 2 \pi \cdot 1  \cdot 2/4}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 1 \cdot 3/4} \\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/4} + x[3]\, e^{-{\bf i} \pi}+ x[4]\, e^{-{\bf i} 2 \pi \cdot 3/4} \\
 &=& x[1] + x[2]\, e^{-{\bf i} 2 \pi \cdot 1/4} - x[3] - x[4]\, e^{-{\bf i} 2 \pi \cdot 1/4} \\
 &=& \left(x[1] - x[3] \right) +  \left( x[2]- x[4]\right) e^{-{\bf i} 2 \pi \cdot 1/4} \\
 &=& \sum_{i = 1}^{2} \left\{\left(x[i]-x[i+2] \right) \, e^{-{\bf i} 2 \pi \cdot (i-1)/4}\right\} e^{-{\bf i} 2 \pi \cdot 0 \cdot (i-1)/2}
\end{eqnarray}

と変形できます (最後の式は結果を知っているので無理矢理書きました).やはり,\displaystyle \sum_{i = 1}^{4}で4項の和だった部分を,変形して\displaystyle \sum_{i = 1}^{2}にできました.

 一般にどんな関係式が成り立つか考えてみてください.

一般形

 Nが2の倍数のとき,離散フーリエ変換

 \displaystyle X(k) = \sum_{i = 1}^{N} x[i]\, e^{-{\bf i} 2 \pi k (i-1)/N }

について,以下の関係が成り立ちます.

 \begin{eqnarray}  X(2k) &=& \sum_{i = 1}^{N/2} \left(x[i]+x[N/2+i] \right) e^{-{\bf i} 2 \pi k \cdot (i-1)/(N/2)} \\
X[2k+1] &=& \sum_{i = 1}^{N/2}\left\{ \left(x[i]-x[N/2+i] \right) e^{-{\bf i} 2 \pi (i-1)/N} \right\} e^{-{\bf i} 2 \pi k \cdot (i-1)/(N/2)}
 \end{eqnarray}

すっきり書くために,W_{N}=e^{-{\bf i} 2 \pi /N}とすれば,

 \begin{eqnarray}  X(2k) &=& \sum_{i = 1}^{N/2} \left(x[i]+x[N/2+i] \right) W_{N/2}^{k(i-1)} \\
X (2k+1) &=& \sum_{i = 1}^{N/2}\left\{ \left(x[i]-x[N/2+i] \right) W_{N}^{i-1} \right\} W_{N/2}^{k(i-1)}
 \end{eqnarray}

となります.

 Nが2のべき乗の長さのとき,この分割を繰り返し行えるので,

 x[i]+x[N/2+i]

と,

 \left(x[i]-x[N/2+i] \right) W_{N}^{i-1}

を繰り返し計算すれば離散フーリエ変換できちゃいます.これが,Cooley-Tukey型FFTです.

Rスクリプトとの対応

 前回,RでCooley-Tukey型FFTスクリプトを書いてみました.
chaos-kiyono.hatenablog.com

 このスクリプトの中の

  # 要素の半分をまとめて指定: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] # 引いて回転因子をかける

の部分が,

 x[i]+x[N/2+i]

と,

 \left(x[i]-x[N/2+i] \right) W_{N}^{i-1}

の計算です.

次回予告:バタフライダイアグラムの読み方

 FFTのポイントはこれだけですが,次回,バタフライ演算の図の読み方と,ビット反転について説明します.

 未来から帰ってきました.つづきの記事は,これです.
chaos-kiyono.hatenablog.com

【Rで高速フーリエ変換】Rのfftは任意の長さで大丈夫だけど,学習用に基数2のFFTを自作

Rとか,Pythonとかを使って,誰でも簡単に高速フーリエ変換 (fast Fourier transform, FFT)を使える時代になりました.最近コロナで,人と会わなくなって,久しぶりに遠隔会議であった方々が老けたなーと感じることが多くなりました.同様に私も年寄りになって,時代の変化を懐かしむ気持ちが強くなっているのだと思います.

 ということで,今回はRでFFTのお話です.北川先生の時系列解析の本に,2のべき乗の長さの時系列を仮定したFFTの説明があります.今の時代に「2のべき乗の長さ」の話かと思いました.なぜなら,Rのfftコマンドは,任意の長さの時系列のフーリエ変換をきちんとやってくれます.つまり,2のべき乗の長さとか,気にする必要はありません.

 そうは言っても,理屈を理解するには2のべき乗の長さのFFT (基数2のFFT)が基本だと思いますので,学習用に基数2のFFTをRで実装してみました.残念ながら私の実力では,高速化はイマイチですので,実用向けではないことに注意してください (高速化の方法が分かれば教えてください).

基数2のFFT

 特に処理速度が速いわけではありませんが,基数2のFFTをRで自作しました.これは,学習用ですので,普段の時系列解析では,Rのfftコマンドを使ってください.時系列の長さが2のべき乗のときは,私のスクリプトでも,fftでも結果は同じです.

 私が作ったRスクリプトは以下です.これは,FFT.rdx2という自作関数の定義です.FFT.rdx2を使う前に,1回だけ,このスクリプトを実行してください.

#######################################
# FFT用のサブルーチン
FFT.sub <- function(x){
 # 時系列の長さNを取得
 N <- length(x)
 # フーリエ変換結果の格納用
 X <- c()
 # 回転因子の計算 (1のN乗根)
 k <- 0:(N/2-1)
 W <- complex(real=cos(2*k*pi/N),imaginary=sin(-2*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.rdx2 <- function(x){
 # 時系列の長さ
 N <- length(x)
 # 2のn乗
 n <- log(N)/log(2)
 # バタフライ演算
 X <- FFT.sub(x)
 # ビットリバースで順番の並び替え
 X[bit.rev(n)+1] <- X
 # 結果を返す 
 return(X)
}


 使い方は,時系列を小文字のxとして,FFT.rdx2(x)を実行するだけです.例えば,

# 時系列の長さ (必ず2のべき乗にしろ)
N <- 2^10
# 正規乱数を生成
x <- rnorm(N)
# 基数2のFFT
X <- FFT.rdx2(x)

を実行してみてください.大文字のXに結果が入ります.


 処理速度はfftと比べれば,かなり遅いです.

N <- 2^19
x <- rnorm(N)
###############
> system.time(X <- FFT.rdx2(x))
   ユーザ   システム       経過  
      2.37       0.00       2.37 
> system.time(X <- fft(x))
   ユーザ   システム       経過  
      0.03       0.00       0.03 

FFT.rdx2(x)では2.37秒,fftでは0.03秒かかったので,fftの方がかなり速いです.

 本日の課題は,「基数2のFFTを使って,2のべき乗の長さではない時系列に,0を詰めて2のべき乗の長さにしたときの影響を調べろ」にしておきます.

Rのfftはいつも速いわけではない

 Rのfftでは,基数2のFFTを一般の基数に拡張したアルゴリズムを使っています.時系列の長さが,たくさんの素数に分解できるときほど速いです.つまり,処理速度は単純に時系列の長さに依存するのではなく,どんだけ細かく素因数分解できるかに依存します.

 実際に,処理速度を測ってみます.時系列は正規乱数を生成して使います.

 私の環境では,時系列の長さがN=2^{18}=262144のとき,

> N <- 2^18
> x <- rnorm(N)
> system.time(X <- fft(x))
   ユーザ   システム       経過  
      0.01       0.00       0.01

0.01 秒でした.

 遅くなるのは,素因数分解できないときです.時系列の長さを素数N=99991にしてみると,

> N <- 99991
> x <- rnorm(N)
> system.time(X <- fft(x))
   ユーザ   システム       経過  
      8.25       0.00       8.25 

8.25秒になりました.時系列の長さは,前の例N=2^{18}=262144よりも短いですが,計算時間は800倍くらい長くかかりました.

 Nが素数のときは,Rが使っているfftアルゴリズムは,FFTとは呼べなくなります.Nが素数の例は極端ですが,2のべき乗じゃなくても,ある程度,素因数分解できれば計算時間はかなり短くなります.

 例えば,時系列の長さを素数99991から1つ除いて,時系列の長さをN=99990 = 2\times 3^2 \times 5 \times 11にすれば,

> N <- 99990
> x <- rnorm(N)
> system.time(X <- fft(x))
   ユーザ   システム       経過  
      0.01       0.00       0.01 

0.01秒になりました.2のべき乗じゃなくても,むちゃくちゃ速くなります.

そもそもFFTすると何が出てくるの?

 FFTというか,離散フーリエ変換の結果の見方が分からない人は,過去の記事を参照してください.
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com

まとめ

 私は,スマートウォッチ用の組み込みアルゴリズムの開発を,企業から依頼されることが結構あります.処理速度とバッテリーの持ちを考えれば,今でも基数2のFFTは,フーリエ変換の選択肢になるので,その理屈と特徴を理解しておくことが大切だと思います.

【Rの処理時間】計算時間短縮のための工夫 ― meanよりもsumが速い ―

プログラムを作って何か処理するとき,処理時間が短い方が良いに決まってます.ということで,今回はRの処理時間についてのお話です.

 処理時間を無駄に長くしないための一般論として,私が説明できることは,「forループは使わずにベクトルのまま計算した方が良い」,「ifをできるだけ避ける」くらいです.私は,Rの内部の仕組みについて知らないので,Rに特有の挙動には精通していません.

 今回は「標本平均の計算」を例として,実験的に処理時間を検証してみます.

サンプルデータの生成

 以下では,xにデータが入っているとします.xの具体例として,長さN,平均3,標準偏差1の正規乱数を用意するRスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8 # 10^8と同じ
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3,1)

3通りの方法で標本平均を計算

 比較する標本平均の計算方法は以下の3通りです.
1. コマンドmeanを使う

mean(x)

2. コマンドsumで和を計算して,xの長さで割る

sum(x)/length(x)

3. forループを使って,プログラムっぽく計算

 s <- 0
 n <- length(x)
 for(i in 1:n){
  s <- s+x[i]
 }
 s/n

計算時間を比較

 標本平均を計算する3通りの方法の処理時間を比べてみます.

 使ったRスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3)
##########################
# 結果比較用の変数を準備
mu <- c() # 標本平均の計算結果
cal.time <- c() # 処理時間
##########################
# 3通りの方法で標本平均を計算
##########################
# 1. コマンドmeanを使用
cal.time[1] <- system.time(mu[1] <- mean(x))[[1]]
# 2. コマンドsumで和を計算して,xの長さで割る
cal.time[2] <- system.time(mu[2] <- sum(x)/length(x))[[1]]
# 3. forループを使って,プログラムっぽく計算 
# forループを使用して関数fを定義
f <- function(x){
 s <- 0
 n <- length(x)
 for(i in 1:n){
  s <- s+x[i]
 }
 s/n
}
cal.time[3] <- system.time(mu[3] <- f(x))[[1]]
########
# 結果の比較
data.frame(method=c("mean","sum/N","for loop"),time=cal.time,mu)

 このRスクリプトを実行すると,処理時間time [s]と標本平均の計算結果muの表が出力されます.私の環境では,以下のようになりました.

    method time       mu
1     mean 0.14 2.999914
2    sum/N 0.08 2.999914
3 for loop 2.11 2.999914

ここで,methodは,それぞれ,上で説明した3通りの方法に対応します.

 当然,標本平均の計算結果muは,どの方法も同じです.しかし,処理時間time [s]は違います.

 意外なのは,meanよりも,sumで計算した後に標本数で割った方が処理時間が短かいことだと思います.半分くらいの時間になるので,

 Rでは,meanよりも,sumの方が速い

ということです.

欠損値NAがあると,処理時間は長くなる

 現実のデータでは,欠損値があることが普通です.データに欠損値NAが含まれると処理時間は長くなります.

 xに欠損値を一つ入れて,上と同じ数値実験をしてみます.ここでは,meansumに,

na.rm=TRUE

のオプションを設定して,欠損値を除きます.

 使ったスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3)
# NAを1つ挿入
x[1] <- NA
##########################
# 結果比較用の変数を準備
mu <- c() # 標本平均の計算結果
cal.time <- c() # 処理時間
##########################
# 3通りの方法で標本平均を計算
##########################
# 1. コマンドmeanを使用
cal.time[1] <- system.time(mu[1] <- mean(x,na.rm=TRUE))[[1]]
# 2. コマンドsumで和を計算して,欠損値を除いたxの長さで割る
cal.time[2] <- system.time(mu[2] <- sum(x,na.rm=TRUE)/(length(x)-sum(is.na(x))))[[1]]
# 3. forループを使って,プログラムっぽく計算 
# forループを使用して関数fを定義
f <- function(x){
 s <- 0
 n <- 0
 for(tmp in x){
  if(!is.na(tmp)){
   s <- s+tmp
   n <- n + 1
  }
 }
 s/n
}
cal.time[3] <- system.time(mu[3] <- f(x))[[1]]
########
# 結果の比較
data.frame(method=c("mean","sum/N","for loop"),time=cal.time,mu)

 このRスクリプトを実行すると,処理時間time [s]と標本平均の計算結果muの表が出力されます.私の環境では,以下のようになりました.

    method  time       mu
1     mean  0.66 2.999827
2    sum/N  0.33 2.999827
3 for loop 12.66 2.999827

ここで,methodは,それぞれ,上で説明した3通りの方法に対応します.

 どの計算方法も,欠損値を気にする必要がなかったケースと比べて計算時間が長くなっています.

 注目すべき点は,今回も,sumを使った標本平均の計算が,最も処理時間が短いということです.meanの半分です.

 おまけで注意しておきますが,否定を使った処理!is.na(x)は,is.na(x)よりも処理時間がかかります.つまり,欠損値がある場合は,

sum(x,na.rm=TRUE)/(sum(!is.na(x)))

ではなく,

sum(x,na.rm=TRUE)/(length(x)-sum(is.na(x)))

とした方が,速いです.

まとめ

 私は18~22歳の間,大学に行かずに,工事現場などで毎日働いていました.その間,学問的なことには全く興味を失っていたので,勉強はしませんでした.その結果,23歳にして中学生レベルの知識しかもたない自分に気がつきました.何の能力もない自分に気づいた後は,勉強して新しいことを知る喜びを感じられたので,誰かに押しつけられるのではなく,自分自身でいろんなことに興味をもって勉強できるようになりました.今となっては人生には無駄な時間はないと感じます.

 人生には無駄な時間はないといっても,Rの処理時間は短い方が良いということが今回のメッセージです.