デジタルフィルタ入門 #2｜ローパスフィルタの設計：移動平均・1 次 IIR・2 次バターワースを同じカットオフで比べる

はじめに

ADC の値がガタガタするので「とりあえず 10 回平均」を入れる。線はきれいになります。ところが、同じ値をモータの速度制御に戻すと、今度はループが揺れ始める。ローパスフィルタは、いちばんよく使うのに、いちばん代償が見えにくいフィルタです。

この記事のゴールは、ローパスフィルタを「勘の平均」から「要件から計算して決めるもの」に変えることです。ローパスは低い周波数を通し、高い周波数を落とすフィルタです。よく使う 3 種類（移動平均、1 次 IIR、2 次バターワース）を、落ち始める周波数を 20 Hz にそろえて同じ土俵に並べ、ノイズがどれだけ減り、どれだけ遅れ、制御ループの余裕がどれだけ減るかを数字で比べます。実機は使わず、すべて計算と数値シミュレーションで示します。

この記事で分かること： 3 種類のローパスの伝達関数がどこから来るのか、同じカットオフでも減衰、遅れ、ノイズの残り方がなぜ違うのか、ループの中に入れると位相余裕がどれだけ減るのかが分かります。

読み終えたときにできること： fs、落としたい周波数、許せる遅れの 3 つから、移動平均の点数、1 次 IIR の α、2 次バターワースの係数を計算できます。計算した値を C のコードに入れ、ボード線図、ステップ応答、白色ノイズで検証できます。

📘 最初に道具を 4 つだけ

この記事を読むのに必要な道具を、4 つに絞って再掲します。導出は基礎編に任せ、ここでは読み方だけを押さえます。

伝達関数 H(z)。 マイコンのフィルタは、今の出力を過去の入力と出力の重み付き和で作る差分方程式で書けます。1 サンプルの遅れを z^{-1} という掛け算に置き換えると、この式は出力と入力の比 H(z) = Y(z)/X(z) という 1 つの分数になります。分数の分母を 0 にする点が極、分子を 0 にする点が零点で、極が半径 1 の円（単位円）の内側にあればフィルタは安定です。極が 1 つもない（過去の出力を戻さない）形が FIR、ある形が IIR です。

ボード線図。 H(z) の z に e^{j\omega T} （ \omega = 2\pi f 、T はサンプリング周期）を入れると、周波数 f の正弦波が何倍になり（振幅）、どれだけずれるか（位相）が分かります。この 2 つを、対数の周波数軸で 2 段に描いたのがボード線図です。振幅の倍率 |H| は 20\log_{10}|H| の dB（デシベル）で書き、−3 dB は約 0.707 倍、−20 dB は 0.1 倍です。−3 dB の点をカットオフ周波数 fc と呼びます。−20 dB/decade は、周波数が 10 倍になると振幅が 0.1 倍になる傾きです。

群遅延。 位相を周波数で微分して符号を変えた値 \tau_g = -d\angle H/d\omega で、その周波数付近の信号が何秒遅れて出てくるかを表します。どの周波数でも同じ値なら、波形の形は崩れません。これを直線位相と呼びます。

PSD。 パワースペクトル密度（PSD）は、周波数 1 Hz あたりのパワーの分布です。白色ノイズはすべての周波数に同じパワーが乗った信号で、フィルタを通すと PSD は周波数ごとに |H|^2 倍になります。分散（ばらつきの 2 乗）と標準偏差 σ、RMS（2 乗平均の平方根。平均が 0 のノイズでは σ と同じ）が、ノイズの大きさを表す数字です。

道具がそろったので、3 種類のローパスを作ります。

🔽 ローパスフィルタの 3 つの形

ローパスは高い周波数を落とし、ゆっくりした成分を残します。ここでは fs = 1 kHz（サンプリング周波数。1 秒間に値を取り込む回数）、−3 dB の点を約 20 Hz にそろえた 3 種類を比べます。

移動平均（FIR）

直近 N 点の平均です。もっとも身近なデジタルフィルタで、伝達関数は等比級数（同じ比で続く数列）の和なので、Σ を使わない閉じた形になります。

z = e^{j\omega T} を入れて整理します。分子は e^{-jN\omega T/2} をくくり出し、オイラーの式から出る e^{j\theta} - e^{-j\theta} = 2j\sin\theta を使うと、次の形になります。

分母も、N を 1 に置き換えた同じ形です。先頭の指数の項どうしの比は e^{-j(N-1)\omega T/2} で、絶対値が 1 です。絶対値を取ると消えて、振幅は sinc 関数（sin x / x の形の波）の離散版になります。ディリクレ核と呼ばれる形です。

