np.fft.fftとnp.fft.rfft
Category: python
概要
Numpyを使ってFFTをするときにこれまで fft を使っていたのだが、少なくとも工学系は多くの場合で rfft で十分なのではないかというお話。
NumpyのFFT
PythonでFFTするときはこれまで np.fft.fft を使っていた。 たとえばsin波をFFTするときはこんな感じ。
import numpy as np
f = 100 # 周波数 [Hz]
fs = 1e3 # サンプリング周波数 [Hz]
T = 100e-3 # 作成する長さ [s]
N = 64 # FFT点数
# 時刻の刻みを作成
t = np.arange(int(T*fs)) / fs
# 各時刻のsinを計算
signal = np.sin(2*np.pi*f*t)
# FFT用の窓
win = np.hamming(N)
# signalの先頭のN点を取り出してFFT
fft_signal = np.fft.fft(win * signal[:N])
signalはsin波で実数なので、この結果はこんな感じになる。
>>> print(fft_signal)
[ 1.74438136e-01+0.00000000e+00j 1.75287567e-01-3.73272090e-03j
1.74864008e-01-4.65222105e-03j 1.54867450e-01+6.33587976e-03j
-6.65125615e-03+8.15314012e-02j -2.64383920e+00+1.13615380e+00j
1.42449339e+01-4.87427558e+00j -1.22519849e+01+3.64618279e+00j
1.42875862e+00-2.86537737e-01j 4.61355027e-02+4.43652316e-02j
-8.10435063e-02+6.14333092e-02j -9.78142067e-02+5.65090513e-02j
-9.49462631e-02+4.95953839e-02j -8.83469778e-02+4.33662615e-02j
-8.16619311e-02+3.80927286e-02j -7.57470422e-02+3.36457356e-02j
-7.07203111e-02+2.98558880e-02j -6.65017794e-02+2.65799365e-02j
-6.29706590e-02+2.37069165e-02j -6.00116496e-02+2.11525123e-02j
-5.75267151e-02+1.88524960e-02j -5.54359559e-02+1.67574133e-02j
-5.36753399e-02+1.48286456e-02j -5.21938770e-02+1.30355650e-02j
-5.09510518e-02+1.13534853e-02j -4.99147118e-02+9.76217367e-03j
-4.90594004e-02+8.24475811e-03j -4.83650717e-02+6.78691140e-03j
-4.78161104e-02+5.37623094e-03j -4.74005959e-02+4.00175686e-03j
-4.71097580e-02+2.65358782e-03j -4.69375890e-02+1.32256471e-03j
-4.68805834e-02+0.00000000e+00j -4.69375890e-02-1.32256471e-03j
-4.71097580e-02-2.65358782e-03j -4.74005959e-02-4.00175686e-03j
-4.78161104e-02-5.37623094e-03j -4.83650717e-02-6.78691140e-03j
-4.90594004e-02-8.24475811e-03j -4.99147118e-02-9.76217367e-03j
-5.09510518e-02-1.13534853e-02j -5.21938770e-02-1.30355650e-02j
-5.36753399e-02-1.48286456e-02j -5.54359559e-02-1.67574133e-02j
-5.75267151e-02-1.88524960e-02j -6.00116496e-02-2.11525123e-02j
-6.29706590e-02-2.37069165e-02j -6.65017794e-02-2.65799365e-02j
-7.07203111e-02-2.98558880e-02j -7.57470422e-02-3.36457356e-02j
-8.16619311e-02-3.80927286e-02j -8.83469778e-02-4.33662615e-02j
-9.49462631e-02-4.95953839e-02j -9.78142067e-02-5.65090513e-02j
-8.10435063e-02-6.14333092e-02j 4.61355027e-02-4.43652316e-02j
1.42875862e+00+2.86537737e-01j -1.22519849e+01-3.64618279e+00j
1.42449339e+01+4.87427558e+00j -2.64383920e+00-1.13615380e+00j
-6.65125615e-03-8.15314012e-02j 1.54867450e-01-6.33587976e-03j
1.74864008e-01+4.65222105e-03j 1.75287567e-01+3.73272090e-03j]
64点でFFTしたので、（0を含んで）正の周波数32点と負の周波数32点の結果が格納されている。 各要素がどの周波数と対応するのかは np.fft.fftfreq で作成できる。
>>> fft_freq = np.fft.fftfreq(N, 1/fs)
>>> print(fft_freq)
[ 0. 15.625 31.25 46.875 62.5 78.125 93.75 109.375
125. 140.625 156.25 171.875 187.5 203.125 218.75 234.375
250. 265.625 281.25 296.875 312.5 328.125 343.75 359.375
375. 390.625 406.25 421.875 437.5 453.125 468.75 484.375
-500. -484.375 -468.75 -453.125 -437.5 -421.875 -406.25 -390.625
-375. -359.375 -343.75 -328.125 -312.5 -296.875 -281.25 -265.625
-250. -234.375 -218.75 -203.125 -187.5 -171.875 -156.25 -140.625
-125. -109.375 -93.75 -78.125 -62.5 -46.875 -31.25 -15.625]
並べて表示してみる。 複素数になってしまうので見づらくなるが、左側が周波数、右側が各周波数のFFT要素である。
>>> print(np.c_[fft_freq, fft_signal])
[[ 0.00000000e+00+0.00000000e+00j 1.74438136e-01+0.00000000e+00j]
[ 1.56250000e+01+0.00000000e+00j 1.75287567e-01-3.73272090e-03j]
[ 3.12500000e+01+0.00000000e+00j 1.74864008e-01-4.65222105e-03j]
[ 4.68750000e+01+0.00000000e+00j 1.54867450e-01+6.33587976e-03j]
[ 6.25000000e+01+0.00000000e+00j -6.65125615e-03+8.15314012e-02j]
[ 7.81250000e+01+0.00000000e+00j -2.64383920e+00+1.13615380e+00j]
[ 9.37500000e+01+0.00000000e+00j 1.42449339e+01-4.87427558e+00j]
[ 1.09375000e+02+0.00000000e+00j -1.22519849e+01+3.64618279e+00j]
[ 1.25000000e+02+0.00000000e+00j 1.42875862e+00-2.86537737e-01j]
[ 1.40625000e+02+0.00000000e+00j 4.61355027e-02+4.43652316e-02j]
[ 1.56250000e+02+0.00000000e+00j -8.10435063e-02+6.14333092e-02j]
[ 1.71875000e+02+0.00000000e+00j -9.78142067e-02+5.65090513e-02j]
[ 1.87500000e+02+0.00000000e+00j -9.49462631e-02+4.95953839e-02j]
[ 2.03125000e+02+0.00000000e+00j -8.83469778e-02+4.33662615e-02j]
[ 2.18750000e+02+0.00000000e+00j -8.16619311e-02+3.80927286e-02j]
[ 2.34375000e+02+0.00000000e+00j -7.57470422e-02+3.36457356e-02j]
[ 2.50000000e+02+0.00000000e+00j -7.07203111e-02+2.98558880e-02j]
[ 2.65625000e+02+0.00000000e+00j -6.65017794e-02+2.65799365e-02j]
[ 2.81250000e+02+0.00000000e+00j -6.29706590e-02+2.37069165e-02j]
[ 2.96875000e+02+0.00000000e+00j -6.00116496e-02+2.11525123e-02j]
[ 3.12500000e+02+0.00000000e+00j -5.75267151e-02+1.88524960e-02j]
[ 3.28125000e+02+0.00000000e+00j -5.54359559e-02+1.67574133e-02j]
[ 3.43750000e+02+0.00000000e+00j -5.36753399e-02+1.48286456e-02j]
[ 3.59375000e+02+0.00000000e+00j -5.21938770e-02+1.30355650e-02j]
[ 3.75000000e+02+0.00000000e+00j -5.09510518e-02+1.13534853e-02j]
[ 3.90625000e+02+0.00000000e+00j -4.99147118e-02+9.76217367e-03j]
[ 4.06250000e+02+0.00000000e+00j -4.90594004e-02+8.24475811e-03j]
[ 4.21875000e+02+0.00000000e+00j -4.83650717e-02+6.78691140e-03j]
[ 4.37500000e+02+0.00000000e+00j -4.78161104e-02+5.37623094e-03j]
[ 4.53125000e+02+0.00000000e+00j -4.74005959e-02+4.00175686e-03j]
[ 4.68750000e+02+0.00000000e+00j -4.71097580e-02+2.65358782e-03j]
[ 4.84375000e+02+0.00000000e+00j -4.69375890e-02+1.32256471e-03j]
[-5.00000000e+02+0.00000000e+00j -4.68805834e-02+0.00000000e+00j]
[-4.84375000e+02+0.00000000e+00j -4.69375890e-02-1.32256471e-03j]
[-4.68750000e+02+0.00000000e+00j -4.71097580e-02-2.65358782e-03j]
[-4.53125000e+02+0.00000000e+00j -4.74005959e-02-4.00175686e-03j]
[-4.37500000e+02+0.00000000e+00j -4.78161104e-02-5.37623094e-03j]
[-4.21875000e+02+0.00000000e+00j -4.83650717e-02-6.78691140e-03j]
[-4.06250000e+02+0.00000000e+00j -4.90594004e-02-8.24475811e-03j]
[-3.90625000e+02+0.00000000e+00j -4.99147118e-02-9.76217367e-03j]
[-3.75000000e+02+0.00000000e+00j -5.09510518e-02-1.13534853e-02j]
[-3.59375000e+02+0.00000000e+00j -5.21938770e-02-1.30355650e-02j]
[-3.43750000e+02+0.00000000e+00j -5.36753399e-02-1.48286456e-02j]
[-3.28125000e+02+0.00000000e+00j -5.54359559e-02-1.67574133e-02j]
[-3.12500000e+02+0.00000000e+00j -5.75267151e-02-1.88524960e-02j]
[-2.96875000e+02+0.00000000e+00j -6.00116496e-02-2.11525123e-02j]
[-2.81250000e+02+0.00000000e+00j -6.29706590e-02-2.37069165e-02j]
[-2.65625000e+02+0.00000000e+00j -6.65017794e-02-2.65799365e-02j]
[-2.50000000e+02+0.00000000e+00j -7.07203111e-02-2.98558880e-02j]
[-2.34375000e+02+0.00000000e+00j -7.57470422e-02-3.36457356e-02j]
[-2.18750000e+02+0.00000000e+00j -8.16619311e-02-3.80927286e-02j]
[-2.03125000e+02+0.00000000e+00j -8.83469778e-02-4.33662615e-02j]
[-1.87500000e+02+0.00000000e+00j -9.49462631e-02-4.95953839e-02j]
[-1.71875000e+02+0.00000000e+00j -9.78142067e-02-5.65090513e-02j]
[-1.56250000e+02+0.00000000e+00j -8.10435063e-02-6.14333092e-02j]
[-1.40625000e+02+0.00000000e+00j 4.61355027e-02-4.43652316e-02j]
[-1.25000000e+02+0.00000000e+00j 1.42875862e+00+2.86537737e-01j]
[-1.09375000e+02+0.00000000e+00j -1.22519849e+01-3.64618279e+00j]
[-9.37500000e+01+0.00000000e+00j 1.42449339e+01+4.87427558e+00j]
[-7.81250000e+01+0.00000000e+00j -2.64383920e+00-1.13615380e+00j]
[-6.25000000e+01+0.00000000e+00j -6.65125615e-03-8.15314012e-02j]
[-4.68750000e+01+0.00000000e+00j 1.54867450e-01-6.33587976e-03j]
[-3.12500000e+01+0.00000000e+00j 1.74864008e-01+4.65222105e-03j]
[-1.56250000e+01+0.00000000e+00j 1.75287567e-01+3.73272090e-03j]]
fftとrfft
FFT結果で注目すべき点は、-500Hzの要素を境に対象な形になっている点である。 例えば484.375Hzと-484.375Hzの要素を取り出すと、
[ 4.84375000e+02+0.00000000e+00j -4.69375890e-02+1.32256471e-03j]
[-4.84375000e+02+0.00000000e+00j -4.69375890e-02-1.32256471e-03j]
であり、FFT要素は複素共役となっていることが分かる。 この関係は他の周波数でも成り立っている。
また、-500HzのFFT要素の虚部はゼロである。 計算されていない500HzのFFT要素は-500HzのFFT要素の複素共役となると考えられることから、虚部がゼロの場合は-500HzのFFT要素と等しい。
なので、0〜484.375Hzの要素と-500Hzの要素、すなわち最初の33（=N/2 + 1）個の要素さえあれば情報量としては十分である。
これは実数をFFTしたときに成り立つ話であり、工学系の多くの用途に当てはまる。 このため、実数のFFTに特化した np.fft.rfft が用意されている。
# signalの先頭のN点を取り出して実数FFT
rfft_signal = np.fft.rfft(win * signal[:N])
rfft_signal の長さは33で、0〜484.375Hzと-500HzのFFT要素、すなわち fft_signal[:33] が格納されている。
実際には同じ計算をするわけではないため計算誤差により要素は完全には一致しない。 しかしその差は極めて小さい。
>>> print(fft_signal[:33] - rfft_signal)
[ 0.00000000e+00+0.00000000e+00j -1.11022302e-16+7.25114413e-16j
7.77156117e-16-1.73472348e-16j 2.22044605e-16-4.11129464e-16j
0.00000000e+00+0.00000000e+00j -4.44089210e-16-4.44089210e-16j
-1.77635684e-15+0.00000000e+00j 1.77635684e-15+4.44089210e-16j
0.00000000e+00-5.55111512e-17j 0.00000000e+00-1.11022302e-16j
-1.11022302e-15+4.44089210e-16j -1.11022302e-16-2.77555756e-16j
0.00000000e+00+0.00000000e+00j 2.08166817e-16-1.38777878e-17j
1.52655666e-16-4.85722573e-16j 7.07767178e-16-6.17561557e-16j
0.00000000e+00+0.00000000e+00j -4.02455846e-16-4.37150316e-16j
0.00000000e+00+1.24900090e-16j -1.31838984e-16+2.25514052e-16j
0.00000000e+00+0.00000000e+00j 5.55111512e-17-1.11022302e-16j
-4.44089210e-16-4.44089210e-16j 4.44089210e-16+7.77156117e-16j
1.11022302e-16+5.55111512e-17j -8.88178420e-16+0.00000000e+00j
8.88178420e-16+8.88178420e-16j 0.00000000e+00-1.11022302e-16j
-6.93889390e-18+0.00000000e+00j -1.52655666e-16+5.20417043e-18j
9.02056208e-16-2.11636264e-16j -8.04911693e-16+4.51028104e-17j
0.00000000e+00+0.00000000e+00j]
実行時間は半分にはならないが、 np.fft.rfft の方が速い。
>>> %timeit np.fft.fft(win * signal[:N])
30.7 µs ± 354 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
>>> %timeit np.fft.rfft(win * signal[:N])
21 µs ± 155 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
逆FFT np.fft.irfft も用意されている。 33点のrfft結果を入力すれば64点の時間領域信号が得られる。
>>> print(np.fft.irfft(rfft_signal).shape)
(64,)
周波数との対応は np.fft.rfftfreq で得られる。 最後の要素は-500Hzではなく500Hzとなっている（どっちも同じだけど）。
>>> rfft_freq = np.fft.rfftfreq(N, 1/fs)
>>> print(rfft_freq)
[ 0. 15.625 31.25 46.875 62.5 78.125 93.75 109.375 125.
140.625 156.25 171.875 187.5 203.125 218.75 234.375 250. 265.625
281.25 296.875 312.5 328.125 343.75 359.375 375. 390.625 406.25
421.875 437.5 453.125 468.75 484.375 500. ]
おわりに
本記事では実数のFFTに特化した np.fft.rfft を紹介した。 N点でFFTしたときになぜ(N/2 + 1)点になるのかが自分には謎だったので、本記事で簡単にまとめた。 実数をFFTして結果の半分を捨てている人は活用すると良いと思う。