相補フィルタの理論とPython実装：ジャイロスコープと加速度計のセンサフュージョン

相補フィルタとは

相補フィルタ（Complementary Filter）は、特性の異なる2つのセンサの出力を周波数領域で最適に統合するフィルタです。一方のセンサにはハイパスフィルタを、もう一方にはローパスフィルタを適用し、両者の和が全帯域で1になるように設計されます。

最も典型的な応用はIMU（慣性計測装置）におけるジャイロスコープと加速度計のセンサフュージョンです。

ジャイロスコープは角速度を積分して角度を求めるため、バイアス誤差がドリフトとして蓄積します（低周波誤差）。一方、加速度計は重力ベクトルから傾斜角を算出できますが、振動や並進加速度の影響を受けやすい（高周波ノイズ）。

相補フィルタはこの特性の相補性を利用し、ジャイロからは高周波成分を、加速度計からは低周波成分を抽出して融合することで、両方の欠点を相殺します。

主な応用分野

• ドローン・マルチコプター: 姿勢推定（ピッチ・ロール角）
• ロボティクス: 二足歩行ロボットの姿勢安定化
• スマートフォン: 画面回転、AR/VRヘッドトラッキング
• 車両制御: 横滑り角の推定

理論

ハイパスフィルタとローパスフィルタの相補性

相補フィルタの核心は、ハイパスフィルタ \(H_{HPF}(z)\) とローパスフィルタ \(H_{LPF}(z)\) が以下の関係を満たすことです。

この性質により、2つのセンサ信号を統合した結果が全帯域にわたって歪みなく再構成されます。

1次相補フィルタの伝達関数

1次ローパスフィルタの伝達関数を次のように定義します。

相補条件 \((1)\) から、ハイパスフィルタは次のようになります。

相補フィルタの出力は、ジャイロ出力 \(X_{gyro}(z)\) と加速度計出力 \(X_{accel}(z)\) に対して次のように表されます。

時間領域での更新式

式 \((4)\) を逆Z変換して時間領域に戻すと、以下の再帰的な更新式が得られます。

ここで、

• \(\theta[n]\) : 現在の推定角度
• \(\dot{\theta}_{gyro}[n]\) : ジャイロスコープの角速度出力
• \(\theta_{accel}[n]\) : 加速度計から算出した角度
• \(\Delta t\) : サンプリング周期
• \(\alpha\) : フィルタ係数（\(0 < \alpha < 1\) ）

式 \((5)\) は直感的に次のように解釈できます。

• 第1項: ジャイロの積分値に重み \(\alpha\) を乗じる（ハイパスフィルタ効果）
• 第2項: 加速度計の角度に重み \((1-\alpha)\) を乗じる（ローパスフィルタ効果）

時定数 \(\tau\) とフィルタ係数 \(\alpha\) の関係

パラメータ \(\alpha\) は時定数 \(\tau\) とサンプリング周期 \(\Delta t\) から決定されます。

時定数 \(\tau\) の物理的意味は次のとおりです。

• \(\tau\) が大きい（\(\alpha \to 1\) ）: ジャイロの寄与が大きい。短期応答に優れるが、ドリフト補正が遅い
• \(\tau\) が小さい（\(\alpha \to 0\) ）: 加速度計の寄与が大きい。ドリフト補正は速いが、振動ノイズに弱い

クロスオーバー周波数 \(f_c\) （LPFとHPFのゲインが等しくなる周波数）は次の関係にあります。

定常カルマンフィルタとの等価性：厳密な導出

「1次の相補フィルタは定常カルマンフィルタの特殊ケースである」という主張を天下り的に述べるのではなく、スカラー（1次元）カルマンフィルタから出発して導出します。ジャイロ角速度 \(\dot\theta_{gyro}[n]\) を予測ステップの制御入力、加速度計角度 \(\theta_{accel}[n]\) を観測値とするスカラー状態空間モデルを考えます。

予測ステップ：

