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
> x <- rnorm(N)
> system.time(X <- fft(x))
ユーザ システム 経過
0.01 0.00 0.010.01 秒でした．
遅くなるのは，素因数分解できないときです．時系列の長さを素数にしてみると，
> N <- 99991
> x <- rnorm(N)
> system.time(X <- fft(x))
ユーザ システム 経過
8.25 0.00 8.25 8.25秒になりました．時系列の長さは，前の例よりも短いですが，計算時間は800倍くらい長くかかりました．
Nが素数のときは，Rが使っているfftのアルゴリズムは，FFTとは呼べなくなります．Nが素数の例は極端ですが，2のべき乗じゃなくても，ある程度，素因数分解できれば計算時間はかなり短くなります．
例えば，時系列の長さを素数から
つ除いて，時系列の長さを
にすれば，
> 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