分子の 1 - z^{-N} の零点は、単位円を N 等分した位置に並びます。z = 1 の零点は分母と打ち消し合うので、残りの N − 1 個が振幅の谷、つまり k f_s/N の周波数を作ります。

−3 dB を 20 Hz に近づける長さは N = 22 で、このとき −3 dB は 20.2 Hz、最初の谷は 45.5 Hz です。係数が左右対称なので直線位相で、遅れは周波数によらず次の値です。

この式は、先ほど絶対値で消えた項の位相 -(N-1)\omega T/2 が ω に比例することから出ます。微分すれば、 \omega の係数がそのまま群遅延です。

N を決める目安も計算しておきます。N が大きいとき、 |H| が 1/√2 になるのは f \approx 0.443\, f_s/N のところです。逆に解くと N \approx 0.443\, f_s/f_c で、fs = 1 kHz、fc = 20 Hz なら 22.15、つまり N = 22 です。

図の位相が途中で跳ねているのは、振幅の符号が谷で反転するためで、群遅延が乱れているわけではありません。一方で弱点もあります。阻止域（落としたい周波数帯）の山が最大 −13 dB 程度までしか下がらず、高い周波数を強く落とす用途には向きません。

移動平均は N 点ぶんのメモリを使います。次の 1 次 IIR は、変数 1 個で似た効果を出します。

1 次 IIR（指数移動平均）

マイコンでいちばん使われるのは、おそらく次の 1 行です。変数 1 個で RC ローパスと同じ働きをします。

入力との差の α 倍だけ追いかける、という素直な式です。基礎編の 連続時間から離散時間へ で、抵抗 R とコンデンサ C の RC ローパスの微分方程式から導いた式と同じ形です。過去の入力の重みが 1 サンプルごとに 1-\alpha 倍ずつ指数関数的に小さくなるので、指数移動平均とも呼びます。z 変換すると伝達関数が出て、極は z = 1-\alpha の 1 個です。

これはアナログの RC ローパス G(s) = 1/(RCs + 1) の親戚です。そこで α を RC から決めます。連続時間で極 s を持つ応答は e^{st} の形で変化し、s が負なら減衰します。サンプリング周期 T ごとに見ると、値は毎回 e^{sT} 倍になります。1 次 IIR の極 z = 1-\alpha も「1 サンプルごとに何倍になるか」を表す数だったので、 z = e^{sT} が離散時間の極に当たります。RC の極 s = -1/RC を入れると z = e^{-T/RC} で、これを 1-\alpha と等しいと置くと、次の α が出ます。

右の形は f_c = 1/(2\pi RC) を入れたものです。連続時間の極を z = e^{sT} で離散時間の極へ移すこの対応は、後の「制御ループに入れると何が起きるか」でも、モータの極を移すときに使います。

−3 dB をぴったり fc に合わせたいときは、 |H|^2 = 1/2 を解きます。上の H(z) に z = e^{j\omega T} を入れ、 \omega = 2\pi f_c 、 c = \cos(2\pi f_c T) と置いて、分母の絶対値の 2 乗を計算します。

これを 1/2 と置いて分母を払うと、α の 2 次方程式になります。

α は正の数なので、解の公式の + の側を取ると次の式です。

fc = 20 Hz、fs = 1 kHz では、どの決め方でも近い値になります。基礎編で導いた後退差分（隣のサンプルとの差で微分を近似する方法）の式はやや低めに出ます。

表の厳密解 α = 0.11794 で特性を描くと、次の図になります。

低周波の群遅延は (1-\alpha)/\alpha サンプルで、この例では 7.5 ms です。この式は、分母を低い周波数で近似すると出ます。 \theta = \omega T が小さいとき e^{-j\theta} \approx 1 - j\theta なので、分母 1 - (1-\alpha)e^{-j\theta} は \alpha + j(1-\alpha)\theta になり、位相は -(1-\alpha)\theta/\alpha です。θ で微分すれば、サンプル数の遅れが出ます。計算はかけ算 1 回、メモリは 1 変数で済みます。

時定数で見ると、fc = 20 Hz は \tau = 1/(2\pi f_c) \approx 8\ \mathrm{ms} に当たります。上の図のステップ応答に 63 % の線は引いていませんが、計算すると入力を変えてから 8 ms で 63 %、40 ms で 99 % まで追いつきます。応答は 1 - e^{-t/\tau} なので、10 % から 90 % までの時間は \tau \ln 9 = 17.5 ms で、次の比較の表の 18 ms はこの値です。

