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

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

【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は,フーリエ変換の選択肢になるので,その理屈と特徴を理解しておくことが大切だと思います.