雑記： numpy.fft の話 - クッキーの日記

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。全体的に離散フーリエ変換の備忘メモであり、高速フーリエ変換の原理の話は全くありません。

参考文献

Discrete Fourier Transform (numpy.fft) — NumPy v1.17 Manual

まとめ

numpy.fft.fft は、長さの入力系列を渡すと、これを式のように分解したときの係数、つまり式の右辺の各列ベクトルの成分の和で表したときの係数を返してくれる。それが実は式で計算できるので実際は式をやってくれている（高速なアルゴリズムで）。
- 入力系列 $\{a_m\}_{m=0}^{n-1}$ を numpy.fft.fft に渡して出た出力系列 $\{A_k\}_{k=0}^{n-1}$ を numpy.fft.ifft に渡すと元の入力系列 $\{a_m\}_{m=0}^{n-1}$ に戻る。単に $(4)$ 式をやってくれているだけである。元に戻ることは計算すればわかる。
- 入力系列 $\{a_m\}_{m=0}^{n-1}$ を numpy.fft.fft に渡して出た出力系列 $\{A_k\}_{k=0}^{n-1}$ は $(4)$ 式の右辺の各列ベクトルの成分の和で表したときの係数なので、周波数 $0, 1/n, \cdots, (n-1)/n$ の波の振幅と解釈できる。ただし、整数でしか値を取らない波なので、周波数には整数を足し引きしても値は変わらない。なので、 $n/2$ 以上の周波数成分からは $1$ を差し引いて周波数の取りうる範囲を $-n/2, \cdots, -1,0,1, \cdots, n/2-1$ （これは $n$ が偶数のとき）と $0$ を中心にしてもよい。こちらの流儀で numpy.fft.fft の出力系列がそれぞれどんな周波数の波の振幅と解釈できるかを取得できるのが numpy.fft.fftfreq である。入力系列のサイズ $n$ を渡すだけである。
- 入力系列 $\{a_m\}_{m=0}^{n-1}$ が時間 $\Delta t$ ごとにサンプリングされたものであるとき、これを反映した周波数にしたいならば、 $\exp(2 \pi i m k/n) = \exp(2 \pi i m \Delta t k/(n \Delta t ))$ より、 $m$ ごとの波を $m \Delta t$ ごとの波にする代わりに周波数 $k/n$ を $k/(n \Delta t)$ にすればよい。numpy.fft.fftfreq に $\Delta t$ も一緒に渡せばそれも踏まえて $\{A_k\}_{k=0}^{n-1}$ はそれぞれどんな周波数の波の振幅なのかを返してくれる。もっとも、numpy.fft.fft や numpy.fft.ifft は $(3)$ 式と $(4)$ 式をやるだけなので $\Delta t$ を気にすることはない。 $\{a_m\}_{m=0}^{n-1}$ や $\{A_k\}_{k=0}^{n-1}$ をグラフにプロットしたときの横軸の目盛にのみ $\Delta t$ が関わってくる。

f:id:cookie-box:20200101101603p:plain:w60

Discrete Fourier Transform (numpy.fft) — NumPy v1.17 Manual によると、numpy.fft における離散フーリエ変換と逆離散フーリエ変換の定義は以下であるとありますね。

$\displaystyle A_k = \sum_{m=0}^{n-1} a_m \exp \left\{ - 2 \pi i \frac{mk}{n} \right\} \qquad k = 0, \cdots, n-1 \tag{1}$ $\displaystyle a_m = \frac{1}{n}\sum_{k=0}^{n-1} A_k \exp \left\{2 \pi i \frac{mk}{n} \right\} \qquad m = 0, \cdots, n-1 \tag{2}$ 行列でかくと以下でしょうか。 $\left( \begin{array}{c} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \end{array} \right) = \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 1 & e^{-2 \pi i 1 \cdot 1 / n} & \cdots & e^{-2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (n-1)/ n} & \cdots & e^{-2 \pi i (n-1) \cdot (n-1)/ n} \end{array} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) \tag{3}$ $\left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) = \displaystyle \frac{1}{n} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 1 & e^{2 \pi i 1 \cdot 1 / n} & \cdots & e^{2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{2 \pi i 1 \cdot (n-1)/ n} & \cdots & e^{2 \pi i (n-1) \cdot (n-1)/ n} \end{array} \right) \left( \begin{array}{c} A_0 \\ A_1 \\ \vdots \\ A_{n-1} \end{array} \right) \tag{4}$ これらは、元の信号をサンプリングした系列 $\{a_m\}_{m=0}^{n-1}$ を周波数成分系列 $\{A_k\}_{k=0}^{n-1}$ にする式と、その逆向きの式ということですよね…これって $(1)$ で変換して $(2)$ で逆変換したら元に戻ってくるんですか？