1 次 IIR は軽い代わりに、fc の上の傾きが −20 dB/decade にとどまります。もっと急に落としたいときの基本形が、次の 2 次バターワースです。

2 次バターワース（双一次変換）

1 次より急に落とすには、次数を 2 に上げます。2 次の形はいくつもありますが、最初の候補になるのがバターワースです。通過域（通したい周波数帯）をできるだけ平らに保つのがバターワース特性で、最大平坦特性とも呼びます。振幅に波打ちを作らず、0 Hz から単調に落ちていきます。設計の元にするアナログの 2 次のフィルタ（プロトタイプ）は、制御工学の標準 2 次系で減衰比を \zeta = 1/\sqrt{2} にしたものと同じです。

標準 2 次系の分母は s^2 + 2\zeta\omega_c s + \omega_c^2 です。 \zeta （減衰比）は振動の収まりやすさを表す値で、小さいほど振動的になります。

このアナログの式は、そのままではマイコンで計算できません。デジタルへ移すのが双一次変換で、 s を次の式で置き換えます。

s と z の関係が分子・分母とも 1 次の分数になっているので、双一次（bilinear）変換と呼びます。この変換は s 平面の虚軸を単位円にそのまま写すので、安定なアナログフィルタは安定なデジタルフィルタになります。代わりに周波数軸が歪みます。アナログの周波数 \omega_a とデジタルの周波数 \omega_d の関係は次のとおりです。

この式の出どころは簡単です。 z = e^{j\theta} （ \theta = \omega_d T ）を変換の式に入れ、分子と分母から e^{-j\theta/2} をくくり出すと、次のようになります。

これに 2/T を掛けた値が s = j\omega_a なので、上の関係が出ます。 \theta が小さいうちは \tan(\theta/2) \approx \theta/2 で \omega_a \approx \omega_d ですが、fs/2 に近づくほど tan が伸びます。

そこで、あらかじめ fc を \tan で引き伸ばしておきます。これがプリワープです。 K = \tan(\pi f_c / f_s) と置くと、 G(s) の中の s/\omega_c は次の W に置き換わります。

G(s) の分子と分母を \omega_c^2 で割ると 1/(W^2 + \sqrt{2}W + 1) です。分子と分母に K^2(1 + z^{-1})^2 を掛けて分数を払います。

z^{-1} のべきごとに集め、定数項の D で全体を割って分母の先頭を 1 にすると、係数は次のようになります。

fc = 20 Hz、fs = 1 kHz で計算の流れを追います。 K = \tan(0.02\pi) = 0.0629147 なので、 K^2 = 0.0039583 、 \sqrt{2}K = 0.0889748 、D = 1.0929330 です。これを上の式に入れると、次の表の値になります。

Q は山や谷の鋭さを表す値で、#4 で定義します。この式の結果は scipy.signal.butter の出力、および Audio EQ Cookbook（W3C が公開している biquad の係数式集）の LPF に Q = 1/√2 を入れた結果と、倍精度（64 ビットの浮動小数点数）の丸め誤差の範囲で一致します。

この係数で計算した特性が次の図です。左上の振幅が fc の手前まで平らなことと、右上の極が 1 組（2 個）あることを見てください。

位相は fc でちょうど −90° を通り、高い周波数へ向かって −180° に近づきます。極は 0.91135 ± 0.08141j で、実軸から離れた複素数の組です。この組がステップ応答に振動の成分を持ち込むので、4.3 % の行き過ぎが出ます。1 次 IIR の極は実軸の上の 1 個だったので、行き過ぎませんでした。

プリワープを省くと、fc を上げるほど狙いとのずれが広がります。

これで同じカットオフのローパスが 3 種類そろいました。次の小節で並べて比べます。

3 つを同じ土俵で比べる

ここまでの 3 種類を、−3 dB を 20 Hz にそろえたまま 1 枚に重ねます。この図で見比べるのは、左の振幅で高い周波数がどこまで落ちるか、中央の群遅延で周波数ごとにどれだけ遅れるか、右のステップ応答で動き方がどう違うか、の 3 点です。

図から読める差を数字にしたのが次の表です。表の状態変数は、次のサンプルまで覚えておく値です。

低周波の群遅延も式で読めます。アナログのプロトタイプの位相は、低い周波数で -\sqrt{2}\,\omega/\omega_c なので、群遅延は \sqrt{2}/\omega_c \approx 0.225/f_c 秒です。fc = 20 Hz で 11.3 ms で、表の 11.2 ms に近い値です。

