実数信号のFFTではrfftで片側の係数を求め、rfftfreqへ標本間隔を渡して周波数軸を作る。振幅は標本数で割り、正負の周波数をまとめる内部の成分だけ2倍する。DC成分と、標本数が偶数のときのNyquist成分を2倍しないことがポイントである。
最小例で確かめる
以下は説明用に作った小さなデータである。実測データや実行速度の測定結果ではない。コード全体をexample.pyとして保存すれば、入力ファイルを別途用意せずに実行できる。assertは、この例で成り立つべき形や値を確認するために入れてある。
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import rfft, rfftfreq, irfft
def amplitude_spectrum(signal, dt):
signal = np.asarray(signal, dtype=float)
if signal.ndim != 1 or signal.size < 2:
raise ValueError("one-dimensional signal with at least two samples")
if not np.isfinite(signal).all() or not np.isfinite(dt) or dt <= 0:
raise ValueError("finite signal and positive sample interval required")
n = signal.size
spectrum = rfft(signal) # 既定のbackward正規化
amplitude = np.abs(spectrum) / n
if n % 2 == 0:
amplitude[1:-1] *= 2
else:
amplitude[1:] *= 2
return rfftfreq(n, d=dt), amplitude, spectrum
fs = 64.
n = 128
t = np.arange(n) / fs
signal = (2. + 1.5 * np.sin(2*np.pi*8*t)
+ 0.5 * np.cos(2*np.pi*20*t) + 0.25 * np.cos(np.pi*np.arange(n)))
freq, amp, spectrum = amplitude_spectrum(signal, 1/fs)
for hz in [0., 8., 20., 32.]:
k = int(np.flatnonzero(np.isclose(freq, hz))[0])
print(f"{hz:4.1f} Hz: {amp[k]:.6f}")
np.testing.assert_allclose(amp[[0,16,40,64]], [2,1.5,0.5,0.25], atol=1e-12)
np.testing.assert_allclose(irfft(spectrum, n=n), signal, atol=1e-12)
# 奇数長では最後のbinにも負周波数側の相手がある。
odd_t = np.arange(15) / 15
odd_signal = 0.3 + np.cos(2*np.pi*3*odd_t) + 0.2*np.cos(2*np.pi*7*odd_t)
odd_freq, odd_amp, _ = amplitude_spectrum(odd_signal, 1/15)
np.testing.assert_allclose(odd_amp[[0,3,7]], [0.3,1.,0.2], atol=1e-12)
print("odd-length final bin:", round(float(odd_amp[-1]), 6))
fig, axes = plt.subplots(2, 1, figsize=(7, 6), layout="constrained")
axes[0].plot(t[:64], signal[:64], color="#26778c", marker=".")
axes[0].set(xlabel="Time [s]", ylabel="Signal [arb. unit]", title="Synthetic sampled signal (first second)")
axes[1].stem(freq, amp, basefmt=" ")
axes[1].set(xlabel="Frequency [Hz]", ylabel="One-sided amplitude", xlim=(-1,33),
title="DC and Nyquist bins are not doubled")
fig.savefig("PY087-fft.png", dpi=160)
plt.close(fig)
実行結果
0.0 Hz: 2.000000
8.0 Hz: 1.500000
20.0 Hz: 0.500000
32.0 Hz: 0.250000
odd-length final bin: 0.2
横軸は標本間隔から作る
rfftの出力の添字はそのままHzではない。rfftfreq(n, d=dt)へ標本数nと秒単位の標本間隔dtを渡すと、周波数がHzで得られる。例では64Hzで128点を測ったので、周波数の刻みは64/128=0.5Hzである。8Hzは添字16に現れる。
tはnp.arange(n)/fsとして等間隔に作った。区間の終点まで含めるlinspaceを無条件に使うと、想定した標本間隔と違う場合がある。実測ログでは時刻差のばらつきや欠測を確認し、不等間隔データをそのまま等間隔FFTとして解釈しないこと。ゼロ埋めで表示点を増やしても、元の観測時間に由来する分解能が増えるわけではない。
片側振幅の2倍を掛ける場所
既定の正規化では、FFT係数をnで割ると両側スペクトルの尺度になる。実数信号の正負周波数の対を片側へまとめるため、内部の正周波数成分は2倍する。一方、DCは一つだけで、偶数長のNyquist成分も正負の端が同じbinに対応するので2倍しない。
奇数長では最後のbinはNyquistではないため、最後も2倍する。コードでは偶数長と奇数長で切り出し方を分け、15点の別の信号でも検査した。また、この例のDCの値2は正の定数オフセットである。absを取った振幅からは、負のオフセットの符号までは読めない。
既知の信号で正規化を検証する
模擬信号には8Hzの振幅1.5、20Hzの振幅0.5、DCの2、Nyquistの振幅0.25を入れた。これらが期待するbinと大きさに現れることを検査している。Nyquist項は標本上の交互符号に相当し、この離散信号としての振幅を示す。任意位相の連続波を標本化したときの情報が全て残るという意味ではない。
irfftへ元のnを渡して逆変換し、入力信号と一致することも確認した。これは係数の往復の検査であり、振幅だけから元の信号を復元したわけではない。位相情報を捨てるabsと、複素係数を使う逆変換は区別しよう。normを変更する場合は、ここで用いた振幅式もそのまま流用しない。
漏れとaliasingを別々に考える
この例は各波の周期が観測窓にちょうど収まるように作っている。一般の信号では境界がつながらず、近くの周波数へ成分が広がるスペクトル漏れが起こる。窓関数は漏れの形を変えるが、振幅の尺度も変わるため、窓の利得補正などが必要になる。このコードは窓なしの整合した例に限定した正規化である。
標本周波数の半分を超える成分は折り返して観測され得る。FFTの後処理だけで元の高周波成分と区別することはできず、測定前の帯域制限が必要になる。また、この片側振幅はパワースペクトル密度ではない。縦軸の単位、窓、正規化、観測時間を併記して、異なるスペクトルを混同しないようにしたい。
動作確認環境と参考資料
Linux・CPython 3.12.14、NumPy 2.3.5、pandas 2.2.3、SciPy 1.17.0、Matplotlib 3.10.8の環境で掲載コードを実行した。使用するライブラリはコード冒頭のimportを参照してほしい。公式資料の最新版と、この実行確認版は区別している。数値の末尾や表の表示幅は環境によって変わることがある。