※ $\exp(\pm 2 \pi i mk /n)$ は $k$ を $k-n$ としても値は変わらないので、 $k$ の動く範囲を $k = 0, \cdots, n-1$ とする代わりに、 $n/2$ 以上については $n$ を差し引いて $k = -n/2, \cdots, -1,0,1, \cdots, n/2-1$ としてもよい（但しこれは $n$ が偶数のとき； $n$ が奇数の場合は $k = -(n-1)/2, \cdots, (n-1)/2$ とできる）。このとき $(1), (2), (3), (4)$ は以下のようになる。 $A_{-n/2}$ はナイキスト周波数（元の信号をサンプリングした周波数の半分の周波数）の周波数成分に対応する振幅である。 $\displaystyle A_k = \sum_{m=0}^{n-1} a_m \exp \left\{ - 2 \pi i \frac{mk}{n} \right\} \qquad k = -n/2, \cdots, -1,0,1, \cdots, n/2-1 \tag{1'}$ $\displaystyle a_m = \frac{1}{n}\sum_{k=-n/2}^{n/2-1} A_k \exp \left\{2 \pi i \frac{mk}{n} \right\} \qquad m = 0, \cdots, n-1 \tag{2'}$ $\left( \begin{array}{c} A_{-n/2} \\ \vdots \\ A_{-1} \\ A_0 \\ A_1 \\ \vdots \\ A_{n/2-1} \end{array} \right) = \left( \begin{array}{cccc} 1 & e^{-2 \pi i 1 \cdot (-n/2) / n} & \cdots & e^{-2 \pi i (n-1) \cdot (-n/2) / n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (-1) / n} & \cdots & e^{-2 \pi i (n-1) \cdot (-1)/ n} \\ 1 & 1 & \cdots & 1 \\ 1 & e^{-2 \pi i 1 \cdot 1 / n} & \cdots & e^{-2 \pi i (n-1) \cdot 1/ n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-2 \pi i 1 \cdot (n/2-1)/ n} & \cdots & e^{-2 \pi i (n-1) \cdot (n/2-1)/ n} \end{array} \right) \left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) \tag{3'}$ $\left( \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{array} \right) = \displaystyle \frac{1}{n} \left( \begin{array}{ccccccc} 1 & \cdots & 1 & 1 & 1 & \cdots & 1 \\ e^{2 \pi i 1 \cdot (-n/2) / n} & \cdots & e^{2 \pi i 1 \cdot (-1) / n} & 1 & e^{2 \pi i 1 \cdot 1 / n} & \cdots & e^{2 \pi i 1 \cdot (n/2 - 1)/ n} \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ e^{2 \pi i (n-1) \cdot (-n/2) / n} & \cdots & e^{2 \pi i (n-1) \cdot (-1) / n} & 1 & e^{2 \pi i (n-1) \cdot 1 / n} & \cdots & e^{2 \pi i (n-1) \cdot (n/2-1)/ n} \end{array} \right) \left( \begin{array}{c} A_{-n/2} \\ \vdots \\ A_{-1} \\ A_0 \\ A_1 \\ \vdots \\ A_{n/2-1} \end{array} \right) \tag{4'}$

f:id:cookie-box:20200101101614p:plain:w60

$(1)$ を $(2)$ に代入したら戻ってくるね。