ここで \(Q\) はジャイロ雑音（標準偏差 \(\sigma_g\) ）が1ステップの積分で角度の不確かさに与える分散で、\(Q \approx (\sigma_g \Delta t)^2\) と近似できます。

更新ステップ：

ここで \(R = \sigma_a^2\) は加速度計の角度換算ノイズ分散です。式 \((9)\) の右辺で \(\alpha := 1 - K[n]\) と置けば、これは式 \((5)\) の更新式と完全に同じ形になります（\(\hat\theta^-[n]\) が第1項のジャイロ積分、\(\theta_{accel}[n]\) が第2項に対応）。つまり相補フィルタとカルマンフィルタは同一の更新式を共有しており、両者の違いはゲイン \(\alpha\) （あるいは \(K\) ）を固定するか時変にするかだけです。

定常ゲインの導出： \(Q, R\) が時不変であれば、共分散 \(P[n]\) は離散リカッチ方程式の不動点 \(P^*\) に収束します。定常状態 \(P[n]=P[n-1]=P^*\) を式 \((8)(9)\) に代入すると、

この2式から \(P^*\) を消去し、\(\alpha^* := 1 - K^*\) と雑音比 \(r := Q/R\) を用いて整理すると、次の2次方程式が得られます（導出：\(P^*=\alpha^*(P^*+Q)\) より \(P^*+Q=P^*/\alpha^*\) 。これを \(K^*\) の式に代入し \(\alpha^*=1-K^*\) を使うと \(P^*=(1-\alpha^*)R\) となり、両者を等置して \(r=Q/R\) で整理する）。

\(0 < \alpha^* < 1\) を満たす解を選ぶと、閉形式の解が得られます。

式 \((12)\) は、相補フィルタの固定パラメータ \(\alpha\) が、ジャイロ・加速度計の雑音分散比 \(r=Q/R\) にちょうど整合するときにのみ、定常カルマンフィルタと数学的に等価になることを示しています。逆に言えば、相補フィルタで \(\alpha\) （や \(\tau\) ）を手動チューニングする作業は、暗黙のうちにこの雑音比を推定していることに相当します。式 \((6)\) の \(\tau = \alpha\Delta t/(1-\alpha)\) と組み合わせれば、雑音統計量 \(\sigma_g, \sigma_a\) から理論的に最適な時定数 \(\tau^*\) を計算できます。この理論値を実データで検証した結果は後述の「パラメータ \(\alpha\) の影響」で数値的に確認します。

なお、この導出は \(Q\) をゼロ平均の過程ノイズとしてモデル化しており、ジャイロのバイアス（非ゼロ平均の系統誤差）は考慮していません。バイアスが存在する実センサでは式 \((12)\) の予測は理論上限としてのみ有効であり、実際の最適 \(\tau\) はバイアスの影響でこれより小さくなる方向にずれます。この点は後述のエッジケースで詳しく検証します。

周波数応答解析

相補フィルタの周波数応答を解析すると、LPFとHPFのゲインがクロスオーバー周波数 \(f_c\) で交差し、全帯域でゲインの合計が0 dBに保たれることを確認できます。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# --- パラメータ ---
fs = 100.0           # サンプリング周波数 [Hz]
dt = 1.0 / fs
tau = 1.0            # 時定数 [s]
alpha = tau / (tau + dt)

BLUE, ORANGE, RED = '#2a78d6', '#eb6834', '#e34948'  # カテゴリカルパレット

# --- 伝達関数の定義 ---
# LPF: H_LPF(z) = (1-alpha) / (1 - alpha*z^{-1})
b_lpf = [1 - alpha]
a_lpf = [1, -alpha]

# HPF: H_HPF(z) = (1 - z^{-1}) / (1 - alpha*z^{-1})
b_hpf = [1, -1]
a_hpf = [1, -alpha]

