はじめに
時系列データを扱っていると、
位置データから速度を求めたい
距離の変化から変化速度を知りたい
蛍光強度の変化速度を見たい
画像トラッキング結果から移動速度を出したい
という場面がよくあります。
このとき自然に使うのが「微分」です。
連続的な数式であれば、
$$
v(t) = \frac{dx(t)}{dt}
$$
と書けます。
ここで、$x(t)$ は時刻 $t$ における位置、$v(t)$ は速度です。
しかし、実験データや観測データでは、$x(t)$ を完全に正確に測れるわけではありません。多くの場合、データには少しノイズが入っています。
そして、この「少しのノイズ」が、微分するとかなり目立つことがあります。
この記事では、ダミーデータを使って、
位置データでは小さく見えるノイズが、速度にすると暴れる理由
をPythonで確認します。
Google Colabでそのまま実行できるコードを使います。
この記事で分かること
この記事では、次のことを扱います。
- 位置、速度、微分、差分、ノイズの基本
- なぜ微分するとノイズが増幅されるのか
- ダミーデータで位置と速度を比較する方法
- 平滑化してから微分するとどうなるか
- 「とりあえず微分する」前に考えるべきこと
先に結論
先に結論を書くと、重要なのは次の3点です。
- 実データでは、微分の代わりに「差分」を使うことが多い
- 差分では、隣り合うデータ点の差を $\Delta t$ で割る
- ノイズも同じように差分されるため、特に高周波のノイズが強調される
差分を式で書くと、例えば次のようになります。
$$
v_i \approx \frac{x_{i+1} - x_i}{\Delta t}
$$
ここで、$\Delta t$ はサンプリング間隔です。
もし位置データ $x_i$ にノイズが入っていると、そのノイズも引き算されて、さらに $\Delta t$ で割られます。
つまり、
微分・差分は、信号の変化を見ていると同時に、ノイズの変化も見てしまう
ということです。
用語の整理
位置とは何か
ここでは、時刻 $t$ におけるある量を $x(t)$ と書きます。
例えば、
- 物体の位置
- 細胞の輪郭位置
- 画像中の特徴点の座標
- センサーで測った変位
- ある測定量の時間変化
などをイメージできます。
この記事では、単位を具体的には決めず、a.u. とします。
a.u. は arbitrary unit の略で、「任意単位」という意味です。
速度とは何か
速度は、位置が時間あたりどれくらい変化したかを表す量です。
直感的には、
$$
速度 = \frac{位置の変化}{時間の変化}
$$
です。
連続的な関数なら、速度は微分で書けます。
$$
v(t) = \frac{dx(t)}{dt}
$$
ただし、実験データは連続的な関数ではなく、離散的な点列として得られます。
つまり、
$$
x_0, x_1, x_2, \cdots
$$
のようなデータです。
そのため、実際には微分そのものではなく、差分を使うことが多いです。
差分とは何か
差分とは、隣り合うデータ点の差を取ることです。
例えば、時刻 $t_i$ における位置を $x_i$ とします。
サンプリング間隔を $\Delta t$ とすると、速度は次のように近似できます。
$$
v_i \approx \frac{x_{i+1} - x_i}{\Delta t}
$$
これは前進差分と呼ばれる方法です。
「次の点との差」を使って変化率を求めています。
ノイズとは何か
ノイズとは、本来見たい信号に混ざっている不要な揺らぎです。
観測された位置を $x_i^{obs}$、本当の位置を $x_i^{true}$、ノイズを $\varepsilon_i$ とすると、
$$
x_i^{obs} = x_i^{true} + \varepsilon_i
$$
と書けます。
ここで、$\varepsilon_i$ が小さければ、位置データだけを見たときには「まあ問題なさそう」に見えるかもしれません。
しかし、速度を求めると話が変わります。
なぜ微分するとノイズが暴れるのか
観測された位置データから差分で速度を求めるとします。
\hat{v}_i=\frac{x_{i+1}^{obs} - x_i^{obs}}{\Delta t}
ここで、
$$
x_i^{obs} = x_i^{true} + \varepsilon_i
$$
を代入すると、
\hat{v}_i=\frac{(x_{i+1}^{true}+\varepsilon_{i+1}) - (x_i^{true}+\varepsilon_i)}{\Delta t}
となります。
整理すると、
\hat{v}_i=\frac{x_{i+1}^{true} - x_i^{true}}{\Delta t}+\frac{\varepsilon_{i+1} - \varepsilon_i}{\Delta t}
前半は、本来知りたい速度です。
$$
\frac{x_{i+1}^{true} - x_i^{true}}{\Delta t}
$$
後半は、ノイズ由来の速度です。
$$
\frac{\varepsilon_{i+1} - \varepsilon_i}{\Delta t}
$$
ここが重要です。
ノイズは単に残るだけではありません。
隣り合うノイズの差を取り、さらに $\Delta t$ で割る
ことになります。
もし各時刻のノイズが独立で、標準偏差が $\sigma_x$ だとすると、前進差分で生じる速度ノイズの標準偏差はおおよそ次のようになります。
$$
\sigma_v \approx \frac{\sqrt{2}\sigma_x}{\Delta t}
$$
つまり、$\Delta t$ が小さいほど、同じ位置ノイズでも速度ノイズは大きく見えます。
これは、
サンプリング周波数を高くすれば常に良い
という単純な話ではありません。
高サンプリング自体が悪いのではなく、微分前のノイズ処理、信号の時間スケール、フィルタ設計をセットで考える必要があるということです。
周波数の視点から見る
もう一つの見方として、微分は高周波成分を強調します。
例えば、
$$
x(t) = \sin(2\pi f t)
$$
を微分すると、
$$
\frac{dx(t)}{dt}=
2\pi f \cos(2\pi f t)
$$
となります。
右辺に $2\pi f$ が出てきます。
つまり、周波数 $f$ が高いほど、微分後の振幅が大きくなります。
ノイズには高周波成分が含まれることが多いため、微分によってノイズが目立ちやすくなります。
Google Colabで実験する
以下のコードは、Google Colabでそのまま実行できます。
今回は、次の流れで実験します。
- なめらかな「真の位置」を作る
- そこに小さなノイズを足して「観測された位置」を作る
- 差分で速度を求める
- ノイズが速度でどれくらい増えるかを見る
- 平滑化してから差分を取るとどうなるかを見る
# ============================================
# Differentiation amplifies noise
# Dummy position data -> velocity by finite difference
# ============================================
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
from IPython.display import display
# -----------------------------
# 1. Settings
# -----------------------------
SEED = 42
rng = np.random.default_rng(SEED)
os.makedirs("fig", exist_ok=True)
fs = 50 # sampling frequency [Hz]
dt = 1 / fs # sampling interval [s]
duration = 10 # duration [s]
t = np.arange(0, duration, dt)
# -----------------------------
# 2. Generate a smooth true position
# -----------------------------
# The true position is a mixture of two sine waves.
# This is just a toy model for educational purposes.
A1, f1 = 1.0, 0.4
A2, f2 = 0.25, 1.2
x_true = (
A1 * np.sin(2 * np.pi * f1 * t)
+ A2 * np.sin(2 * np.pi * f2 * t)
)
# Analytical true velocity
def true_velocity(tt):
return (
A1 * 2 * np.pi * f1 * np.cos(2 * np.pi * f1 * tt)
+ A2 * 2 * np.pi * f2 * np.cos(2 * np.pi * f2 * tt)
)
# -----------------------------
# 3. Add small measurement noise
# -----------------------------
sigma_x = 0.03
x_obs = x_true + rng.normal(0, sigma_x, size=t.size)
# -----------------------------
# 4. Finite difference
# -----------------------------
def forward_diff(y, dt):
return np.diff(y) / dt
# The velocity estimated by forward difference is located
# between t[i] and t[i+1], so we use the midpoint time.
t_v = t[:-1] + dt / 2
v_true = true_velocity(t_v)
v_clean_diff = forward_diff(x_true, dt)
v_noisy_diff = forward_diff(x_obs, dt)
# -----------------------------
# 5. Smoothing before differentiation
# -----------------------------
# Savitzky-Golay filter locally fits a polynomial.
# Here we use it as a simple smoothing method.
window_short = 15 # 15 points = 0.30 s when fs = 50 Hz
window_long = 51 # 51 points = 1.02 s when fs = 50 Hz
polyorder = 3
x_smooth_short = savgol_filter(
x_obs,
window_length=window_short,
polyorder=polyorder
)
x_smooth_long = savgol_filter(
x_obs,
window_length=window_long,
polyorder=polyorder
)
v_smooth_short = forward_diff(x_smooth_short, dt)
v_smooth_long = forward_diff(x_smooth_long, dt)
# -----------------------------
# 6. Quantify errors
# -----------------------------
def rmse(y, y_ref):
return np.sqrt(np.mean((y - y_ref) ** 2))
summary = pd.DataFrame({
"method": [
"clean position -> diff",
"noisy position -> diff",
"Savitzky-Golay 15 points -> diff",
"Savitzky-Golay 51 points -> diff",
],
"velocity_RMSE": [
rmse(v_clean_diff, v_true),
rmse(v_noisy_diff, v_true),
rmse(v_smooth_short, v_true),
rmse(v_smooth_long, v_true),
]
})
summary["velocity_RMSE"] = summary["velocity_RMSE"].round(4)
theory_velocity_noise_std = np.sqrt(2) * sigma_x / dt
actual_velocity_noise_std = np.std(forward_diff(x_obs - x_true, dt))
print(f"dt = {dt:.3f} s")
print(f"position noise std = {np.std(x_obs - x_true):.4f}")
print(f"theoretical velocity noise std ≈ {theory_velocity_noise_std:.3f}")
print(f"actual velocity noise std ≈ {actual_velocity_noise_std:.3f}")
display(summary)
# -----------------------------
# 7. Figure 1: Position data
# -----------------------------
plt.figure(figsize=(10, 4))
plt.plot(t, x_true, label="true position")
plt.plot(t, x_obs, alpha=0.65, label="observed position with noise")
plt.title("Position data: small noise is visible but not disastrous")
plt.xlabel("time [s]")
plt.ylabel("position [a.u.]")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig1_position_noise.png", dpi=160)
plt.show()
# -----------------------------
# 8. Figure 2: Velocity from finite difference
# -----------------------------
plt.figure(figsize=(10, 4))
plt.plot(t_v, v_true, label="true velocity")
plt.plot(t_v, v_clean_diff, linestyle="--", label="diff of clean position")
plt.plot(t_v, v_noisy_diff, alpha=0.65, label="diff of noisy position")
plt.title("Velocity from finite difference: noise is strongly amplified")
plt.xlabel("time [s]")
plt.ylabel("velocity [a.u./s]")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig2_velocity_noise.png", dpi=160)
plt.show()
# -----------------------------
# 9. Figure 3: Smoothing before differentiation
# -----------------------------
plt.figure(figsize=(10, 4))
plt.plot(t_v, v_true, label="true velocity")
plt.plot(t_v, v_noisy_diff, alpha=0.35, label="no smoothing")
plt.plot(t_v, v_smooth_short, label="Savitzky-Golay 15 points")
plt.plot(t_v, v_smooth_long, label="Savitzky-Golay 51 points")
plt.title("Smoothing before differentiation: noise reduction and distortion")
plt.xlabel("time [s]")
plt.ylabel("velocity [a.u./s]")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig3_smoothing_velocity.png", dpi=160)
plt.show()
# -----------------------------
# 10. Figure 4: RMSE comparison
# -----------------------------
plt.figure(figsize=(8, 4))
plt.bar(summary["method"], summary["velocity_RMSE"])
plt.title("Velocity RMSE for each method")
plt.ylabel("RMSE")
plt.xticks(rotation=25, ha="right")
plt.tight_layout()
plt.savefig("fig/fig4_rmse_bar.png", dpi=160)
plt.show()
# -----------------------------
# 11. Figure 5: Noise amplification by finite difference
# -----------------------------
fs_candidates = np.array([10, 20, 50, 100, 200, 500])
dt_candidates = 1 / fs_candidates
std_candidates = np.sqrt(2) * sigma_x / dt_candidates
plt.figure(figsize=(8, 4))
plt.plot(fs_candidates, std_candidates, marker="o")
plt.title("Theoretical velocity noise grows with sampling frequency")
plt.xlabel("sampling frequency [Hz]")
plt.ylabel("std of differentiated noise [a.u./s]")
plt.tight_layout()
plt.savefig("fig/fig5_noise_amplification.png", dpi=160)
plt.show()
実行結果
次のような値になります。
dt = 0.020 s
position noise std = 0.0288
theoretical velocity noise std ≈ 2.121
actual velocity noise std ≈ 1.930
速度推定のRMSEは、次のようになります。
method velocity_RMSE
clean position -> diff 0.0013
noisy position -> diff 1.9301
Savitzky-Golay 15 points -> diff 0.2627
Savitzky-Golay 51 points -> diff 0.5921
ここでRMSEは、真の速度と推定された速度のずれの大きさを表します。
RMSEは root mean squared error の略で、日本語では二乗平均平方根誤差と呼ばれます。
式で書くと、
$$
RMSE =
\sqrt{
\frac{1}{N}
\sum_{i=1}^{N}
(\hat{v}_i - v_i)^2
}
$$
です。
ここで、$\hat{v}_i$ は推定された速度、$v_i$ は真の速度です。
結果1：位置データではノイズは小さく見える
最初の図では、青線が真の位置、オレンジ色の線がノイズを含む観測位置です。
位置データだけを見ると、ノイズはそこまで大きく見えません。
確かに少しギザギザしていますが、全体の動きはかなりよく追えています。
この段階では、
このくらいなら速度もそれなりに求められそう
と思うかもしれません。
しかし、実際に差分を取ると状況が変わります。
結果2：速度にするとノイズが大きく見える
2枚目の図では、速度を比較しています。
-
true velocity：真の速度 -
diff of clean position：ノイズなし位置から差分で求めた速度 -
diff of noisy position：ノイズあり位置から差分で求めた速度
ノイズなし位置から求めた速度は、真の速度とほぼ重なります。
つまり、差分そのものが必ず悪いわけではありません。
問題は、ノイズあり位置から差分を取った場合です。
位置データでは小さく見えたノイズが、速度にすると大きなギザギザとして現れます。
これは、
$$
\frac{\varepsilon_{i+1} - \varepsilon_i}{\Delta t}
$$
という項が出てくるためです。
今回の設定では、位置ノイズの標準偏差は約 $0.03$ です。
一方、サンプリング間隔は、
$$
\Delta t = 0.02 \ \mathrm{s}
$$
です。
そのため、速度ノイズの標準偏差は、おおよそ
$$
\frac{\sqrt{2}\times 0.03}{0.02}
\approx 2.12
$$
になります。
位置ノイズとしては小さく見えても、速度にすると無視しにくい大きさになります。
結果3：平滑化してから微分すると改善するが、やりすぎには注意
3枚目の図では、ノイズを含む位置データに平滑化をかけてから差分を取っています。
ここでは、Savitzky-Golayフィルタを使いました。
Savitzky-Golayフィルタは、局所的に多項式を当てはめることで、波形の形をある程度保ちながら滑らかにする方法です。
今回の結果では、
- 平滑化なし：かなりギザギザ
- 15点の平滑化：真の速度にかなり近づく
- 51点の平滑化：滑らかにはなるが、細かい動きが丸まる
という結果になりました。
ここで大事なのは、
平滑化すれば必ず正しくなるわけではない
ということです。
平滑化はノイズを減らしますが、同時に本来の信号成分も削ることがあります。
今回のダミーデータには、$1.2$ Hzの成分が含まれています。
$1.2$ Hzの周期はおよそ、
$$
\frac{1}{1.2} \approx 0.83 \ \mathrm{s}
$$
です。
一方、51点の平滑化窓は、サンプリング周波数 $50$ Hz なので、
$$
\frac{51}{50} = 1.02 \ \mathrm{s}
$$
です。
これは、信号の周期よりも長い窓です。
そのため、細かい変化まで一緒に滑らかにしてしまいます。
結果4：RMSEで比較する
RMSEで見ると、結果はさらに分かりやすくなります。
今回の例では、
- ノイズなし位置から差分：RMSEはほぼゼロ
- ノイズあり位置から差分：RMSEが大きい
- 15点で平滑化してから差分：大きく改善
- 51点で平滑化してから差分：改善はするが、15点より悪い
という結果でした。
つまり、
微分前に平滑化することは有効だが、平滑化の強さは信号の時間スケールに合わせる必要がある
ということです。
結果5：同じ位置ノイズでも、サンプリング間隔で速度ノイズは変わる
最後の図では、理論式
$$
\sigma_v \approx \frac{\sqrt{2}\sigma_x}{\Delta t}
$$
に基づいて、サンプリング周波数と速度ノイズの関係を描いています。
サンプリング周波数を $f_s$ とすると、
$$
\Delta t = \frac{1}{f_s}
$$
です。
したがって、
$$
\sigma_v \approx \sqrt{2}\sigma_x f_s
$$
となります。
この式だけを見ると、サンプリング周波数が高いほど、差分で見える速度ノイズは大きくなります。
ただし、これは
各サンプルに同じ大きさの独立な位置ノイズが乗る
という単純化された条件での話です。
実際の測定では、センサーの性質、露光時間、フィルタ、信号帯域、ノイズの相関なども関係します。
そのため、この図から言いたいことは、
高サンプリングが悪い
ではありません。
そうではなく、
速度を求めるときは、サンプリング間隔とノイズの大きさを一緒に考える必要がある
ということです。
なぜ「微分してから考える」は危ないのか
微分は、変化を強調する操作です。
そのため、うまく使えば非常に強力です。
一方で、ノイズも「変化」として扱われます。
特に、ランダムにギザギザしたノイズは、隣り合う点で値が変わりやすいため、差分を取ると大きくなりやすいです。
その結果、
元データでは小さく見えたノイズが、微分後には大きな変動に見える
ということが起きます。
これは、位置から速度を求める場合だけではありません。
例えば、
- 速度から加速度を求める
- 蛍光強度から立ち上がり速度を求める
- 濃度変化から反応速度を求める
- 画像トラッキングから移動速度を求める
- センサー波形から変化率を求める
といった場面でも同じ問題が出ます。
実データで考えるべきこと
実データから速度や変化率を求める前には、少なくとも次の点を確認した方がよいです。
1. まず元データを見る
いきなり微分するのではなく、まず元の時系列を見ます。
元データに、
- 外れ値
- 欠損
- ドリフト
- 急なジャンプ
- 明らかな測定ミス
があると、微分後に大きなアーティファクトになります。
2. サンプリング間隔を確認する
速度の単位は、
$$
\frac{位置の単位}{時間の単位}
$$
です。
例えば、位置が $\mu m$、時間が秒なら、速度は $\mu m/s$ になります。
$\Delta t$ を間違えると、速度のスケールそのものが間違います。
3. 信号の時間スケールを考える
平滑化の窓幅は、信号の時間スケールより長すぎると、見たい変化まで消してしまいます。
例えば、0.5秒くらいで起きる変化を見たいのに、2秒の移動平均をかけると、その変化はかなり丸まってしまいます。
4. 平滑化のパラメータを変えて確認する
平滑化の窓幅を1つだけ選んで終わりにするのではなく、複数の値で結果を比較するのがおすすめです。
今回のように、
- 平滑化なし
- 弱い平滑化
- 強い平滑化
を並べると、解析結果がどれくらい処理条件に依存しているかが見えます。
5. できれば既知のダミーデータで試す
実データでは、真の速度は分かりません。
そのため、今回のように真の答えが分かっているダミーデータを作り、解析手法のクセを見ることは有用です。
これは、研究データ解析でもかなり重要です。
まとめ
この記事では、位置データから速度を求めるときに、なぜノイズが増幅されるのかを見ました。
重要なポイントは次の通りです。
- 実データでは、微分の代わりに差分を使うことが多い
- 差分は、隣り合う値の差を $\Delta t$ で割る
- ノイズも差分されるため、速度ではノイズが目立ちやすい
- 微分は高周波成分を強調する
- 平滑化してから微分すると改善することがある
- ただし、平滑化しすぎると本来の信号も削ってしまう
- 速度を求める前には、ノイズ、サンプリング間隔、信号の時間スケールを考える必要がある
一言でまとめると、
微分は便利だが、ノイズにも敏感な操作である
ということです。
おまけ：自分で試すと理解が深まる変更点
以下を変えて実行すると、挙動がかなり分かりやすくなります。
1. 位置ノイズを変える
sigma_x = 0.01
sigma_x = 0.05
sigma_x = 0.10
ノイズを大きくすると、速度がどれくらい暴れるかを確認できます。
2. サンプリング周波数を変える
fs = 20
fs = 100
fs = 200
同じ位置ノイズでも、$\Delta t$ が変わると速度ノイズの見え方が変わります。
3. 平滑化窓を変える
window_short = 7
window_short = 15
window_short = 31
窓を大きくするとノイズは減りますが、波形の細かい変化も失われやすくなります。
4. 真の信号の周波数を変える
f2 = 0.8
f2 = 2.0
f2 = 4.0
速い変化を含む信号では、強い平滑化がより危険になることがあります。
最後に
微分は、時系列データ解析の基本操作です。
しかし、実データに対しては、
微分する前に、ノイズを見る
微分する前に、サンプリング間隔を見る
微分する前に、信号の時間スケールを見る
という姿勢が大切です。
位置から速度を求めるだけなら簡単そうに見えます。
でも、その一歩手前にある「測定データをどう扱うか」が、解析結果の信頼性を大きく左右します。