$\displaystyle \begin{split} a_m &= \frac{1}{n}\sum_{k=0}^{n-1} \left( \sum_{m'=0}^{n-1} a_{m'} \exp \left\{ - 2 \pi i \frac{m' k}{n} \right\} \right) \exp \left\{2 \pi i \frac{mk}{n} \right\} \\ &= \frac{1}{n} \sum_{k=0}^{n-1} \left( a_0 + a_1 e^{-2\pi i k/n} + a_2 e^{-2\pi i (2k)/n} + \cdots + a_{n-1} e^{-2\pi i (n-1)k/n} \right) e^{2 \pi i mk/n} \\ &= \frac{1}{n} \sum_{k=0}^{n-1} \left( a_0 e^{-2\pi i (-m)k/n} + a_1 e^{-2\pi i (1-m)k/n} + a_2 e^{-2\pi i (2-m)k/n} + \cdots + a_{n-1} e^{-2\pi i (n-1-m)k/n} \right) \\ &= \frac{a_0}{n} \frac{1 - e^{-2\pi i (-m)}}{1 - e^{-2\pi i (-m)/n}} + \frac{a_1}{n} \frac{1 - e^{-2\pi i (1-m)}}{1 - e^{-2\pi i (1-m)/n}} + \cdots + \frac{a_m}{n} n + \cdots + \frac{a_1}{n} \frac{1 - e^{-2\pi i (n-1-m)}}{1 - e^{-2\pi i (n-1-m)/n}} \\ &= a_m \end{split}$ というか、 $(4)$ の左辺のベクトル（入力系列）を、 $(4)$ の右辺の行列の各列の基底ベクトルで表現したときの係数が $A_0, \, A_1, \, \cdots, \, A_{n-1}$ だからね。 $(1)$ をして $(2)$ をしたら元に戻ったというより、 $(2)$ のような分解の係数を与えるのが $(1)$ なんだよね。

「元に戻ったというより元々そういう係数を出していたのだ」といわれても、 $(1)$ で係数が出てくること自体が非自明なんですがそれは…。まあ御託はいいです。とりあえず numpy.fft を使ってみましょう。とりあえず正弦波からサンプリングした系列を numpy.fft.fft して、再び numpy.fft.ifft してみればよいでしょう。

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 9, 3
rcParams['font.family']='Ume Hy Gothic O5'
rcParams['font.size'] = 15

T = 4
f = 1 / T
dt = 0.5

x = np.arange(-2, 2, dt)
y = np.sin(2 * np.pi * f * x)

X = np.fft.fftfreq(y.size, d=dt)
Y = np.fft.fft(y)

y_ = np.fft.ifft(Y)

print('周期   {}'.format(T))
print('周波数 {}'.format(f))
print('サンプル数             {}'.format(y.size))
print('サンプリング時間間隔   {}'.format(dt))
print('サンプリングレート     {}'.format(1/dt))

fig, ax = plt.subplots(5, 1, figsize=(7, 11))

ax[0].scatter(x, y)
ax[0].set_ylabel('元の信号')

ax[1].scatter(X, Y.real)
ax[1].set_ylabel('fftの実部')

ax[2].scatter(X, Y.imag)
ax[2].set_ylabel('fftの虚部')

ax[3].scatter(x, y_.real)
ax[3].set_ylabel('ifftの実部')

ax[4].scatter(x, y_.imag)
ax[4].set_ylabel('ifftの虚部')

plt.subplots_adjust(hspace=0.4)
plt.show()

周期   4
周波数 0.25
サンプル数             8
サンプリング時間間隔   0.5
サンプリングレート     2.0

f:id:cookie-box:20200128123638p:plain:w400

…確かに ifft の実部に元の正弦波が復元されたようにみえますね。みえますが、しかし、何が起きたのかよくわかりません…。numpy.fft さんは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの周波数成分の振幅を出してくださったようですが、そもそも私はそのようにお願いした覚えがなく…。

みづらいけど、-0.25 ヘルツの振幅が $-4i$ 、+0.25 ヘルツの振幅が $4i$ になっているね。確かに周期4の正弦波が復元されるね。

$\displaystyle \begin{split} & \frac{1}{8} \left( -4i e^{2 \pi i m (-0.25)} + 4i e^{2 \pi i m (0.25)}\right) \\ &= \frac{1}{8} \left( -4i \cos \left(-\frac{m \pi}{2} \right) + 4 \sin \left(-\frac{m \pi}{2} \right) + 4i \cos \left(\frac{m \pi}{2} \right) - 4 \sin \left(\frac{m \pi}{2} \right) \right) \\ &= - \sin \left(\frac{m \pi}{2} \right) \end{split}$ 部長が最初につくった正弦波と符号がずれているけど、部長は元の信号の左端を $-2$ にしているから、部長の座標の取り方に直せばこれを2だけマイナス側にずらして $\displaystyle - \sin \left(\frac{(m+2) \pi}{2} \right) = \sin \left(\frac{m \pi}{2} \right)$ だね。

ああ、numpy.fft さんは私が渡した系列の1つ目を時刻ゼロにおける信号の値と受けとっているのですか？ …ん？それでは、2つ目の値は何であると受け取っているのでしょう？ numpy.fft.fft にはサンプリング系列のみしか渡しておらず、どの時刻の値であるかなどは渡していません。

…numpy.fft.fft は、入力系列のベクトルを $(4)$ の右辺の行列の各列の成分ごとに分解するだけだよ。だから、1つ目の値はベクトルの第1成分、2つ目の値はベクトルの第2成分、…くらいにしか思ってないだろう。それで出てくる系列は $(4)$ の右辺の行列の各列ベクトルの成分の係数になる。各列ベクトルは周波数が $0, 1/n, \cdots, (n-1)/n$ である波の $m = 0, \cdots, n-1$ での値になっているから、出力系列は周波数 $0, 1/n, \cdots, (n-1)/n$ の周波数成分の振幅と解釈してもいい。

では、上の例は $n=8$ ですから、周波数 0, 1/8=0.125, …, 7/8=0.875 ヘルツの成分の振幅になって…いませんね。出ているのは -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツの成分の振幅ですよ？

2点あって、まず、 $m$ が整数しか取らない世界では $\exp(2 \pi i m k/n)$ の $k/n$ に整数を足し引きしても値は変わらない。だから、0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツの成分を、4/8, 5/8, 6/8, 7/8, 0, 1/8, 2/8, 3/8 ヘルツと並び替えて、さらに -4/8, -3/8, -2/8, -1/8, 0, 1/8, 2/8, 3/8 ヘルツだと考えていい。描画するならこう並べ替えた方が低周波数成分が真ん中に集まってみやすいだろう。加えてもう1点は、もし元の入力系列に対して「これは時間 $\Delta t$ ごとにサンプリングしたものだ」という思いを込めたいなら、 $\exp(2 \pi i m \Delta t k/(n \Delta t ))$ と、周波数側にも $\Delta t$ を押し付けることができる。部長は $\Delta t = 0.5$ にしたから、これを反映して周波数を1/0.5倍=2倍にする。つまり、-8/8, -6/8, -4/8, -2/8, 0, 2/8, 4/8, 6/8 ヘルツになる。-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75 ヘルツと一致するね。これまで単位時間に起きていたと思っていたことを「やっぱり時間 0.5 の間に起きていたことにして」としたんだから、周波数は2倍になるよね。でも、それらは人間がグラフの軸目盛にかく数字の問題でしかない。入力の点数が同じ8点なら、 $\Delta t$ がどうあったって基底ベクトルは同じだからね（記事の下部の2つのケース；横軸の目盛と何ヘルツと解釈するかが違うだけで、絵は全く同じ）。まとめると、numpy.fft.fft の返り値は素朴には 0, 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8 ヘルツに対応する振幅と受け取れるんだけど、この2点を反映すると何ヘルツに対応する振幅と受け取るべきかを取得できるのが、部長が上のスクリプトでよんでいる np.fft.fftfreq だね。

つづかない

numpy.fft.fft に8点の系列を渡したときの基底（左列が実部、右列が虚部）

各行が $(4')$ の各列ベクトルに対応。

x = np.linspace(0.0, 7.0, 8)
x_ = np.arange(0.0, 8.0, 0.01) # 点だけだとよくわからないのでcos, sinカーブを描画する用
vec_f = np.sort(np.fft.fftfreq(8, d=1)) # 8点の系列を渡したときの周波数成分

fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
    ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
    ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
    ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()

f:id:cookie-box:20200127231000p:plain:w720

numpy.fft.fft に8点の系列を渡してΔt=0.5と指定したときの基底（左列が実部、右列が虚部）

各行が $(4')$ の各列ベクトルに対応。

x = np.linspace(0.0, 3.5, 8)
x_ = np.arange(0.0, 4.0, 0.01) # 点だけだとよくわからないのでcos, sinカーブを描画する用
vec_f = np.sort(np.fft.fftfreq(8, d=0.5)) # 8点の系列を渡してかつΔt=0.5としたときの周波数成分

fig, ax = plt.subplots(8, 2, figsize=(11, 13))
for i, f in enumerate(vec_f):
    ax[i, 0].plot(x_, np.cos(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 1].plot(x_, np.sin(2 * np.pi * f * x_), linestyle='dotted')
    ax[i, 0].scatter(x, np.cos(2 * np.pi * f * x))
    ax[i, 1].scatter(x, np.sin(2 * np.pi * f * x))
    ax[i, 0].text(-0.5, 0.75, '{} Hz'.format(f), transform=ax[i, 0].transAxes)
plt.subplots_adjust(hspace=0.4)
plt.show()

f:id:cookie-box:20200128124655p:plain:w720