# --- 周波数応答 ---
w_lpf, h_lpf = signal.freqz(b_lpf, a_lpf, worN=4096, fs=fs)
w_hpf, h_hpf = signal.freqz(b_hpf, a_hpf, worN=4096, fs=fs)

# 合計の周波数応答
h_sum = h_lpf + h_hpf

fc = 1.0 / (2 * np.pi * tau)  # クロスオーバー周波数

# --- プロット ---
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# ゲイン特性
ax1.semilogx(w_lpf, 20 * np.log10(np.abs(h_lpf) + 1e-10),
             label='LPF (Accelerometer)', linewidth=2, color=BLUE)
ax1.semilogx(w_hpf, 20 * np.log10(np.abs(h_hpf) + 1e-10),
             label='HPF (Gyroscope)', linewidth=2, color=ORANGE)
ax1.semilogx(w_lpf, 20 * np.log10(np.abs(h_sum) + 1e-10),
             label='Sum (HPF + LPF)', linewidth=2, linestyle='--', color='black')
ax1.axvline(fc, color=RED, linestyle=':', alpha=0.8,
            label=f'$f_c$ = {fc:.3f} Hz')
ax1.axhline(-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
ax1.set_xlabel('Frequency [Hz]')
ax1.set_ylabel('Gain [dB]')
ax1.set_title(f'Complementary Filter Frequency Response (τ={tau} s, α={alpha:.4f})')
ax1.legend()
ax1.grid(True, which='both', alpha=0.3)
ax1.set_xlim([0.001, fs / 2])
ax1.set_ylim([-40, 5])

# 位相特性
ax2.semilogx(w_lpf, np.angle(h_lpf, deg=True), label='LPF', linewidth=2, color=BLUE)
ax2.semilogx(w_hpf, np.angle(h_hpf, deg=True), label='HPF', linewidth=2, color=ORANGE)
ax2.axvline(fc, color=RED, linestyle=':', alpha=0.8)
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Phase [degrees]')
ax2.set_title('Phase Response')
ax2.legend()
ax2.grid(True, which='both', alpha=0.3)
ax2.set_xlim([0.001, fs / 2])

plt.tight_layout()
plt.savefig('complementary_frequency_response.png', dpi=150, bbox_inches='tight')
plt.show()

実行結果は \(\tau = 1.0\) s のとき \(\alpha = 0.9901\) 、クロスオーバー周波数 \(f_c = 0.159\) Hz となり、LPFとHPFのゲインはこの周波数で \(-3\) dB に交差し、合計（黒破線）は全帯域で 0 dB（= 1）に保たれていることが確認できます。これが「相補（complementary）」の名前の由来です。位相特性を見ると、クロスオーバー周波数付近でLPFとHPFの位相が逆符号かつほぼ対称に振れており、振幅だけでなく位相の観点でも打ち消し合いが生じていることが分かります。

Pythonによる実装

IMUセンサのシミュレーション

実際のIMUセンサをシミュレーションし、相補フィルタの効果を検証します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)  # 再現性のため乱数シードを固定

BLUE, ORANGE, GREEN, RED = '#2a78d6', '#eb6834', '#1baf7a', '#e34948'

# --- シミュレーションパラメータ ---
fs = 100.0           # サンプリング周波数 [Hz]
dt = 1.0 / fs
T = 30.0             # シミュレーション時間 [s]
t = np.arange(0, T, dt)
N = len(t)

# --- 真の角度（テスト信号）---
# ゆっくりした正弦波運動 + ステップ変化
theta_true = (
    15.0 * np.sin(2 * np.pi * 0.1 * t) +     # 0.1 Hz の揺動
    5.0 * np.sin(2 * np.pi * 0.5 * t) +       # 0.5 Hz の揺動
    10.0 * (t > 15)                             # t=15s でステップ変化
)

# --- 真の角速度 ---
omega_true = np.gradient(theta_true, dt)

