B. フーリエ解析 と スペクトル解析¶
本章では、周波数解析の代表的手段であるフーリエ解析およびスペクトル解析について解説する。 なお、１次元の実数データについて考えることにする。
B-1. フーリエ（調和）解析¶
ある決まった周期成分の振幅や位相を求めることをフーリエ解析（調和解析）という。年周変動や日周変動への適用がその代表例である。また空間方向について特定の波数成分を抽出して解析する際にも用いられる（プラネタリー波の解析など）。
(1) 三角関数を用いた表現¶
時系列データは、様々な周期の波成分の重ね合わせと考えることができる． まずは時系列データ$f(x_k)\ (x_k = 2\pi/N; k=0,...N-1)$を以下のように三角関数で展開することを考えよう：
$N$が偶数: $N=2n$ $$f(x_k) = \frac{A_0}{2} + \sum_{j=1}^{n-1} \left( A_j \cos{\frac{2\pi jk}{N}} + B_j \sin{\frac{2\pi jk}{N}} \right) + \frac{A_n}{2} \cos{\frac{2\pi nk}{N}} \ \ (k=0,..,N-1) \tag{1}$$
$N$が奇数: $N=2n+1$ $$f(x_k) = \frac{A_0}{2} + \sum_{j=1}^{n} \left( A_j \cos{\frac{2\pi jk}{N}} + B_j \sin{\frac{2\pi jk}{N}} \right) \tag{2}$$
このとき、波数$j$成分の係数は
\begin{align} A_j &= \frac{2}{N} \sum_{k=0}^{N-1} f(x_k) \cos{\frac{2\pi jk}{N}} \tag{3} \\ B_j &= \frac{2}{N} \sum_{k=0}^{N-1} f(x_k) \sin{\frac{2\pi jk}{N}} \tag{4} \end{align}として求められる（証明は下）。
（証明）Nが偶数の場合の$A_j$についてのみ示す. 式(3)の右辺(r.h.s.)に、式(1)の表現$f(x_k)=...$を代入すると \begin{align} r.h.s. &= \frac{2}{N}\sum_{k=0}^{N-1} \bigl (\frac{A_0}{2} + \sum_{l=1}^{n-1} A_l \cos{\frac{2\pi lk}{N}} + \frac{A_n}{2}\cos{\frac{2\pi n k}{N}}\bigr )\cos{\frac{2\pi jk}{N}} \\ &=\frac{2}{N} \bigl (\frac{A_0}{2} \sum_{k=0}^{N-1} \cos{\frac{2\pi jk}{N}} + \sum_{l=1}^{n-1} A_l \sum_{k=0}^{N-1} \cos{\frac{2\pi lk}{N}}\cos{\frac{2\pi jk}{N}} + \frac{A_n}{2}\sum_{k=0}^{N-1} \cos{\frac{2\pi nk}{N}}\cos{\frac{2\pi jk}{N}} \bigr ) \end{align} ここで \begin{align} \sum_{k=0}^{N-1} \cos{ \frac{2\pi lk}{N} } \cos{\frac{2\pi jk}{N}} &= \sum_{k=0}^{N-1} \left( \cos{ \frac{2\pi(l+j)k}{N} } + \cos{ \frac{2\pi(l-j)k}{N}} \right) \\ &= \begin{cases} 0 & \text{$l+j \ne mN $ かつ $ l-j \ne m'N$} \\ \frac{N}{2} & \text{$l+j = mN $ もしくは $ l-j = m'N$} \\ N & \text{$l+j = mN $ かつ $ l-j = m'N$} \end{cases} \end{align} を用いると($m, m'$は整数）、左辺に等しいことが示せる（途中、$j \ne N/2$ では $l+j<N$であることに注意する）。
（注）上式から$j=0$、$j=n$の場合のみ特別に考えねばならない理由が理解できると思う．すなわち、式(1)、(2)において$A_0, A_n$だけ"ハズレ"て定義されれている．
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
#from scipy.fftpack import rfft
#from scipy.fftpack import irfft
from scipy.fftpack import fft
from scipy.fftpack import ifft
pi = np.pi
N = 1000 #データ数
dT = 1.0 #間隔
A = 2.0; B = 1.0
Ts = 100 # 周期 = 100 day のサンプルデータを作る
x = np.linspace(0.0,(N-1)*dT,N)
y = A*np.cos(2.0*pi*x/Ts) + B*np.sin(2.0*pi*x/Ts)
#サンプル時系列の描画
fig = plt.figure(figsize=(10,3))
ax = fig.add_subplot(111)
ax.plot(x, y)
plt.show()
式(3)-(4)に従ってこの時系列にフーリエ解析を施し、周期$T_s$のcos, sinの係数を推定してみる。
fc=np.cos(2.0*pi*x/Ts); fs=np.sin(2.0*pi*x/Ts)
Ainv = np.sum(y*fc)/N*2.0
Binv = np.sum(y*fs)/N*2.0
print("Calculated amplitude: A, B = ",Ainv, Binv)
print("True amplitude: A, B= ", A, B)
(2) 複素数を用いた表現¶
上記の$\cos$、$\sin$関数の係数をまとめて扱うために、複素数による表現を用いるのが便利である。以下では、$N$:偶数のときを考える（奇数のときも同様に示せる）。
まず \begin{align} \cos \theta &= \frac{e^{i\theta} + e^{-i\theta}}{2} \tag{9} \\ \sin \theta &= \frac{e^{i\theta} - e^{-i\theta}}{2i} \tag{10} \end{align}
を用いて
\begin{align} f(x_k) &= \frac{A_0}{2} + \sum_{j=1}^{n-1} \left( A_j \frac{e^{i2\pi j k/N} + e^{-i2\pi j k/N}}{2} + B_j \frac{e^{i2\pi j k/N} - e^{-i2\pi j k/N}}{2i} \right) + \frac{A_n}{2}\frac{e^{i2\pi n k/N} + e^{-i2\pi n k/N}}{2} \\ &= \frac{A_0}{2} + \sum_{j=1}^{n-1} \left( \frac{A_j - i B_j}{2} e^{i2\pi j k/N} + \frac{ A_j + i B_j}{2} e^{-i2\pi j k/N} \right) + \frac{A_n}{2}\frac{e^{i2\pi n k/N} + e^{-i2\pi n k/N}}{2} \\ &= c_0 + \sum_{j=1}^{n-1} \left( c_j e^{i2\pi j k/N} + c_j^* e^{-i2\pi j k/N} \right ) + \frac{1}{2} (c_ne^{i2\pi n k/N} + c_n^*e^{-i2\pi n k/N})\tag{11} \end{align}と書き直す.ここで \begin{align} c_j = \frac{A_j - iB_j}{2} \tag{12} (B_0 = 0,\ B_n = 0) \end{align} と置いた．
また、$\cos$, $\sin$関数の対称性と周期性に注目すると、 \begin{align} A_{j} &= A_{-j} = A_{N-j} \tag{13} \\ B_{j} &= -B_{-j} = -B_{N-j} \tag{14} \\ \end{align} であるから \begin{align} c_j^* = \frac{A_j + iB_j}{2} = \frac{A_{N-j} - iB_{N-j}}{2} = c_{N-j} \tag{15} \end{align}
とかける。したがって、
\begin{align}
f(x_k) &= c_0 + \sum_{j=1}^{n-1} \left( c_j e^{i2\pi j k/N} + c_j^* e^{-i2\pi j k/N} \right) + \frac{1}{2} (c_ne^{i2\pi n k/N} + c_n^*e^{-i2\pi n k/N})\\
&= c_0 + \sum_{j=1}^{n-1} \left( c_j e^{i2\pi j k/N} + c_{N-j} e^{i2\pi (N-j) k/N} \right) + \frac{1}{2} (c_ne^{i2\pi n k/N} + c_{N-n}e^{-i2\pi (N-n) k/N})\\
&= c_0 + \sum_{j=1}^{n-1} \left( c_j e^{i2\pi j k/N} + c_{N-j} e^{i2\pi (N-j) k/N} \right) + c_ne^{i2\pi n k/N} \\
%&= c_0 + \sum_{j=1}^{N-1} \left( c_j e^{-i2\pi j k/N} \right) \ \ \\
&= \sum_{j=0}^{N-1} \left( c_j e^{i2\pi j k/n} \right) \ \ \tag{16}
\end{align}
と表せる。これをフーリエ逆変換と言う。
一方、フーリエ正変換は
$$c_j = \frac{1}{N}\sum_{k=0}^{N-1} f(x_k) e^{-i2\pi jk/N} \ \ \ (j=0,..,N-1) \tag{17}$$
であり、これにより各波数の振幅（$c_j$: 複素数であることに注意！）が求まる（証明は以下）。
（証明: 右辺が左辺に等しいことを示す） \begin{align} \frac{1}{N}\sum_{k=0}^{N-1} f(x_k) e^{-i2\pi jk/N} &= \frac{1}{N}\sum_{k=0}^{N-1} (\sum_{m=0}^{N-1} c_m e^{i2\pi mk/N}) e^{-i2\pi jk/N} \\ &= \sum_{m=0}^{N-1} c_m (\frac{1}{N}\sum_{k=0}^{N-1} e^{-i2\pi (j-m)k/N}) \\ &= \sum_{m=0}^{N-1} c_m \delta_{j,m} \end{align} ここで以下の関係式（等比級数の公式から示せる）を用いた: \begin{equation} \frac{1}{N} \sum_{k=0}^{N-1} e^{i2\pi(j-m)k/N} = \begin{cases} 1 & j \equiv m \ \ (\text{mod}\ N) \\ 0 & j \not\equiv m \ \ (\text{mod}\ N) \end{cases} \end{equation}
まとめ¶
以上まとめると,
フーリエ正変換:
\begin{align} c_j = \frac{1}{N}\sum_{k=0}^{N-1} f(x_k) e^{-i2\pi jk/N} \ \ \ (j=0,..,N-1) \tag{18} \end{align}フーリエ逆変換: \begin{align} f(x_k) = \sum_{j=0}^{N-1} c_j e^{i2\pi j k/N} \tag{19} \end{align}
（注意） 正変換・逆変換の定義は（前にかかるファクターを含めて）いろいろと流儀があるので注意されたい（今回用いるPythonのパッケージの定義もこれと異なる；式(29)-(30)を参照）。
（補足）¶
式(16)の第一式に注目すると$j$成分の変動は \begin{align} A_j \cos({2\pi jk/N}) + B_j \sin({2\pi jk/N}) = \text{Re} \left( 2c_j e^{i2\pi j k/N} \right) = \text{Re} \left( (A_j-B_j) e^{i2\pi j k/N} \right) \end{align} と表せることがわかる。
したがって、$j=0...n(=N/2)$の和で \begin{align} f(x_k) &= \frac{A_0}{2} + \sum_{j=1}^{n} \left( c_j e^{i2\pi j k/N} + c_j^* e^{-i2\pi j k/N} \right) \\ &= \text{Re} \sum_{j=0}^{n} \left( C_j e^{i2\pi j k/N} \right)\ \ (C_j = 2c_j (j\ne 0), C_j = c_j (j=0)) \end{align} とも表せる。
（実践編）¶
サンプルデータを使って、複素数を用いたフーリエ解析を実践してみる. 式(12)に注目すると、周波数$j$成分の$\cos$、$\sin$の係数は
\begin{align} A_j &= 2\ \text{Real}(c_j) \tag{20} \\ B_j &= - 2\ \text{Imag}(c_j) \tag{21} \end{align}から求められる。
img = 1j #pythonでの複素数の表現 i
fexp=np.exp(-img*2.0*pi*x/Ts)
cinv = np.sum(y*fexp)/N
Ainv, Binv = 2*cinv.real, -2*cinv.imag
print("Calculated amplitude = ",Ainv, Binv)
print("True amplitude = ", A, B)
B-2. スペクトル解析¶
（概要）¶
実際のデータ解析において、予め現象の周期$(T_s)$が分かっていて（日変化、季節変化等）その振幅・位相を求めたい場合には、 前章のように$T_s$についてフーリエ解析を施せば良い．一方、あるデータにおいて卓越する周波数帯を見い出したい場合には、全ての 周波数についてフーリエ解析を行い、各成分の大きさを比べることになる。これをスペクトル解析と言う。
（スペクトルの表現）¶
各周波数成分のの大きさとしては、各周波数成分のエネルギーを基準とするのが合理的である。
まずはその準備にあたってパーセバルの式を紹介する。
これは、データの分散（エネルギーと考えて良い）が、波数の二乗和に等しいことを示すものである。
パーセバルの式¶
\begin{equation} \frac{1}{N}\sum_{k=0}^{N-1}|f(x_k)|^2 = \sum_{j=0}^{N-1}|c_j|^2 \tag{22} \end{equation}（証明） \begin{align} \frac{1}{N}\sum_{k=0}^{N-1}f(x_k) f(x_k)^* & = \frac{1}{N} \sum_{k=0}^{N-1}f(x_k) \left( \sum_{l=0}^{N-1} c_l^* e^{-i2\pi lk/n} \right) \\ & = \sum_{l=0}^{N-1}c_l^* \left( \sum_{k=0}^{N-1} f(x_k) e^{-i2\pi lk/N} \right) \\ & = \sum_{l=0}^{N-1}c_l^* c_l \end{align}
なお、$c_j$の定義から、調和関数の振幅($A_j, B_j$）との間には \begin{align} |c_j|^2 = \frac{A_j^2 + B_j^2}{4} \tag{23} \end{align} の関係があることにも注意しておく。
パワースペクトル密度関数¶
実際のデータ解析では($|c_j|^2$ではなく）以下のように定義したパワースペクトル密度(Power Spectrum Density: PSD)を図示することが多い。
\begin{align}
P_s(\omega_j) = \frac{|c_j|^2}{\Delta \omega} \tag{24}
\end{align}
ここで$\Delta \omega$は周波数分解能で \begin{align} \Delta \omega = \frac{2\pi}{N dT} (\text{rad}\ s^{-1}) \tag{25} \end{align} と表せる($N$:データの長さ, $dT$:データの時間分解能)。
元のスペクトルを最小周波数で割っているだけだが、PSDは(元の時系列)$^{-1}$の次元を持ち(すなわり、"単位周波数あたり"）
\begin{align}
\sum_{j=0}^{N-1}|c_j|^2 = \sum_{j=0}^{N-1} P_s(\omega_j) \Delta \omega = \int_{-\omega_{max}}^{\omega_{max}} P_s(\omega) d\omega \tag{26}
\end{align}
となって、周波数空間でのPSDの積分は元データの分散に一致する（むしろ、そうなるようにPSD($P_s(\omega)$)を定義した）。
すなわち、$P_s(\omega)\Delta \omega $は全分散（平均パワー）への"$\omega \pm \Delta \omega/2$"の周波数帯からの寄与率を表す。
両側スペクトルと片側スペクトル¶
ここで、上式において$c_j = c_{N-j}^*$であることに注意しよう。
すなわち、$|c_j|^2$(あるいはPSD)を図示すると、$\frac{N}{2}$に関して対称になっている（両側スペクトルと言う）。
そこで、結果を図示するにあたってはj = 0..N/2 の範囲で以下のように定義した片側スペクトルを使うことが多い。
このとき、パーセバルの式は \begin{align} \sum_{j=0}^{N-1}|c_j|^2 = \int_{0}^{\omega_{max}} G_s(\omega) d\omega \tag{28} \end{align} となる．
（実践）¶
では、サンプル時系列データを使って実際にスペクトル解析を実践してみよう。 効率よく計算するアルゴリズムとしてはFFT（高速フーリエ変換）がある。ここではSciPyのFFTライブラリ関数を用いてこれを実装する。
Python(Scipy)のライブラリ関数¶
フーリエ正変換
(ライブラリ関数fft): \begin{align} c_j = \sum_{k=0}^{N-1} f(x_k) e^{-i2\pi jk/N} \ \ \ (j=0,..,N-1) \tag{29} \end{align}フーリエ逆変換
(ライブラリ関数ifft): \begin{align} f(x_k) = \frac{1}{N}\sum_{j=0}^{N-1} c_j e^{i2\pi j k/N} \tag{30} \end{align}
と定義されている．
ここで、これらを式(18)-(19)と比較すると1/Nのファクターが定義とは逆の変換に付いているので注意が必要．
(参考) ScipyのFFTライブラリの説明 (https://docs.scipy.org/doc/scipy/reference/fftpack.html)
# フーリエ(逆)変換
c = fft(y)/N
pw = abs(c)**2 #abs(c[0:int(N/2)])
# Power Spectrum Density (PSD)
domg = 2.0*pi/(N*dT)
psd = pw/domg
psd_oneside = psd.copy()*2
psd_oneside[int(N/2)] /= 2.0
psd_oneside[0] /= 2.0
# 0,N/2は二倍しなくて良いので元に戻す（ほとんど寄与しないが）
# 周波数スペクトル
xfft = np.linspace(0.0,domg*(N/2),int(N/2)+1) # (frequency times s-1)
xfft = 2*pi/xfft # (period s)
fig = plt.figure(figsize=(10,3))
ax = fig.add_subplot(111)
print(xfft.shape, pw.shape)
ax.plot(xfft[1:int(N/2)+1],psd_oneside[1:int(N/2)+1]) #バグ修正(pw->psd_oneside): 2019/04/24
#（注意）a[start:end]としたとき、endは含まない
plt.show()
# 元データの分散 (パーセバルの式を確かめる)
var = (y**2).sum()/N
# スペクトルの周波数積分が全分散に一致することを確かめる
print("元データの分散 = ", var)
print("スペクトルの和 = ", pw.sum())
print("PSDの周波数積分 = ", (psd*domg).sum())
（補足）フーリエ変換/逆変換で元に戻ることの確認¶
逆変換 -> 順変換 -> 逆変換 として元の時系列に戻ることの確認。
なお、この性質を用いてフィルターを作ることもできる（周波数空間で、適当な周波数帯のパワーのみ残したものを逆変換すればよい）。
# フーリエ逆変換
c = fft(y)/N
# フーリエ順変換
yinv = ifft(c)*N
# plot
fig = plt.figure(figsize=(10,3))
ax = fig.add_subplot(111)
#元の時系列
ax.plot(x, y)
#フーリエ逆変換->順変換 として戻ってきた時系列（元の時系列と一致するはず)
ax.plot(x,yinv,color='r')
plt.show()