バターワースの群遅延は fc の少し下で 13.6 ms まで膨らみます。カットオフ付近の成分だけが余計に遅れるので、これが波形の形の変化として現れます。移動平均は常に 10.5 ms です。

3 種類がそろったところで、次の章では、制御ループに入れたときの影響を別の図で見ます。

📊 ナイキスト線図で見る

フィルタを制御ループに入れたときの安定性は、ボード線図だけでは読み取りにくいものです。ボード線図は振幅と位相を別々のグラフに分けて描きます。ナイキスト線図は、同じ情報を複素平面の 1 本の曲線にまとめた図です。 複素平面は、複素数を横軸が実部、縦軸が虚部の点で表す平面です。 周波数を 0 から f_s/2 まで動かしながら、 H(e^{j\omega T}) の点を描いていきます。原点からの距離が振幅、実軸からの角度が位相です。

ローパス自身のナイキスト線図

まず、ローパス 1 本だけをこの図に描きます。この図で見るのは、周波数を上げたときに曲線がどこから出発してどこへ向かうかです。灰色の円は半径 1 の単位円で、内側へ入るほど振幅は小さくなります。黒丸に周波数を添え、矢印で周波数が上がる向きを示しました。

ローパスは 0 Hz で 1 から出発し、位相を遅らせながら原点へ巻き込みます。20 Hz の点が虚軸の上にあるのは、2 次バターワースの位相が fc でちょうど −90° だからです。

形の読み方が分かったところで、次の小節ではこの図を制御ループに使います。

制御ループに入れると何が起きるか

センサのフィルタは、制御ループの中に入ると安定性を変えます。この小節では、その変化を位相余裕とゲイン余裕の数字で確かめます。制御ループは、出力を測り、目標との差から次の操作を決める、を繰り返す仕組みです。モータの速度制御を例に、速度センサへ LPF を入れた場合を計算します。モータは、ゲイン K、時定数 τ の 1 次遅れ G(s) で表します。

これを 1 ms 周期の PI 制御で動かします。PI 制御は、目標との差に比例した量と、差を積分した量を足して操作量（モータへの指令）を決める方式です。マイコンの操作量は、次の更新まで同じ値に保たれます。この階段状の保持がゼロ次ホールドです。

ゼロ次ホールドを前提に、プラント（制御対象）を離散化（差分方程式の形に変換）します。1 次遅れの極 s = -1/\tau は、1 次 IIR の節と同じ z = e^{sT} の対応で z = p = e^{-T/\tau} に移ります。操作量 u が 1 サンプルの間一定なら、出力はその間に最終値 Ku へ向かって、残りの差の 1-p の割合だけ進みます。残りの差が p 倍になるので、次の差分方程式です。

今のサンプルの u[n] が出力に現れるのは、次のサンプルの y[n+1] です。両辺を z 変換して比を取ると、この 1 サンプルが分子の z^{-1} になります。

z = 1 を入れると分子と分母の 1-p が消えて K になるので、一定の入力に対する倍率は連続時間のモータと同じです。操作量の計算に手間取り、出すのがさらに 1 サンプル遅れる場合は z^{-1} がもう 1 個増えますが、この記事の計算には入れていません。2 次以上のプラントも含めた離散化の一般形は、参考に挙げた制御の教科書に譲ります。

次は PI 制御そのものです。積分は基礎編の 連続時間から離散時間へ の置き換え \int \to T\sum をそのまま使います。積分した値を I[n] 、目標との差を e[n] と書くと、1 サンプル前までの和に今の値を足す形の I[n] = I[n-1] + T\,e[n] です。z 変換すると遅れの項が z^{-1} の掛け算になるので、積分器は T/(1 - z^{-1}) です。これに比例の項を足したものが、PI 制御の伝達関数 C(z) です。

比例のゲインは K_p = 3 、積分のゲインは K_i = 60 です。 K_p/K_i を \tau に合わせて、プラントの極を PI のゼロでほぼ打ち消しています。センサの LPF を F(z) とすると、ループを一周した伝達関数（ループを 1 か所で切って見る形なので、開ループ伝達関数と呼ぶ）は L = C P F です。フィードバックをつないだ全体を閉ループと呼び、その分母は 1 + L なので、 L = -1 になる周波数があると、信号が一周して同じ形で戻り、減衰しなくなります。

フィルタなし、LPF 20 Hz、LPF 10 Hz の 3 通りについて、開ループ L を描いたのが次の図です。左のナイキスト線図は −1 点の近くを拡大したもので、灰色の破線が単位円です。