# --- ジャイロスコープ出力（角速度 + バイアスドリフト + ノイズ）---
gyro_bias = 0.5      # バイアス [deg/s]（ドリフトの原因）
gyro_noise_std = 2.0  # ノイズ標準偏差 [deg/s]

gyro_output = omega_true + gyro_bias + gyro_noise_std * np.random.randn(N)

# --- 加速度計出力（角度 + 振動ノイズ）---
accel_noise_std = 5.0  # ノイズ標準偏差 [deg]

accel_output = theta_true + accel_noise_std * np.random.randn(N)

# --- ジャイロのみによる角度推定（積分）---
theta_gyro = np.zeros(N)
theta_gyro[0] = theta_true[0]
for i in range(1, N):
    theta_gyro[i] = theta_gyro[i - 1] + gyro_output[i] * dt

# --- 相補フィルタ ---
tau = 1.0            # 時定数 [s]
alpha = tau / (tau + dt)

theta_comp = np.zeros(N)
theta_comp[0] = theta_true[0]
for i in range(1, N):
    # 式 (5) の実装
    theta_comp[i] = alpha * (theta_comp[i - 1] + gyro_output[i] * dt) \
                  + (1 - alpha) * accel_output[i]

# --- 比較プロット ---
fig, axes = plt.subplots(4, 1, figsize=(12, 14), sharex=True)

