畳庵〜tatamiya practice〜

機械学習・統計モデリング練習帳。読書・source code読みの記録や，実装など。

NumPyでフーリエ変換を計算するためのメモ〜sin, cosで展開したときの係数をどう読み取るか〜

時系列分析 信号処理 Python

はじめに

内容

  • numpy.fft.fft/rfftの返り値は何を表しているか？
  • numpy.fft.fft/rfftの出力値から、\sin, \cosでFourier変換した表現に書き直すにはどうすればいいか？
    • その際、データ数の偶奇で出てくる項・和の範囲が異なるのは何故か？

背景

  • Fourier変換は時系列データの周期性を調べるときによく使われる
    • 時系列を三角関数の重ね合わせとして表現する
  • PythonでFourier変換を計算してみようとすると、引っかかる点が多い
    • numpy.fftは指数関数表現をもとにしている
    • 指数関数表現だと複素数が出てきてイメージがしにくい
    • 指数関数表現と三角関数表現の対応がわかりにくい
      • 周波数成分ごとに和をとる際、両者で和の範囲が異なる
      • 三角関数表現では、入力データ数の偶奇で処理の仕方が異なる
  • 三角関数表現をイメージをしながらnumpy.fftを触れるように、対応関係を整理したい。

TL;DR

  • 長さTの時系列データ\left \lbrace y_t \right \rbrace_{t=0}^{T-1}をnumpy.fft.fftでFourier変換した際の返り値の中身は、以下に示す指数関数表現の係数である：
{ \begin{aligned} c_k &= \sum_{t=0}^{T-1} y_t \exp \left(- i \frac{2\pi k}{T} t \right) \\ \end{aligned} }
  • 三角関数を用いて展開した際の係数との対応は、以下のようになっている：
{ \begin{aligned} y_t &= c + \sum_{k=1}^{M} a_k \cos\left( \frac{2\pi k}{T} t\right) + \sum_{k=1}^{M} b_k \sin\left( \frac{2\pi k}{T} t\right) \\ \end{aligned} }
{ \begin{aligned} c = \frac{1}{T} (c_0 + c_{T/2}) \end{aligned} }
{ \begin{aligned} a_k = \frac{2}{T} \,{\rm Re}\,c_k, \quad b_k = - \frac{2}{T} \,{\rm Im}\,c_k\\ \end{aligned} }

ただし、

{ \begin{aligned} M= \begin{cases} (T-1)/2 \quad & T: {\rm 奇数}\\ T/2 -1 \quad & T: {\rm 偶数}\\ \end{cases} \end{aligned} }
{ \begin{aligned} c_{T/2}= \begin{cases} 0 \quad & T: {\rm 奇数}\\ \sum_{t=0}^{T-1}(-1)^ty_t \quad & T: {\rm 偶数}\\ \end{cases} \end{aligned} }
  • 三角関数表現にてデータ長Tの偶奇で計算方法が異なるのは、和\sum _ {k=0}^{T-1}をとる際に複素共役の関係c _ {T-k} = c_k^*を利用してT/2で折り返しc_kとc_{T-k}でペアを作るためである。
  • 入力データが実数であれば、numpy.fft.rfftを使って計算するのが便利である

序章：Fourier変換とはどのようなものか？

長さTの時系列データ\left \lbrace y_t \right \rbrace_{t=0}^{T-1}が与えられた時、 （離散）Fourier変換では以下のように時系列を三角関数を使って書き表します1：

{ \begin{aligned} y_t &= c + \sum_{k=1}^{M} a_k \cos\left( \frac{2\pi k}{T} t\right) + \sum_{k=1}^{M} b_k \sin\left( \frac{2\pi k}{T} t\right) \\ \end{aligned} \tag{1} \label{sincos_fourier} }

\sin\left( \frac{2\pi k}{T} t\right), \cos\left( \frac{2\pi k}{T} t\right)はともに周期T/kの関数ですので、ざっくりいうといかなる時系列も周期関数に分解できる、ということを意味します。

この性質を利用して、Fourier変換はある時系列がどのような周期を持っているかを調べるときに使われます。

例えば、以下のような長さT=100の時系列があったとします：

f:id:tatamiyatamiyatatatamiya:20201226145437p:plain

全体的に周期50の\sinカーブを描いていますが、その周りで細かい動きも乗っています。 この細かい動きは単なるノイズなのか、またはなんらかの周期を持っているかまではよくわかりません。

ここでFourier変換を行い、さらに周期T/kを持つ項の係数の大きさに着目するために、パワースペクトルと呼ばれる\frac{1}{2} (a _ {k}^2 + b _ {k}^2)を見てやると、以下のようになります：

f:id:tatamiyatamiyatatatamiya:20201226145504p:plain

