【Rで時系列解析】FFTを使った離散フーリエ変換の実際（各成分の周波数について）

RでFFT（Fast Fourier Transform）を計算するコマンドfftの結果と，離散フーリエ変換の関係，特に，周波数について説明します．

実数の時系列を $x[i$ ] ( $i=1, 2, 3, \cdots, N$ )とします．時系列解析では， $i=0, 1, 2, \cdots, N-1$ とすることが多い（C言語やpythonも0から）ですが，ここではRに合わせます．つまり，

x[1], x[2], x[3] , ..., x[N]

に実数の時系列データが入っています．

このとき，離散フーリエ変換（サンプリング間隔を1としています）は，
$\displaystyle X(k) = \sum_{j=0}^{N-1} x \left[ j +1 \right] \exp \left( -i 2 \pi j \frac{k}{N} \right)$
あるいは，
$\displaystyle X(k) = \sum_{j=1}^{N} x \left[ j \right] \exp \left( -i 2 \pi (j-1) \frac{k}{N} \right)$
です．私は，物理学科出身なので，虚数単位は $i$ を使っています．

周波数で表したいので， $k$ を周波数に対応させると，
$\displaystyle f_k = \frac{k}{N}$
となります．多くの場合，逆フーリエ変換の式で， $k$ は0から $N-1$ の値をとっているので
$0 \le f_k \le \frac{N-1}{N}$
だね，となります． $\displaystyle f_k$ を使うと離散フーリエ変換は，
$\displaystyle X(f_k) = \sum_{j=0}^{N-1} x \left[ j +1 \right] \exp \left( -i 2 \pi j f_k \right)$
です．しかし，実際はナイキスト周波数というのがあって，サンプリング周波数が1 Hzの場合，意味のある周波数は，
$\displaystyle -\frac{1}{2} \le f_k \le \frac{1}{2}$
です（単位はHzです）．ということで，
$\frac{1}{2} < f_k \le \frac{N-1}{N}$
のときは，1を引いて，
$\frac{1}{2} -1 < f_k \le \frac{N-1}{N} -1$
$-\frac{1}{2} < f_k \le -\frac{1}{N}$
と解釈する必要があります（単位はHzです）．

Rを使って，離散フーリエ変換をしたければ，

X <- fft(x)

を実行するだけです．結果は複素数のベクトルです（Rの世界のベクトルで，一般的には配列です）．

結果は以下のようになります．

$\displaystyle k$	$\displaystyle f_k$	$X(f_k)$ ( $\alpha_k, \beta_k$ は実数)	Rのベクトルの要素
$\displaystyle 0$	$\displaystyle 0$	$\displaystyle X(0) = \alpha_0 + 0 i$	X[1]
$\displaystyle 1$	$\displaystyle \frac{1}{N}$	$\displaystyle X(1) = \alpha_1 + \beta_1 i$	X[2]
$\displaystyle 2$	$\displaystyle \frac{2}{N}$	$\displaystyle X(2) = \alpha_2 + \beta_2 i$	X[3]
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\displaystyle N-1$	$\displaystyle \frac{N-1}{N}$	$\displaystyle X(N-1) = \alpha_{N-1} + \beta_{N-1} i$	X[N]

が，この理解ではダメ（不十分で間違いにつながる）です．

上の表で，フーリエ変換の各要素に対応する周期成分は，
$\displaystyle \begin{eqnarray} x^{(k)} [i] &=& A_k \cos(2 \pi f_k (i-1) + \theta_k ) \quad (i = 1, \cdots, N)\end{eqnarray}$
です．ここで， $A_k$ は $X(f_k)$ の絶対値をNで割ったもの， $\theta_k$ は， $X(f_k)$ の偏角です．

$k=N-m$ とすれば， $\displaystyle f_{N-m} = \frac{N-m}{N}=1-\frac{m}{N}$ なので，
$\displaystyle \begin{eqnarray} x^{(k)} [i] &=& A_k \cos \left(2 \pi \frac{N-m}{N} (i-1) + \theta_k \right) \quad (i = 1, \cdots, N) \\ &=& A_k \cos \left(2 \pi \left( 1 - \frac{m}{N} \right) (i-1) + \theta_k \right) \\ &=& A_k \cos \left(2 \pi (i - 1) + 2 \pi \left( - \frac{m}{N} \right) (i-1) + \theta_k \right) \\ &=& A_k \cos \left(2 \pi \left( - \frac{m}{N} \right) (i-1) + \theta_k \right) \\ \end{eqnarray}$
となります．したがって，周波数は $- \frac{m}{N}$ と解釈することもできます．
$\frac{N-m}{N}$ と $- \frac{m}{N}$ のどちらを選ぶべきかは
$\displaystyle -\frac{1}{2} \le f_k \le \frac{1}{2}$
の範囲に収まっていることで判断します．