fc を下げるほど曲線は −1 点へ近づき、右のボード線図では位相が −180° を通る周波数が下がっていきます。LPF なしを含む 5 通りで余裕を計算すると、次の表になります。

このループのゲイン交差は約 9.6 Hz です。fc = 20 Hz の LPF はゲイン交差の 2 倍の周波数にあり、振幅はほとんど減らしません。それでも位相余裕は 40° 減ります。交差周波数 fgc での 2 次バターワースの位相遅れが、次の式で効くからです。

r = 0.47 を入れると 40° で、表の差と合います。ゲイン余裕の 8.6 dB も時間領域で検算できます。ゲインを 8.6 dB 上げたループをシミュレーションすると振幅がほぼ一定の振動になり、0.3 dB 下げると減衰、0.3 dB 上げると発散しました。

滑らかにするほど制御は不安定に寄ります。 これを時間波形で見たのが次の図です。

上段では、LPF 10 Hz の揺れが 600 ms たっても残っています。下段でフィルタなしの線だけがガタつくのは、σ = 0.02 のセンサノイズが比例のゲイン 3 倍で操作量へ出るためです。この差を数字で比べます。次の表のノイズ RMS は、センサに σ = 0.02 のガウスノイズ（平均の近くに多い釣り鐘形の分布に従うノイズ）を足したときの定常状態の値です。

LPF は操作量のガタつきを 1/4 にしますが、代わりに行き過ぎが 0 % から 23.5 % になります。どこでバランスさせるかが設計です。

次の章では、ループから離れて、ノイズそのものがフィルタでどれだけ減るかを計算します。

🎲 白色ノイズを通す

実際のセンサノイズは、あらゆる周波数が混ざった信号です。その代表である白色ノイズを 3 つのローパスに通します。条件は、fs = 1 kHz、N = 2^16 点、σ = 1 のガウス白色ノイズです（乱数列を決める初期値のシードは固定）。理論はこの記事の冒頭のとおりで、PSD は |H|^2 倍になり、分散はインパルス応答の 2 乗和倍になります。導出と測り方は基礎編の 白色ノイズと PSD と Welch 法 にまとめました。

時間波形

次の図は、同じ 1 本の白色ノイズを 3 つのフィルタに通し、過渡が収まったあとの 0.5 秒ぶんを切り出したものです（8000 サンプル目から）。縦軸は 4 段共通なので、高さの差がそのまま振幅の差です。見比べるのは、波の高さと、線の細かさ（どれくらいの速さで動いているか）の 2 つです。

いちばん上の生の白色ノイズは、毎サンプル独立に値が決まります。前の値が次の値の手がかりにならない、つまり相関（2 つの値がどれだけ連動するかを −1〜1 で表す量）がありません。隣のサンプルどうしの相関は 0.003 で、ほぼ 0 です。

下の 3 段では細かい動きが消え、ゆっくりした揺れだけが残ります。これがノイズの中の 20 Hz 以下の成分です。出力は過去の入力の重み付き和なので隣のサンプルと似た値になり、隣どうしの相関は移動平均で 0.954、1 次 IIR で 0.881、2 次バターワースで 0.993 まで上がります。1 次 IIR の線に細かいギザギザが残るのは、20 Hz より上の落ち方がゆるやかだからです。残った揺れは入力より遅れて出てきます。「3 つを同じ土俵で比べる」で見た群遅延です。

図の中の数字は出力の標準偏差で、RMS 比はそれを生のノイズの標準偏差 0.999 で割った値です。分散比（出力の分散を入力の分散で割った値）とあわせてまとめると次のとおりです。

表の値はこの乱数列のシミュレーションの値で、理論値との比較は「分散は何分の 1 になるか」の表で行います。時間で見えた違いは、どの周波数がどれだけ残ったかの違いです。次の小節で、周波数ごとに確かめます。

Welch 法の PSD と理論値

時間波形で見えた違いを、周波数ごとの数字にします。シミュレーションの出力から PSD を推定し、理論の計算と重なるかを確かめます。推定には、短い区間の FFT（信号を周波数成分に分ける計算）を平均する Welch 法を使います。1 回の長い FFT ではノイズのギザギザが減らないからです。条件はハン窓（区間の端を 0 に絞る重み）、4096 点、50 % オーバーラップ（区間を半分ずつ重ねて取ること）で、理由は基礎編の Welch 法 で説明しています。

理論値は、入力の PSD に |H|^2 を掛けて計算した値です。σ = 1 の生のノイズの水準は 10\log_{10}(2/f_s) = -27.0 dB です。次の図では、シミュレーションの結果を実線、理論値を破線で重ねました。見るのは、実線が破線の上に乗るかどうかです。