この結果ではk=2の強いピークに加えて、k=10のところにもピークが出ていますので、周期100/2=50の成分に加えて、周期100/10=10の成分も持っている、と理解できます。

後出しになりますが、先ほどの時系列は、以下のように周期50のsin波と周期10のcos波を重ね合わせたものにノイズを加えて生成しました。したがって、上記のパワースペクトルの結果は、周期性をうまく捉えられているといえそうです。

t = np.arange(100)
y = 2.0 * np.sin(2*np.pi * t / 50) + 1.0 * np.cos(2 * np.pi * t / 10) + 0.8 * np.random.randn(t.shape[0])

なお、式(\ref{sincos_fourier})の中で、\sumの範囲でMという変数を用いましたが、こちらは以下のように入力データの長さTの偶奇によって変わります：

{ \begin{aligned} M= \begin{cases} (T-1)/2 \quad & T: {\rm 奇数}\\ T/2 -1 \quad & T: {\rm 偶数}\\ \end{cases} \end{aligned} }

どちらの場合も「T/2未満の最大の整数」となっています。 何故このようになるかについては、後ほど解説します。

NumPyでFourier変換を行う方法

Fourier変換を実際にプログラミングで計算してみたい、という場合には、PythonのNumPyパッケージにあるnumpy.fftを使うのが便利です。

ただし、使う際にいくつか気を付けるべきことがあります。

試しに先ほどの例で示した例を、関数numpy.fft.fftにそのまま突っ込んでみます。

fy = np.fft.fft(y)
print(fy.shape)
print(fy)

"""
(100,)
[-5.09497885e+00+0.00000000e+00j -1.47588178e+00-7.84995285e+00j
  2.48806875e+00-9.64778465e+01j  3.05813510e+00-1.08640360e+01j
 -1.88409511e+00-4.69607698e+00j -4.18655969e-01+8.53439687e+00j
  4.24523418e+00+8.25281748e+00j -5.00154635e+00+3.56548344e+00j
  ...
"""

このように、複素数のよくわからない配列が返ってくるのです。

これの意味を理解し、先ほどの式(\ref{sincos_fourier})で示した三角関数を用いた表現と対応付けるには、

  • Fourier変換の指数関数表現とは何か
  • 指数関数表現と三角関数表現はどのように対応付くか

を知る必要があります。

numpy.fft.fftの出力と指数関数を使った表現

Fourier変換には、三角関数を用いた式(\ref{sincos_fourier})とは別に指数関数を用いた表現方法もあります。

そして、先ほどのnp.fft.fftで出てきた複素数の配列の正体は、下記の指数関数表現の係数c_kがそのまま対応しています：

{ \begin{aligned} y_t &= \frac{1}{T} \sum_{k=0}^{T-1} c_k \exp \left( i \frac{2\pi k}{T} t \right)\\ c_k &= \sum_{t=0}^{T-1} y_t \exp \left(- i \frac{2\pi k}{T} t \right) \end{aligned}\tag{2} \label{exp_fourier} }

いきなり指数関数\expと虚数iが出てきましたが、\exp(\pm i\frac{2\pi k}{T}t)の部分はオイラーの公式を用いて

{ \begin{aligned} \exp \left(\pm i\frac{2\pi k}{T}t\right) = \cos\left(\frac{2\pi k}{T}t\right) \pm i\sin\left(\frac{2\pi k}{T}t\right) \end{aligned} \tag{3} \label{euler} }

のように表すことができます（オイラーの公式の簡単な導出は、最後に補足します）。

これを用いれば三角関数を使って書き直すことができるものの、複素数が含まれているため、 実数である式(\ref{sincos_fourier})の係数a_k, b_kとは このままでは対応がつきません。 さらに、\sumで和をとる範囲も2つの表現で異なっています。

指数関数表現と三角関数表現の対応

では、冒頭の式(\ref{sincos_fourier})で示した三角関数による表現と、先ほどの式(\ref{exp_fourier})の指数関数表現とは、どう対応しているのでしょうか？

結論からいくと、以下のように対応します：

{ \begin{aligned} c = \frac{1}{T} (c_0 + c_{T/2}) \end{aligned} }
{ \begin{aligned} c_{T/2}= \begin{cases} 0 \quad & T: {\rm 奇数}\\ \sum_{t=0}^{T-1}(-1)^ty_t \quad & T: {\rm 偶数}\\ \end{cases} \end{aligned} }
{ \begin{aligned} a_k = \frac{2}{T} \,{\rm Re}\,c_k, \quad b_k = - \frac{2}{T} \,{\rm Im}\,c_k\\ \end{aligned} }

これを以下で示します。