下の図は， $N=8$ の場合について，離散フーリエ変換で分解される周期成分を模式的に描いたものです．左側の列に，周波数 $\frac{N-m}{N}$ の周期成分を水色，周波数 $- \frac{m}{N}$ の周期成分を桃色で描いてあります．水色の波形は，サンプリング間隔1とくらべて振動が早すぎる（細かすぎるので），そんな振動があっても離散サンプリングでは正確に捉えきれません．サンプリング結果から，我々の感覚でとらえられる振動は桃色の波形でしょう．図の右側が，離散の点のみを描いたものです．とはいえ，右上図の点の配置から，波を直感的にイメージ（再構成）するのは難しいと思います．

時系列の長さが，奇数と偶数の場合に分けてもう一度まとめると以下のようになります．

奇数のとき（ここでは， $N=2M+1$ とする）

$\displaystyle k$	$\displaystyle f_k$	$X(f_k)$ ( $\alpha_k, \beta_k$ は実数)	Rのベクトルの要素
$\displaystyle 0$	$\displaystyle 0$	$\displaystyle X(0) = \alpha_0 + 0 i$	X[1]
$\displaystyle 1$	$\displaystyle \frac{1}{N}$	$\displaystyle X(1) = \alpha_1 + \beta_1 i$	X[2]
$\displaystyle 2$	$\displaystyle \frac{2}{N}$	$\displaystyle X(2) = \alpha_2 + \beta_2 i$	X[3]
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\displaystyle M-1$	$\displaystyle \frac{M-1}{N}$	$\displaystyle X(M-1) = \alpha_{M-1} + \beta_{M-1} i$	X[M]
$\displaystyle M$	$\displaystyle \frac{M}{N}$	$\displaystyle X(M) = \alpha_M + \beta_M i$	X[M+1]
$\displaystyle M+1$	$\displaystyle -\frac{M}{N}$	$\displaystyle X(M+1) = \alpha_M - \beta_M i$	X[M+2]
$\displaystyle M+2$	$\displaystyle -\frac{M-1}{N}$	$\displaystyle X(M+2) = \alpha_{M-1}- \beta_{M-1} i$	X[M+3]
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\displaystyle N-2$	$\displaystyle -\frac{2}{N}$	$\displaystyle X(N-2) = \alpha_2 - \beta_2 i$	X[N-1]
$\displaystyle N-1$	$\displaystyle -\frac{1}{N}$	$\displaystyle X(N-1) = \alpha_1 - \beta_1 i$	X[N]

偶数のとき（ここでは， $N=2M$ とする）

$\displaystyle k$	$\displaystyle f_k$	$X(f_k)$ ( $\alpha_k, \beta_k$ は実数)	Rのベクトルの要素
$\displaystyle 0$	$\displaystyle 0$	$\displaystyle X(0) = \alpha_0 + 0 i$	X[1]
$\displaystyle 1$	$\displaystyle \frac{1}{N}$	$\displaystyle X(1) = \alpha_1 + \beta_1 i$	X[2]
$\displaystyle 2$	$\displaystyle \frac{2}{N}$	$\displaystyle X(2) = \alpha_2 + \beta_2 i$	X[3]
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\displaystyle M-1$	$\displaystyle \frac{M-1}{2M}$	$\displaystyle X(M-1) = \alpha_{M-1} + \beta_{M-1} i$	X[M+1]
$\displaystyle M$	$\displaystyle \frac{1}{2}$	$\displaystyle X(M) = \alpha_{M} + 0 i$	X[M+2]
$\displaystyle M+1$	$\displaystyle -\frac{M-1}{2 M}$	$\displaystyle X(M+1) = \alpha_{M-1}- \beta_{M-1} i$	X[M+3]
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\displaystyle N-2$	$\displaystyle -\frac{2}{N}$	$\displaystyle X(N-2) = \alpha_2 - \beta_2 i$	X[N-1]
$\displaystyle N-1$	$\displaystyle -\frac{1}{N}$	$\displaystyle X(N-1) = \alpha_1 - \beta_1 i$	X[N]

私が伝えたいことが分かりますか．注目してほしいのは， $\displaystyle k$ 番目の成分の周波数，同じ周波数 $\displaystyle f_k$ がどの位置に何回登場するか， $X(f_k)$ の虚部の符号の関係です．
この関係を理解しておかないと，FFTを使いこなすことは難しいと思います．

実際にRのfftを使って説明した，以下の記事も参考にしてください．
chaos-kiyono.hatenablog.com

例えば，以下の記事で使っているRスクリプトを理解するには，ここで説明したことを理解する必要があります．
chaos-kiyono.hatenablog.com