シミュレーションと理論の差は、1 次 IIR と 2 次バターワースの平均で約 0.7 dB でした（谷の深い移動平均は対数の差が発散するので除いています）。実線が破線に乗るということは、ボード線図の形がそのままノイズのスペクトルの形になるということです。

PSD の形が理論どおりだと分かったので、次の小節では、その面積である分散が入力の何分の 1 になるかを数字にします。

分散は何分の 1 になるか

ノイズがどれだけ減るかを 1 つの数字で言うために、出力の分散を比べます。あわせて、帯域の広さを 1 つの数字で比べるのが等価雑音帯域幅（ENBW）です。同じ出力パワーになる理想の長方形フィルタの幅を表します。白色ノイズでは分散比が \sum h[n]^2 なので、全帯域の幅 f_s/2 にこの比を掛けて、 B_n = (f_s/2)\sum h[n]^2 （ピークのゲインが 1 のとき）で求まります。

理論の列は手で計算できます。移動平均のインパルス応答は 1/N が N 個なので、2 乗和は N \cdot (1/N)^2 = 1/N です。N = 22 で 0.0455 です。

1 次 IIR のインパルス応答は h[n] = \alpha(1-\alpha)^n です。2 乗和は比が (1-\alpha)^2 の等比級数なので、閉じた形になります。

α = 0.11794 を入れると 0.0627 で、表の値です。

ENBW はカットオフとの比で覚えられます。アナログの 1 次は |H|^2 = 1/(1 + (f/f_c)^2) 、2 次バターワースは |H|^2 = 1/(1 + (f/f_c)^4) で、0 から無限大まで積分するとそれぞれ fc の \pi/2 倍と \pi/(2\sqrt{2}) \approx 1.11 倍です。同じ −3 dB でも、1 次 IIR は裾が広いぶん ENBW が 31.3 Hz と大きく、ノイズが残ります。

RMS で見ると、20 Hz の 2 次バターワースでノイズは 0.209 倍です。ENBW は fc に比例し、RMS はその平方根に比例します。一方で遅れは 1/fc に比例します。

ノイズと遅れの関係が数字になりました。次の章では、ここまでの式を設計の手順にまとめます。

🧪 設計手順と実装

ローパスの設計は、4 段の手順で進めます。要件を書き出し、式で候補を計算し、係数にして、検証する。この順番を守ると、「なんとなく 10 回平均」から抜け出せます。

1. 要件：fs、残したい信号の上限、落としたいノイズの周波数と減衰量、許せる遅れを書き出す
2. 式：移動平均なら N \approx 0.443\, f_s/f_c 、1 次 IIR なら \alpha = 1 - e^{-2\pi f_c T} 、2 次バターワースなら K と D から係数を計算する
3. 係数：C の初期化子や alpha の定数にする
4. 検証：ボード線図で減衰、ステップ応答で動き、白色ノイズで分散比を確かめる

設計例 1：傾き表示

加速度センサで傾きを表示する場面で、手順を通します。要件は、fs = 200 Hz、手の振動の 10 Hz を −14 dB（1/5）以下に落とす、ゆっくりした傾きの変化は 1 Hz まで残す、遅れは 100 ms 以内、とします。fc は残す上限より少し上の 2 Hz を試します。

1 次 IIR の α は次のとおりです。

2 次バターワースは K = \tan(\pi \cdot 2/200) = 0.03143 、D = 1.04543 から、b0 = 0.00094469、a1 = −1.9111971、a2 = 0.9149758 です。2 つの候補を要件と突き合わせると次の表になります。遅れの式は、1 次 IIR が (1-\alpha)/\alpha サンプル、2 次バターワースが FAQ で触れる約 0.225/fc 秒です。

10 Hz は 1 次 IIR でぎりぎり要件を満たし、遅れも収まります。2 次バターワースは減衰に余裕がありますが、遅れが要件を超えます。この要件なら 1 次 IIR です。

次の図で検証します。左の振幅が計算どおりかを見て、右で実際の波形を見ます。右の入力は、0.5 秒で傾きが 0 から 1 へ変わる段差に、10 Hz の振動と小さなノイズを足したものです。

右の図で、1 次 IIR の線には振幅が約 1/5 になった 10 Hz の揺れが残っています。揺れをさらに小さくしたくなったら、要件の遅れを見直す番です。これが「どこでバランスさせるか」の中身です。

設計例 2：エンコーダの速度