先ほど式(\ref{euler})でお見せしたオイラーの公式を、式(\ref{exp_fourier})の指数関数表現に代入して、\sin, \cosに分けて行きます。

その際に、以下の関係に着目します：

{ \begin{aligned} \cos \left(\frac{2\pi (T-k)}{T} t \right) &=\cos \left(-\frac{2\pi k}{T} t \right)\\ &=\cos \left(\frac{2\pi k}{T} t \right)\\ \sin \left(\frac{2\pi (T-k)}{T} t \right) &= \sin \left(-\frac{2\pi k}{T} t \right)\\ &= -\sin \left(\frac{2\pi k}{T} t \right). \end{aligned} }

これを用いて \sum_{k=0}^{T-1}和をとる際にT/2で折り返し、kとT-kをペアにすると、

{ \begin{aligned} y_t &= \frac{1}{T} \sum_{k=0}^{T-1} c_k \exp \left( i \frac{2\pi k}{T} t \right)\\ &= \frac{1}{T} \sum_{k=0}^{T-1} c_k \cos \left(\frac{2\pi k}{T} t \right) + i\frac{1}{T} \sum_{k=0}^{T-1} c_k \sin \left( \frac{2\pi k}{T} t \right)\\ &= \frac{1}{T} c_0 + \frac{1}{T}c_{T/2}\\ &\qquad + \frac{1}{T} \sum_{k=1}^{M} (c_k + c_{T-k} )\cos \left(\frac{2\pi k}{T} t \right)\\ &\qquad + i\frac{1}{T} \sum_{k=1}^{M} (c_k - c_{T-k}) \sin \left( \frac{2\pi k}{T} t \right) \end{aligned} }

のように書き表すことができます。

このとき、k=0についてはペアを作れないため、c_0がそのまま残ります。

さらに、時系列の長さTが偶数の際にも、k=T/2がペアを作れないためc_{T/2}がそのまま出てきます。

この偶奇による違いを理解するには、円周上で(\cos \left( \frac{2\pi k}{T} t \right), \sin \left( \frac{2\pi k}{T} t \right))を表示してみるとわかりやすいです：

f:id:tatamiyatamiyatatatamiya:20201226145543p:plain

以上から、

{ \begin{aligned} c = \frac{1}{T} (c_0 + c_{T/2}) \end{aligned} }
{ \begin{aligned} a_k = \frac{1}{T}(c_k + c_{T-k}), \quad b_k = i\frac{1}{T} (c_k - c_{T-k}) \end{aligned} \tag{4} \label{sincos_coefficients} }

と置くことで、式(\ref{sincos_fourier})のように表せることがわかりました。

依然として係数a_k, b_kが複素数のように見えますが、入力時系列\left \lbrace y_t \right \rbrace_{t=0}^{T-1}が実数（たいていの場合実数ですが）であれば、a_k, b_kも実数となります。

何故なら、c_kとc_{T-k}とは以下のように複素共役の関係にあるからです：

{ \begin{aligned} c_{T-k} &=\sum_{t=1}^{T-1}y_t\exp \left( -i \frac{2\pi (T-k)}{T} t \right)\\ &= \sum_{t=1}^{T-1}y_t\exp \left( +i \frac{2\pi k}{T} t \right) \\ &= c_k^*. \end{aligned}\tag{4} \label{conjugate} }

途中、tが整数なので\exp(\theta -i2\pi t) = \exp(\theta)という関係を使いました。

この関係を用いると、a_k, b_kは以下のように書き直せます、両者が実数となることがわかります：

{ \begin{aligned} a_k &= \frac{1}{T} 2 \,{\rm Re}\,c_k \\ &=\frac{2}{T}\sum_{t=0}^{T-1} y_t \cos\left( \frac{2\pi k}{T} t\right)\\ b_k &= - \frac{1}{T} 2\,{\rm Im}\,c_k\\ &=\frac{2}{T}\sum_{t=0}^{T-1} y_t \sin\left( \frac{2\pi k}{T} t\right). \end{aligned} }

（係数1 / Tに2倍かかることに気をつけてください。）

なお、c_kとa_k, b_kの間には、

{ \begin{aligned} c_k = \frac{T}{2} (a_k - i b_k) \end{aligned} }

という関係があることもわかります。

numpy.fftを使って三角関数表現のFourier変換を行う

以上をもとに、numpy.fftを使って時系列をを\sin, \cosでFourier変換します。

前章ではnumpy.fft.fftという関数を使いましたが、 今回はnumpy.fft.rfftという別の関数を使います。

これは入力データが実数のときに使える関数で、 式(\ref{conjugate})に示したc_{T-k}=c_k^*という複素共役の関係に着目し、numpy.fft.fftの出力のうち前半部分k=1, 2, ..., M (, T/2)だけを取り出したものです。