# 真の角度
axes[0].plot(t, theta_true, 'k-', linewidth=2, label='True angle')
axes[0].set_ylabel('Angle [deg]')
axes[0].set_title('True Angle')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ジャイロのみ
axes[1].plot(t, theta_gyro, color=ORANGE, alpha=0.9, label='Gyro only (integration)')
axes[1].plot(t, theta_true, 'k--', alpha=0.5, label='True')
axes[1].set_ylabel('Angle [deg]')
axes[1].set_title('Gyroscope Only (drift accumulation)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 加速度計のみ
axes[2].plot(t, accel_output, color=BLUE, alpha=0.6, label='Accelerometer only')
axes[2].plot(t, theta_true, 'k--', alpha=0.5, label='True')
axes[2].set_ylabel('Angle [deg]')
axes[2].set_title('Accelerometer Only (noisy)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

# 相補フィルタ
axes[3].plot(t, theta_comp, color=GREEN, linewidth=2, label=f'Complementary (α={alpha:.3f})')
axes[3].plot(t, theta_true, 'k--', alpha=0.5, label='True')
axes[3].set_ylabel('Angle [deg]')
axes[3].set_xlabel('Time [s]')
axes[3].set_title('Complementary Filter')
axes[3].legend()
axes[3].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('complementary_filter_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

# --- RMSE の計算 ---
rmse_gyro = np.sqrt(np.mean((theta_gyro - theta_true) ** 2))
rmse_accel = np.sqrt(np.mean((accel_output - theta_true) ** 2))
rmse_comp = np.sqrt(np.mean((theta_comp - theta_true) ** 2))

print(f"RMSE (Gyro only):          {rmse_gyro:.2f} deg")
print(f"RMSE (Accelerometer only): {rmse_accel:.2f} deg")
print(f"RMSE (Complementary):      {rmse_comp:.2f} deg")

seed=42 で実行すると、次の実測値が得られます。

RMSE (Gyro only):          9.79 deg
RMSE (Accelerometer only): 5.05 deg
RMSE (Complementary):      0.53 deg

ジャイロのみの推定（RMSE 9.79°）はバイアスドリフトにより時間とともに真値から大きく乖離し、加速度計のみの推定（RMSE 5.05°）は振動ノイズで大きく変動します。相補フィルタ（RMSE 0.53°）はこれら両方の問題を軽減し、単独センサ利用時と比べてジャイロ比で約18分の1、加速度計比で約9.5分の1にRMSEを低減できることが確認できます。図の4段目を見ると、相補フィルタの出力（緑線）は真値（グレー破線）にほぼ完全に重なっており、\(t=15\) s のステップ変化にも大きなオーバーシュートなく追従しています。

パラメータ \(\alpha\) の影響

\(\alpha\) （時定数 \(\tau\) ）の値が推定精度に与える影響を分析します。

# --- 異なる α（τ）での比較 ---
tau_values = [0.1, 0.5, 1.0, 5.0, 10.0]
colors = [BLUE, ORANGE, GREEN, '#4a3aa7', RED]

fig, axes = plt.subplots(2, 1, figsize=(12, 10))

for tau_val, c in zip(tau_values, colors):
    alpha_val = tau_val / (tau_val + dt)
    theta_est = np.zeros(N)
    theta_est[0] = theta_true[0]
    for i in range(1, N):
        theta_est[i] = alpha_val * (theta_est[i - 1] + gyro_output[i] * dt) \
                      + (1 - alpha_val) * accel_output[i]

    rmse = np.sqrt(np.mean((theta_est - theta_true) ** 2))
    axes[0].plot(t, theta_est, color=c, label=f'τ={tau_val}s (α={alpha_val:.4f}, RMSE={rmse:.2f}°)')

axes[0].plot(t, theta_true, 'k--', linewidth=2, alpha=0.7, label='True')
axes[0].set_ylabel('Angle [deg]')
axes[0].set_xlabel('Time [s]')
axes[0].set_title('Effect of Time Constant τ on Complementary Filter')
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)

# τ vs RMSE のプロット
tau_range = np.logspace(-2, 2, 100)
rmse_list = []
for tau_val in tau_range:
    alpha_val = tau_val / (tau_val + dt)
    theta_est = np.zeros(N)
    theta_est[0] = theta_true[0]
    for i in range(1, N):
        theta_est[i] = alpha_val * (theta_est[i - 1] + gyro_output[i] * dt) \
                      + (1 - alpha_val) * accel_output[i]
    rmse_list.append(np.sqrt(np.mean((theta_est - theta_true) ** 2)))

axes[1].semilogx(tau_range, rmse_list, color=BLUE, linewidth=2)
axes[1].set_xlabel('Time Constant τ [s]')
axes[1].set_ylabel('RMSE [deg]')
axes[1].set_title('RMSE vs Time Constant τ')
axes[1].grid(True, which='both', alpha=0.3)

# 数値探索による最適τ
optimal_idx = np.argmin(rmse_list)
optimal_tau = tau_range[optimal_idx]
optimal_rmse = rmse_list[optimal_idx]
axes[1].axvline(optimal_tau, color=RED, linestyle='--',
                label=f'Empirical optimal τ={optimal_tau:.2f}s (RMSE={optimal_rmse:.2f}°)')

# 式(12)による解析的τ*（ゼロ平均ノイズのみを仮定、バイアスは無視）
Q = (gyro_noise_std * dt) ** 2   # ジャイロ雑音の1ステップ積分分散
R = accel_noise_std ** 2         # 加速度計の角度換算分散
r = Q / R
alpha_star = 1 + r / 2 - np.sqrt(r ** 2 / 4 + r)
tau_star = alpha_star * dt / (1 - alpha_star)
axes[1].axvline(tau_star, color='#4a3aa7', linestyle=':',
                label=f'Analytic τ*={tau_star:.2f}s (noise-only, bias ignored)')
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.savefig('alpha_effect.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"Q={Q:.6f} deg^2, R={R:.3f} deg^2, r=Q/R={r:.8f}")
print(f"Empirical optimal tau = {optimal_tau:.4f} s (RMSE={optimal_rmse:.3f} deg)")
print(f"Analytic alpha* = {alpha_star:.6f}, tau* = {tau_star:.4f} s")

seed=42 の実行結果は次の通りです。

Q=0.000400 deg^2, R=25.000 deg^2, r=Q/R=0.00001600
Empirical optimal tau = 0.7221 s (RMSE=0.493 deg)
Analytic alpha* = 0.996008, tau* = 2.4950 s

数値探索で見つかった最適時定数は \(\tau \approx 0.72\) s（RMSE 0.49°）ですが、式 \((12)\) が予測する解析的最適値は \(\tau^* \approx 2.50\) s であり、両者は一致しません。これは理論の誤りではなく、式 \((8)\) 〜\((12)\) の導出がゼロ平均の過程ノイズのみを仮定しており、シミュレーションに含まれるジャイロの系統バイアス \(0.5\) deg/s を考慮していないためです。実際、gyro_bias = 0 として同じ探索をやり直すと、数値最適 \(\tau\) は \(2.97\) s まで上昇し、解析値 \(2.50\) s に近づくことを確認しました（同条件・seed=42）。バイアスが存在すると \(\tau\) を大きくする（ジャイロを信頼する）ほど積分によりバイアス起因の誤差が蓄積するため、雑音統計だけから決まる理論値よりも小さい \(\tau\) が実用上有利になります。この結果は次のエッジケースの節で扱う「ジャイロバイアスは相補フィルタ／定常カルマンフィルタの前提を破る」という論点と直接つながっています。

\(\alpha\) の選択指針:

最適な \(\tau\) はセンサのノイズ特性に依存します。ジャイロのバイアスが大きい場合は \(\tau\) を小さく、加速度計のノイズが大きい場合は \(\tau\) を大きく設定します（式 \((12)\) はバイアスがない理想条件での上限的な指針として使い、実データでは必ず数値的な検証を併用してください）。

カルマンフィルタとの比較

相補フィルタとカルマンフィルタはどちらもセンサフュージョンに用いられますが、設計思想が異なります。

選択の指針:

• 計算リソースが限られるマイコン（Arduino等）→ 相補フィルタ
• ノイズ特性が既知で最適推定が必要 → カルマンフィルタ
• プロトタイピングや教育目的 → 相補フィルタ（理解しやすい）
• 多センサ統合（GPS + IMU等） → カルマンフィルタ

実際には、1次の相補フィルタは定常カルマンフィルタの特殊ケースと見なすこともできます。カルマンフィルタのゲインが収束した定常状態では、更新式が式 \((5)\) と同形になることを「理論」の節の式 \((8)\) 〜\((12)\) で厳密に導出し、雑音統計量から解析的最適時定数 \(\tau^*\) を計算できることを示しました。ただし数値実験（式 \((12)\) の節）で確認した通り、ジャイロバイアスが存在する現実的な条件では解析値と実測最適値がずれるため、実運用では必ず数値的なチューニングを併用してください。

エッジケースと実装上の注意

相補フィルタは実装が単純である一方、次のような落とし穴や適用限界があります。

ジャイロバイアスは理論の前提を破る

式 \((8)\) 〜\((12)\) の定常カルマン等価性の導出は、過程ノイズ \(Q\) をゼロ平均と仮定しています。しかし実際のMEMSジャイロには温度依存の系統バイアスがあり、これは積分されると時間に比例して増大する誤差（ドリフト）になります。「パラメータ \(\alpha\) の影響」の数値実験で見た通り、バイアス \(0.5\) deg/s が存在すると数値的な最適時定数（\(\tau\approx0.72\) s）は理論値（\(\tau^*\approx2.50\) s、バイアス無視）より大幅に小さくなります。バイアスが未知・時変の場合は、相補フィルタ単体ではこれを打ち消せないため、事前キャリブレーションや、バイアスを追加の状態変数として推定するバイアス補償型（bias-augmented）カルマンフィルタ／相補フィルタの利用を検討する必要があります。

可変サンプリング周期（ジッタ）への対応

式 \((6)\) の \(\alpha = \tau/(\tau+\Delta t)\) は固定 \(\Delta t\) を前提としています。組み込みシステムではタスクスケジューリングの遅延で \(\Delta t\) が実行のたびに変動する（ジッタが乗る）ことが珍しくありません。この場合、\(\alpha\) を起動時に1回だけ計算して固定すると、瞬間的な \(\Delta t\) の増減で相補特性（\(H_{HPF}+H_{LPF}=1\) ）が崩れます。実装上は、ループのたびに実測 \(\Delta t_n\) を計測し、毎回 \(\alpha_n = \tau/(\tau+\Delta t_n)\) を再計算するのが安全です。

角度のラップアラウンド（±180°境界）

ヨー角のように \(\pm 180°\) で不連続に折り返す量に単純な線形補間の式 \((5)\) をそのまま適用すると、真値が \(179°\) から \(-179°\) へ移った瞬間に \(358°\) 分の誤差として誤認識され、フィルタが暴走します。対策としては、角度差分を atan2(sin(Δθ), cos(Δθ)) のように \((-180°, 180°]\) に正規化してから式 \((5)\) を適用するか、後述のようにクォータニオン表現に切り替える方法があります。今回のシミュレーションは \(\pm 30°\) 程度の範囲に収まるようテスト信号を設計しており、この問題は意図的に回避しています。

加速度計は「重力＋並進加速度」を区別できない

加速度計から傾斜角を計算する式は、センサが検出する加速度ベクトルが重力のみであることを仮定しています。しかし車両やドローンが並進加速度を伴う運動（急加速・急旋回）をすると、加速度計出力は重力方向からずれ、傾斜角の推定に系統誤差が生じます。この誤差はローパスフィルタで完全には除去できません（並進加速度の周波数成分がクロスオーバー周波数 \(f_c\) 以下に重なる場合、姿勢角の低周波誤差として通過してしまう）。高ダイナミクスの用途では、GPS速度や対気速度センサから並進加速度を推定し補正する拡張が必要です。

ロール・ピッチ2軸を超える3次元姿勢への拡張とジンバルロック

本記事の式 \((5)\) は1軸（ロールまたはピッチ）のスカラー角を対象としています。3軸姿勢（ロール・ピッチ・ヨー）をオイラー角のまま相補フィルタで扱うと、ピッチが \(\pm 90°\) に近づいた際にロールとヨーが区別できなくなるジンバルロックが発生します。実用上は、姿勢をクォータニオンや回転行列で表現し、ジャイロの角速度をクォータニオンの時間微分方程式で積分し、加速度計（と磁力計）から得られる姿勢との球面線形補間（SLERP）またはMahony/Madgwickフィルタのような非線形フィードバック則で相補的に補正する設計が標準的です（参考文献のMahony et al. 2008が代表例）。

ヨー角は加速度計だけでは観測できない

加速度計は重力方向（鉛直軸）しか検出できないため、鉛直軸まわりの回転であるヨー角には低周波成分を供給できません。ヨー角のドリフト補正には磁力計（コンパス）またはGPS方位、視覚オドメトリなど別の絶対基準源が必要です。「実用例」の表にある9軸IMU構成はこの制約への対処です。

浮動小数点演算による長時間ドリフト

式 \((5)\) の再帰計算を単精度（float32）で長時間（数時間〜数日）実行すると、丸め誤差が蓄積し、理論上のRMSEでは説明できない緩やかなドリフトが生じることがあります。長時間運用する組み込み実装では、可能な限り倍精度（float64）を使うか、定期的に既知の基準姿勢でリセットする設計が推奨されます。

最新研究動向

固定パラメータ \(\alpha\) の限界（特に上記のジャイロバイアス問題や運動状態依存の最適 \(\tau\) ）に対し、近年の研究は主に「ゲインを状況に応じて適応的に変える」方向で発展しています。

• Jiang, P., Liu, C., Cong, H., & Zhang, F. (2024). A fuzzy adaptive complementary filter for attitude estimation based on norm judgment. Measurement and Control. 加速度計・磁力計・ジャイロそれぞれのノルムの規範からの逸脱度をファジィ推論に入力し、運動状態（静止・低ダイナミクス・高ダイナミクス）に応じて相補フィルタのゲインをオンラインで切り替える手法を提案しています。本記事の数値実験で確認した「バイアス・運動状態によって最適 \(\tau\) が理論値からずれる」問題に対する実践的な解決アプローチです。
• Yamagishi, S., & Jing, L. (2026). Quaternion-Averaging-Based Adaptive Complementary Filter for Pedestrian Dead Reckoning With a Foot-Mounted AHRS. arXiv preprint arXiv:2607.05451. 足装着型AHRS（Attitude and Heading Reference System）のための相補フィルタで、線形補間より厳密なクォータニオン平均化による姿勢統合と、歩行フェーズ・磁気外乱の大きさに応じた適応的な重み付けを組み合わせています。カルマンフィルタ系の手法に匹敵する低誤差を、より低い計算コストで達成したと報告しており、本記事で述べたジンバルロック回避（クォータニオン化）と適応ゲインの両方を実運用レベルで統合した最新事例といえます。

いずれも「相補フィルタの計算コストの低さ」という利点を維持しながら、「固定 \(\alpha\) では運動状態やセンサ特性の変化に追従できない」という本記事で数値的に確認した弱点を補う方向性で一致しています。

実用例

関連記事

• ローパスフィルタの設計と比較：移動平均・バターワース・チェビシェフ - 相補フィルタのLPF部分の理論的背景を解説しています。
• ハイパスフィルタの設計とPython実装 - 相補フィルタのHPF部分の理論的背景を解説しています。
• カルマンフィルタの理論とPython実装 - 相補フィルタの上位互換として位置づけられるカルマンフィルタを解説しています。
• 指数移動平均（EMA）フィルタの周波数特性 - 相補フィルタのLPF部分はEMAと等価であり、その周波数特性を解説しています。
• 適応フィルタ（LMS/RLS）の理論とPython実装 - 環境に応じてパラメータを自動調整する適応フィルタを解説しています。
• バンドパスフィルタの設計とPython実装 - 相補フィルタのLPF+HPF分離を1段にした帯域選択フィルタの設計手順を解説しています。
• ベッセルフィルタの理論とPython実装 - 相補フィルタで位相歪みを抑えたい場合に LPF/HPF のプロトタイプとして利用できる、群遅延が最大限平坦な IIR 設計です。
• ウィーナーフィルタの理論とPython実装 - 相補フィルタは SNR の周波数依存性に基づくウィーナーの直感的特殊ケースとみなせます。
• デジタルフィルタ設計指針ハブ：選定の3軸と特性比較 — 本記事の相補フィルタ（センサフュージョン）を他のフィルタ群と比較しながら選定する際の指針として活用できます。

参考文献

• Mahony, R., Hamel, T., & Pflimlin, J.-M. (2008). Nonlinear Complementary Filters on the Special Orthogonal Group. IEEE Transactions on Automatic Control, 53(5), 1203-1218.
• Higgins, W. T. (1975). A Comparison of Complementary and Kalman Filtering. IEEE Transactions on Aerospace and Electronic Systems, AES-11(3), 321-325.
• Jiang, P., Liu, C., Cong, H., & Zhang, F. (2024). A fuzzy adaptive complementary filter for attitude estimation based on norm judgment. Measurement and Control. https://doi.org/10.1177/00202940241227069
• Yamagishi, S., & Jing, L. (2026). Quaternion-Averaging-Based Adaptive Complementary Filter for Pedestrian Dead Reckoning With a Foot-Mounted AHRS. arXiv preprint arXiv:2607.05451. https://arxiv.org/abs/2607.05451
• Arduino and MPU6050 Accelerometer and Gyroscope Tutorial

関連ツール

• DevToolBox - 開発者向け無料ツール集 - JSON整形、正規表現テスターなど85種類以上の開発者向けツール
• CalcBox - 暮らしの計算ツール - 統計計算、周波数変換など61種類以上の計算ツール