もう 1 つ、移動平均の例です。エンコーダの速度を 1 ms ごとに計算し、カウントの刻みを 10 点の移動平均でならします。N = 10、fs = 1 kHz を入れると、−3 dB は 44.5 Hz、最初の谷は fs/N = 100 Hz、遅れは (N-1)T/2 = 4.5 ms、白色ノイズの分散比は 1/N = 0.1 です。分解能の計算は、活用事例の章で続けます。

C のコードに入れる値

基礎編の C の biquad の biquad_t に、次の係数を入れます。biquad_step() と biquad_reset() はそのまま使えます。

// 2 次バターワース LPF（fc = 20 Hz, fs = 1 kHz）
static biquad_t lpf = {
    0.0036216815f, 0.0072433630f, 0.0036216815f,
    -1.8226949252f, 0.8371816513f, 0.0f, 0.0f
};

1 次 IIR は 1 行で書けます。α の値は表のとおりです。

y += alpha * (x - y);   // alpha は次の表から

FPU のないマイコンでは、基礎編の 固定小数点の 1 次 IIR を使います。shift を 1 つ決めると \alpha = 2^{-\mathrm{shift}} なので、α は 2 のべき乗にしか選べません。狙いの α にいちばん近いものを選び、実際の −3 dB を計算し直します。fs = 1 kHz で 20 Hz を狙うと α = 0.11794 なので、近いのは shift = 3（α = 1/8）で、−3 dB は 21.3 Hz です。shift = 6（α = 1/64）なら −3 dB は約 2.5 Hz です。

設計の手順がそろいました。次の章で、よくある場面に当てはめます。

💡 活用事例

ローパスを、電子工作でよくある場面に当てはめます。まず表で全体を見て、そのあと 1 つずつ理由を説明します。

加速度センサの傾き表示。 KXR94-2050 のデジタル水平器 のように、重力の向きから角度を出す用途です。ほしいのは静かな重力成分で、手の振動は不要です。fs = 200 Hz、fc = 2 Hz の 1 次 IIR なら α = 0.061、遅れは約 80 ms で、目で見る表示には十分です。ジャイロと組み合わせる話は #5 のカルマンフィルタにつながります。

PID の D 項。 PID 制御は、比例（P）、積分（I）、微分（D）の 3 項で操作量を決める方式で、D 項はその微分の項です。微分は高い周波数を増幅します。1 サンプルの差分をとるとノイズの標準偏差は √2 倍になり、さらに 1/T が掛かります。fs = 1 kHz なら 1,414 倍です。そのため D 項には LPF が付き物です。ただし LPF の遅れは位相余裕を削るので、先ほどの表の考え方で fc を決めます。

エンコーダの速度。 1000 パルス/回転を 4 逓倍（A 相と B 相の両方のエッジを数えて分解能を 4 倍にすること）し、1 ms ごとに差分をとると、1 カウントが 15 rpm に相当します。移動平均を 10 点とると分解能は 1.5 rpm です。カウント c の差分の移動平均は (c[n] - c[n-10]) / (10 * T) と同じなので、リングバッファ（古い値から順に上書きする配列）1 本で書けます。遅れは 4.5 ms です。

温度センサのドリフト。 「ドリフトだから HPF（ハイパスフィルタ）」は誤用の代表です。温度そのものがゆっくり動く信号なので、HPF を通すと測りたい値が消えてしまいます（#3 で扱います）。ここで使うのは LPF で、ADC のノイズを落とします。電源投入時のドリフトは、フィルタではなく校正や待ち時間の問題です。

どの場面でも、何を残して何を捨てるかからフィルタを選んでいます。次の章で、この考え方を手順にまとめます。

🧭 選び方

フィルタの選び方を、判断の順に並べます。最初に問うのは、消したいものが特定の周波数かどうかです。

図の分かれ道のほかに、遅れやメモリの条件で決まることもあります。条件ごとの向き不向きは次のとおりです。

firwin は SciPy の FIR 設計関数で、基礎編の 直線位相が要るなら FIR で扱っています。ループの中に入れるフィルタは、ボード線図だけでなく位相余裕を計算して決めます。逆に、表示やログのようにループの外で使うフィルタは、遅れの心配が少ないので強めにかけられます。

フィルタを選んで動かしたあと、思いどおりにならないときは、次の章の表から原因を探します。

🛠️ トラブルシューティング

実装で起きやすい症状を、原因と対策の組にしてまとめます。

float の丸めだけでも、fc が低いと効いてきます。fs = 1 kHz の 2 次バターワース LPF で、係数を float に丸めたときの DC ゲイン（0 Hz でのゲイン）を計算すると次のようになります。