これを使うと、冒頭の例でお見せしたデータは以下のようにFourier変換することができます：

T = y.shape[0]

if T %2 == 0:
    c = (ck[0].real + ck[-1].real) / T
    ak = ck[1:-1].real * 2 / T
    bk = -ck[1:-1].imag * 2 / T
else:
    c = ck[0].real / T
    ak = ck[1:].real * 2 / T
    bk = -ck[1:].imag * 2 / T

y_sum = c * np.ones_like(t)

for k, (a, b) in enumerate(zip(ak, bk), 1):
    
    y_k = a * np.cos(2*np.pi*k / T * t) + b * np.sin(2*np.pi*k / T * t)
    y_sum += y_k

各周波数成分と、あるkまでの累積和をプロットすると以下のようになります。

f:id:tatamiyatamiyatatatamiya:20201226145358p:plain

おわりに

本記事では、numpy.fftの出力値の読み取りかたと、三角関数表現への変換方法を解説しました。

冒頭で触れたパワースペクトルの計算方法には、時系列\lbrace y_t\rbrace_{t=0}^{T-1}そのものをFourier変換するほかに、自己共分散をFourier変換する方法もあります。

また、実際にパワースペクトルを計算する際には、1からFourier変換を行わずに、SciPyの専用の関数2を使うことができます。

その辺りについても、また記事として書きたいと思います。

最後になりますが、今回の計算・作図を行った際のコードは、こちらに公開しています。

補足

係数1/Tのかけかたのパターン

Fourier変換の際に、1/Tをどこでかけるかについては、以下の3パターンの流儀があります。

{ \begin{aligned} {\rm Pattern \,1.}\\ y_t &= \frac{1}{T} \sum_{k=0}^{T-1} c_k \exp \left( i \frac{2\pi k}{T} t \right), \, & c_k &= \sum_{t=0}^{T-1} y_t \exp \left(- i \frac{2\pi k}{T} t \right) \\ {\rm Pattern \,2.}\\ y_t &= \frac{1}{\sqrt{T}} \sum_{k=0}^{T-1} c_k \exp \left( i \frac{2\pi k}{T} t \right), \, & c_k &= \frac{1}{\sqrt{T}} \sum_{t=0}^{T-1} y_t \exp \left(- i \frac{2\pi k}{T} t \right) \\ {\rm Pattern \,3.}\\ y_t &= \sum_{k=0}^{T-1} c_k \exp \left( i \frac{2\pi k}{T} t \right), \, & c_k &= \frac{1}{T} \sum_{t=0}^{T-1} y_t \exp \left(- i \frac{2\pi k}{T} t \right) \\ \end{aligned} }

いずれの場合も最終的に得られる結果は同じです。

本記事では、NumPyの仕様に合わせてPattern 1で書きました。

人や分野によってどのパターンを用いるかが異なってくるので3、専門書を読み比べる時や複数のライブラリを扱う際にはお気をつけください。

Eulerの公式の簡易な導出

式(\ref{euler})で用いたEulerの公式

{ \begin{aligned} e^{i\theta} = \cos \theta + i \sin \theta \end{aligned} }

について、Tayler展開を用いた簡単な導出を示します。

{ \cos\theta = 1 - \frac{1}{2!}\theta^2 + \frac{1}{4!}\theta^4-\frac{1}{6!}\theta^6+... }
{ \sin \theta = \theta - \frac{1}{3!}\theta^3 + \frac{1}{5!}\theta^5-\frac{1}{7!}\theta^7+... }

さらに，e^{x}のTayler展開はe^{x}=\sum_{k=0}^\infty x^{k}/k!であることを利用し，

{ \begin{aligned} e^{i\theta} &= 1 + i\theta + \frac{1}{2!}(i\theta)^2 + \frac{1}{3!}(i\theta)^3 + \frac{1}{4!}(i\theta)^4+\frac{1}{5!}(i\theta)^5+...\\ &= \left( 1 - \frac{1}{2!}\theta^2 + \frac{1}{4!}\theta^4-\frac{1}{6!}\theta^6+...\right)\\ &\qquad +i\left(\theta - \frac{1}{3!}\theta^3 + \frac{1}{5!}\theta^5-\frac{1}{7!}\theta^7+... \right)\\ &= \cos \theta + i\sin\theta \end{aligned} }

参考文献

  1. 詳細は省きますが、これは近似式ではなく、厳密かつ一意に成り立ちます。↩

  2. 例えばscipy.signal.periodogramやscipy.signal.welch↩

  3. 物理分野では時空間ドメインと周波数ドメインを対称に扱うためにPattern 2を用いる人が多いようです。↩