0.1 Hz では DC ゲインが 5 % 以上ずれます。 1 + a_1 + a_2 が 0 に近い小さな数になり、float の分解能では表せなくなるためです。fc/fs が極端に小さいときは、double を使うか、先に移動平均で間引いて fs を下げてから設計します。

✅ まとめ

• 同じ fc でも、移動平均は谷と直線位相、1 次 IIR は軽さ、2 次バターワースは切れ味が持ち味
• 2 次バターワースの係数は、双一次変換とプリワープから K と D の式で求まる
• ナイキスト線図で見ると、ループに入れた LPF は開ループを −1 点へ近づける。fc = 20 Hz で位相余裕は 88.4° から 48.0° へ下がった
• 白色ノイズの PSD は理論の |H|² をなぞり、分散はインパルス応答の 2 乗和で決まる。RMS を半分にする代償は遅れ 4 倍
• 設計は、要件、式、係数、検証の順に進める

ここまでで、要件からローパスの係数を計算し、代償を数字で説明できるようになりました。次の一手は、手元のセンサで要件を書き出し、この記事の表と同じ 4 行を自分の値で埋めてみることです。#3 ハイパスフィルタの設計 では逆に、低い周波数を落とすフィルタを同じ手順で設計します。

よくある質問（FAQ）

Q. 移動平均と 1 次 IIR、どちらを使えばいいですか？

A. 迷ったら 1 次 IIR が手軽です。変数 1 個で済み、同じ −3 dB なら遅れも小さくなります。特定の周波数に谷を置きたいときや、遅れを周波数によらず一定にしたいときは移動平均を選びます。ノイズを減らす力は移動平均のほうが上で、同じ −3 dB で分散比は 0.0455 対 0.0627 です。

Q. カットオフ周波数はどう決めればいいですか？

A. 残したい信号の一番速い成分より高く、消したいノイズより低い位置が候補です。そのうえで、遅れが許せるかを群遅延で見ます。2 次バターワースの低周波の遅れは、約 0.225/fc 秒です（fc = 20 Hz で 11.2 ms）。制御ループの中で使うなら、位相余裕の計算まで行います。

Q. 位相余裕を減らさずに LPF を入れられますか？

A. 減らさないことはできませんが、減りを小さくできます。fc をゲイン交差の 5〜10 倍に置くと、この記事の例では位相余裕の減少が約 16°〜8° で済みます。fc を下げるほど急に減り、ゲイン交差の 2 倍では 40° 減ります。

Q. ノイズの RMS を半分にするにはどうすればいいですか？

A. 白色ノイズが相手なら fc を 1/4 にします。RMS は ENBW の平方根に比例し、ENBW は fc に比例するからです。代わりに遅れは 4 倍になります。遅れが許せないときは、ADC のオーバーサンプリング（必要より高い fs で読むこと）など、フィルタの外の対策も考えます。

関連記事

この連載：

• #1 基礎編 — 伝達関数、ボード線図、白色ノイズの測り方、C の骨格
• #3 ハイパスフィルタの設計 — DC とドリフトを消す
• #4 バンドパスとノッチの設計 — 狙った周波数を扱う

フィルタの計算：

• RC フィルタ カットオフ周波数 計算ツール — ADC の前に置くアンチエイリアス RC の定数を求める

センサの入口：

• 加速度センサ KXR94-2050 で作るデジタル水平器 — LPF を入れたくなる代表例。ADC 値から角度を計算する
• Arduino の AD 変換とシリアル通信 — ADC の値を読む基本

参考

取得日はいずれも 2026-11-05 です。図と数値はすべて筆者が Python（NumPy 2.4.6、SciPy 1.17.1、Matplotlib 3.11.2）で計算したものです。

理論

• Julius O. Smith III「Introduction to Digital Filters with Audio Applications」（CCRMA、オンライン版）
• 同「Group Delay」
• Julius O. Smith III「Physical Audio Signal Processing」Bilinear Transformation
• Alan V. Oppenheim, Ronald W. Schafer「Discrete-Time Signal Processing」3rd ed.（Pearson）
• Karl J. Åström, Björn Wittenmark「Computer-Controlled Systems: Theory and Design」3rd ed.（Prentice Hall）

係数の式

• W3C「Audio EQ Cookbook」（Working Group Note、2021 年 6 月 8 日）

設計ツール（SciPy v1.18.0 Manual）

• scipy.signal.butter
• scipy.signal.bilinear
• scipy.signal.freqz
• scipy.signal.group_delay
• scipy.signal.welch