はじめに
マイクの信号を ADC で読むと、0 を中心に揺れてほしい波が 2048 付近を中心に揺れています。脈波センサの値は、大きな一定の値の上に小さな拍動が乗り、そのうえ全体がじわじわ動きます。測りたいのは揺れの部分だけで、下駄は要りません。
この記事のゴールは、この下駄を要件から計算して外せるようになることです。使うのはハイパスフィルタで、低い周波数を落として高い周波数を通します。DC オフセット(信号全体に乗った一定のずれ)やドリフト(時間とともに値がじわじわ動くこと)の除去が主な仕事です。便利な一方で、段差のある信号では偽の山と谷を作り、温度のように平均値そのものを測りたい信号では測りたい値まで消してしまいます。どこまで使ってよいかも、式で確かめます。実機は使わず、すべて計算と数値シミュレーションで示します。
- #1 基礎編 — サンプリング、伝達関数、ボード線図、効き目の測り方、C の骨格
- #2 ローパス — 移動平均、1 次 IIR、2 次バターワースを同じカットオフで比べる
- #3 ハイパス(この記事) — DC とドリフトを消す
- #4 バンドパスとノッチ — Q と Audio EQ Cookbook で狙った周波数を扱う
- #5 カルマンフィルタ(予定)
- #6 EKF・UKF(予定)
- #7 パーティクルフィルタ(予定)
- #8 適応フィルタ(予定)
- #9 以降 用途別のフィルタ(予定)— 平滑化(Savitzky–Golay・メディアン・Hampel)、検出(Goertzel・整合フィルタ・ヒルベルト変換)、レート変換(CIC・コム・オールパス)
この記事で分かること: 2 次バターワース HPF の係数がどこから来るのか、1 次の簡易形とどう違うのか、LPF とどんな関係にあるのか、段差でなぜ負に振れるのか、どんな信号に使ってはいけないのかが分かります。
読み終えたときにできること: 残したい信号の周波数から fc を決め、1 次と 2 次のどちらで足りるかを判断して係数を計算できます。計算した値を C のコードに入れ、ボード線図、ステップ応答、白色ノイズで検証できます。
マイコンで ADC の値を読んだことがある方を想定しています。伝達関数やボード線図は次の章で手短に説明するので、#1 基礎編 を読んでいなくても進められます。詳しくは基礎編の z 変換と伝達関数、周波数応答、群遅延、白色ノイズと PSD で導出しています。設計の手順は #2 ローパス と同じなので、先に読んでおくと比べやすくなります。
📘 最初に道具を 4 つだけ
この記事を読むのに必要な道具を、4 つに絞って再掲します。導出は基礎編に任せ、ここでは読み方だけを押さえます。
伝達関数 H(z)。 マイコンのフィルタは、今の出力を過去の入力と出力の重み付き和で作る差分方程式で書けます。1 サンプルの遅れを z^{-1} という掛け算に置き換えると、この式は出力と入力の比 H(z) = Y(z)/X(z) という 1 つの分数になります。分母を 0 にする点が極、分子を 0 にする点が零点です。 z = 1 は DC(0 Hz)に当たるので、そこに零点を置けば DC は完全に消えます。
ボード線図。 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 と呼びます。
群遅延。 位相を周波数で微分して符号を変えた値で、その周波数付近の信号が何秒遅れて出てくるかを表します。周波数ごとに値が違うと、波形の形が変わります。ハイパスが段差の形を変えるのも、この位相の影響を含みます。
PSD。 パワースペクトル密度(PSD)は、周波数 1 Hz あたりのパワーの分布です。白色ノイズはすべての周波数に同じパワーが乗った信号で、フィルタを通すと PSD は周波数ごとに |H|^2 倍になります。出力の分散(ばらつきの 2 乗)を入力の分散で割った値が分散比で、ノイズがどれだけ残るかを表します。RMS(2 乗平均の平方根)は、平均が 0 のノイズでは標準偏差 σ と同じ値です。
時間波形を読む言葉も先に置きます。ステップ応答は、入力を 0 から 1 へ一度に切り替えたときの出力の動きです。アンダーシュートは、動きと反対の向きへ出力が振れることです。過渡は入力を変えた直後の出力が動いている間で、落ち着いたあとが定常状態です。
道具がそろったので、ハイパスを作ります。
🔼 ハイパスフィルタを作る
ハイパスは低い周波数を落とします。使いどころは DC オフセットの除去と、ゆっくりしたドリフトの除去です。ここでは fs = 1 kHz(サンプリング周波数。1 秒間に値を取り込む回数)、fc = 20 Hz で設計します。
2 次バターワース HPF
DC やドリフトをしっかり落としたいときの基本形です。バターワース特性は、通過域(通したい周波数帯)をできるだけ平らに保つ特性です。ローパスと同じく、アナログの 2 次のフィルタ(プロトタイプ)から出発します。
ローパスのプロトタイプは s/\omega_c の分数で書けます。この s/\omega_c を \omega_c/s に置き換えると、低い周波数と高い周波数の役割が入れ替わり、ハイパスになります。分子と分母に (s/\omega_c)^2 を掛けて整理すると、次の形です。
G_{HP}(s) = \frac{s^2}{s^2 + \sqrt{2}\,\omega_c s + \omega_c^2}デジタルへ移すには、#2 ローパス と同じ双一次変換を使います。 s を (2/T)(1 - z^{-1})/(1 + z^{-1}) で置き換える変換で、安定なアナログフィルタを安定なデジタルフィルタへ移します。周波数軸が tan の形で歪むので、 K = \tan(\pi f_c/f_s) と置いて fc を合わせます(プリワープ)。すると s/\omega_c は次の W に置き換わります。
W = \frac{1}{K}\cdot\frac{1 - z^{-1}}{1 + z^{-1}}G_{HP} は W^2/(W^2 + \sqrt{2}W + 1) です。分子と分母に K^2(1 + z^{-1})^2 を掛けると、分子は (1 - z^{-1})^2 だけが残ります。分母は LPF と同じで、 z^{-1} のべきでまとめて定数項の D で割ります。
D = 1 + \sqrt{2}\,K + K^2 a_1 = \frac{2(K^2-1)}{D} a_2 = \frac{1 - \sqrt{2}\,K + K^2}{D}できあがった伝達関数は次の形です。
H(z) = \frac{1}{D}\cdot\frac{1 - 2z^{-1} + z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}分子は (1 - z^{-1})^2 です。z = 1 に零点が 2 個重なるので、DC は完全に 0 になります。係数は b_0 = b_2 = 1/D 、 b_1 = -2/D です。
fc = 20 Hz、fs = 1 kHz で計算します。K = 0.0629147、D = 1.0929330 なので、1/D = 0.9149691 です。
| 係数 | 値 | 係数 | 値 |
|---|---|---|---|
| b0 | 0.9149691 | a1 | −1.8226949 |
| b1 | −1.8299383 | a2 | 0.8371817 |
| b2 | 0.9149691 | 極 | 0.91135 ± 0.08141j |
a1、a2、極は同じ fc の LPF と同じ値です。違うのは分子だけで、分母が決める応答の速さは LPF と変わりません。
LPF とパワーを分け合う
同じ fc の 2 次バターワース LPF と HPF は、パワー(振幅の 2 乗に比例する量)を分け合う関係にあります。全周波数で次の式が成り立ち、数値でも倍精度の誤差の範囲で一致します。
|H_{LP}|^2 + |H_{HP}|^2 = 1式の出どころを確かめます。プリワープした周波数を \Omega = \tan(\pi f/f_s)/K と置くと、z を単位円に置いたとき W は j\Omega になります。分母の絶対値の 2 乗は次のとおりです。
\begin{aligned} & |1 - \Omega^2 + j\sqrt{2}\,\Omega|^2 \\ & = (1 - \Omega^2)^2 + 2\Omega^2 = 1 + \Omega^4 \end{aligned}LPF の分子の絶対値の 2 乗は 1、HPF は \Omega^4 なので、2 つの比を足すと分母と分子が同じになって 1 です。入力のパワーが、fc を境に LPF の出力と HPF の出力へ、取りこぼしなく分かれます。
2 次の計算を減らしたいときは、次の小節の 1 次の形も使えます。
1 次の簡易形(差分+漏れ積分)
もっと簡単な形もあります。入力の差分を、少しずつ漏らしながら積分する形です。
y[n] = a\left(y[n-1] + x[n] - x[n-1]\right)伝達関数は次のとおりです。 a = 1-\alpha と置くと「入力から 1 次 LPF の出力を引いたもの」と完全に同じになります。1 次 LPF は H_{LP}(z) = \alpha/(1 - a z^{-1}) で、#2 で扱った y[n] = y[n-1] + \alpha(x[n] - y[n-1]) の伝達関数です。
H(z) = \frac{a\,(1 - z^{-1})}{1 - a\,z^{-1}} = 1 - H_{LP}(z)右の等号は通分で確かめられます。 1 - H_{LP} の分子は 1 - a z^{-1} - \alpha で、 1 - \alpha = a を入れると a(1 - z^{-1}) です。
実装するなら lpf += alpha * (x - lpf); hpf = x - lpf; の 2 行で済みます。#2 で fc = 20 Hz に合わせた α = 0.11794 をそのまま使うと a = 0.88206 で、−3 dB は 22.9 Hz になります。LPF の fc とは一致しない点に注意してください。
一致しない理由は、1 次の組が分け合うのがパワーではなく振幅だからです。 H_{LP} + H_{HP} = 1 は複素数の和で成り立ちますが、2 乗の和は 1 になりません。fs/2 では z = -1 を入れて、ゲインは次の値です。
|H(-1)| = \frac{2a}{1 + a} = 0.937fs/2 でも 1 に届かず −0.56 dB です。
2 次バターワースと 1 次の簡易形を重ねると、次の図になります。
低域の傾きは 2 次が +40 dB/decade、1 次が +20 dB/decade。ステップ応答の 2 次はいったん −0.21 まで負に振れてから 0 へ戻る
左上の振幅を見ると、1 次は fc の下でゆるやかにしか落ちません。2 Hz の減衰は、2 次の −40.0 dB に対して 1 次は −20.6 dB です。右下のステップ応答では、1 次は 0 へ単調に戻り、2 次だけが負の側へ振れます。
ステップで負に振れる理由
右下の図の違いは、式から説明できます。ステップ入力の z 変換は 1/(1 - z^{-1}) なので、出力は H(z)/(1 - z^{-1}) です。分子の (1 - z^{-1}) が 1 個打ち消されます。
1 次では a/(1 - a z^{-1}) が残り、出力は y[n] = a^{n+1} です。a が 0〜1 の間にあるので、0 へ単調に減ります。fc = 20 Hz の a = 0.88206 なら、最初の値は 0.88 です。
2 次では次の式が残ります。
Y(z) = \frac{b_0\,(1 - z^{-1})}{1 + a_1 z^{-1} + a_2 z^{-2}}分母の極は 0.91135 ± 0.08141j の複素数の組です。実軸から離れた極の組は、減衰しながら振動する成分を作ります。振動なので、0 を横切って負の側へ行きます。半径は 0.915 で、1 サンプルごとに振幅が 0.915 倍になります。 n = 0 の値は b_0 で、最終値は H(1) = 0 です。fc = 20 Hz の計算では、9 ms で 0 を横切り、17 ms で −0.21 の底に達します。
2 次の HPF は、入力が 0 から 1 へ跳ねると出力が 0.91 まで跳ねたあと、−0.21 までアンダーシュートします。ドリフトだけを消すつもりが、段差のある信号では偽の山と谷を作ります。段差そのものが意味を持つ信号には使えません。
この −0.21 は fc によらずほぼ同じ値です。fc を変えても、時間軸が伸び縮みするだけで形はほとんど変わらないからです。次の章では、ハイパスを別の図で見直します。
📊 ナイキスト線図で見る
ボード線図は振幅と位相を別々のグラフに分けて描きます。ナイキスト線図は、同じ情報を複素平面の 1 本の曲線にまとめた図です。 複素平面は、複素数を横軸が実部、縦軸が虚部の点で表す平面です。 周波数を 0 から f_s/2 まで動かしながら、 H(e^{j\omega T}) の点を描いていきます。原点からの距離が振幅、実軸からの角度が位相です。
この図で見るのは、周波数を上げたときに曲線がどこから出発してどこへ向かうかです。灰色の円は半径 1 の単位円で、内側へ入るほど振幅は小さくなります。黒丸に周波数を添え、矢印で周波数が上がる向きを示しました。
#2 のローパス とは逆に、原点から出発して 1 へ向かう。20 Hz の点は虚軸の上にある
ハイパスは 0 Hz で原点から出発し、周波数を上げると 1 へ向かいます。ローパスの曲線を逆向きにたどる形です。20 Hz の点が虚軸の上にあるのは、2 次バターワースの位相が fc でちょうど +90° だからです。 \Omega = 1 を W^2/(W^2 + \sqrt{2}W + 1) に入れると、 W = j なので -1/(\sqrt{2}\,j) = j/\sqrt{2} で、大きさ 0.707、角度 +90° です。
形の読み方が分かったところで、次の章では、低い成分が消える様子を時間波形で確かめます。
🌊 低い成分が消える様子を波形で見る
ボード線図の倍率が実際の波形で何を意味するのかを、時間の軸で確かめます。入力は 3 つの成分を足した信号です。大きさ 1 の DC、振幅 1 で 2 Hz のゆっくりしたうねり、振幅 0.3 で 100 Hz の細かい振動です。センサで言えば、オフセットとゆっくりした揺れの上に、測りたい振動が乗った形です。これを fs = 1 kHz、fc = 20 Hz の 2 次バターワース HPF に通します。
次の図の上の 2 段は、同じ信号の HPF の前と後です。縦軸をそろえたので、高さの差がそのまま成分の差です。オレンジの線は、入力のうち低い成分だけを描いたものです。下の段は出力の一部を時間方向に拡大し、入力に足した 100 Hz と重ねました。
上段の大きなうねりは中段で消え、0 を中心にした 100 Hz の振動だけが残る。下段の拡大では高さが元の 100 Hz とそろい、位相だけが 0.44 ms 進んでいる
入力は −0.29〜2.29 の範囲を上下していますが、出力は過渡のあと ±0.31 に収まります。オレンジの線で描いた低い成分だけが抜け、100 Hz の振動は高さを保ったまま 0 の線のまわりへ移りました。成分ごとの倍率は次の表のとおりです。出力の振幅は、1〜2 秒の区間の波形を成分に分解した値です。
| 成分 | 入力の振幅 | HPF の倍率 | 出力の振幅 |
|---|---|---|---|
| DC | 1 | 0 | 0 |
| 2 Hz | 1 | 0.00997(−40.0 dB) | 0.00997 |
| 100 Hz | 0.3 | 0.9993(−0.006 dB) | 0.2998 |
3 つの倍率は、ここまでの話で説明できます。DC は z = 1 の零点があるので、倍率は厳密に 0 です。2 Hz は fc の 1/10 で、+40 dB/decade の傾きの上にあるので、振幅は約 1/100 になります。100 Hz は fc の 5 倍で、ほとんどそのまま通ります。
中段の左端の揺れは過渡です。フィルタの内部状態を 0 から始めると、入力が 0 から 1 へ跳ねたのと同じことになるので、先ほどのステップ応答と同じ振れが出ます。過渡の分だけを取り出すと、18 ms で −0.18 まで振れ、41 ms で ±0.01 以内に収まります。2 Hz と 100 Hz の立ち上がりも重なるので、ステップ単体の −0.21 とは少し違う値です。起動直後の約 40 ms は、出力をそのまま判定に使わないほうが安全です。
下段で青い線が破線より少し先を走っているのは、HPF の位相の進みです。100 Hz では +15.9° で、周期 10 ms の 15.9/360 にあたる 0.44 ms だけ時刻が前へずれます。複数のチャネルの時刻を比べるときは、全チャネルに同じフィルタを通すと、このずれがそろいます。
ここまでは周波数がはっきりした成分だけの信号でした。次の章では、すべての周波数が均等に混ざった白色ノイズを通します。
🎲 白色ノイズを通す
実際のセンサノイズは、あらゆる周波数が混ざった信号です。その代表である白色ノイズを HPF に通します。条件は、fs = 1 kHz、N = 2^16 点、σ = 1 のガウス白色ノイズです(乱数列を決める初期値のシードは固定)。導出と測り方は基礎編の 白色ノイズと PSD と Welch 法 にまとめました。
時間波形
次の図は、同じ 1 本の白色ノイズを 2 次バターワース HPF(20 Hz)に通し、過渡が収まったあとの 0.5 秒ぶんを切り出したものです(8000 サンプル目から)。縦軸は 2 段共通なので、高さの差がそのまま振幅の差です。
縦軸は 2 段共通。HPF の出力は生の波形とほとんど同じに見える
先ほどの波形とは違い、HPF の出力は生の波形とほとんど同じに見えます。RMS は 0.978 倍、分散比(出力の分散を入力の分散で割った値)は 0.9562 です。白色ノイズのパワーは 0〜500 Hz に均等に広がっているので、20 Hz より下にあるのは 20/500 = 4 % だけです。HPF が捨てるのは、ほぼこのゆっくりした揺れの分だけで、95.6 % は通ります。#2 の LPF の分散比 0.0438 と足すと 1 になるのは、先ほどのパワーの分け合いの結果です。
ここから分かるのは、HPF はノイズを減らすフィルタではないことです。低い側を捨てるだけで、高い周波数のノイズはそのまま通します。
Welch 法の PSD と理論値
どの周波数がどれだけ残ったかを、周波数ごとの数字にします。推定には、短い区間の FFT(信号を周波数成分に分ける計算)を平均する Welch 法を使います。条件はハン窓(区間の端を 0 に絞る重み)、4096 点、50 % オーバーラップ(区間を半分ずつ重ねて取ること)で、理由は基礎編の Welch 法 で説明しています。
理論値は、入力の PSD に |H|^2 を掛けて計算した値です。σ = 1 の生のノイズの水準は 10\log_{10}(2/f_s) = -27.0 dB です。次の図では、シミュレーションの結果を実線、理論値を破線で重ねました。見るのは、低い周波数で実線が破線に沿って下がるかどうかです。
20 Hz より下で +40 dB/decade の傾きに沿って落ち、上では生のノイズに重なる
シミュレーションと理論の差は平均で約 0.7 dB でした。図の傾き +40 dB/decade は、周波数が 10 倍になると振幅が 100 倍になる傾きで、z = 1 の 2 重の零点が作っています。
分散比も理論と比べます。ENBW(等価雑音帯域幅)は、同じ出力パワーになる理想の長方形フィルタの幅で、分散比に f_s/2 を掛けた値です。
| フィルタ | 分散比(理論) | 分散比(シミュレーション) | ENBW |
|---|---|---|---|
| 2 次バターワース HPF 20 Hz | 0.9557 | 0.9562 | 477.8 Hz |
1 次の簡易形の分散比も手で計算できます。インパルス応答は h[0] = a 、 n \ge 1 で h[n] = -\alpha a^n です。2 乗和は等比級数(同じ比で続く数列)の和で、整理すると次の式になります。
\sum_{n=0}^{\infty} h[n]^2 = \frac{2a^2}{1 + a}a = 0.88206 なら 0.827 です。同じ α の 1 次 LPF の分散比は 0.0627 なので、足しても 0.889 で 1 になりません。1 次の組がパワーを分け合わないことが、ここにも表れています。
測り方まで確かめたので、次の章では、設計の手順にまとめます。
🧪 設計手順と実装
ハイパスの設計は、#2 ローパス と同じ 4 段で進めます。
- 要件:fs、残したい信号のいちばん低い周波数、消したい成分の周波数、許せる過渡の長さを書き出す
- 式:1 次なら \alpha = 1 - e^{-2\pi f_c T} 、2 次バターワースなら K と D から係数を計算する
- 係数:C の初期化子や
alphaの定数にする - 検証:ボード線図で残したい周波数の減衰、ステップ応答で過渡、白色ノイズで分散比を確かめる
ハイパスで大事なのは 1 段目です。fc の基準になるのは消したいノイズではなく、残したい信号の下限です。
設計例 1:脈波の DC 除去
要件は、fs = 1 kHz、拍動の 1 Hz を残す、DC とゆっくりしたドリフトを消す、切れ味は要らない、です。切れ味が要らないので 1 次の簡易形を選びます。
音の用途なら fc = 5 Hz が候補です。 \alpha = 1 - e^{-2\pi \cdot 5/1000} = 0.0309 、a = 0.9691 で、−3 dB は 5.16 Hz です。ところが、このままでは 1 Hz が −14.3 dB 落ちます。要件を満たしません。
fc を 0.5 Hz まで下げると、 \alpha = 1 - e^{-2\pi \cdot 0.5/1000} = 0.00314 で、1 Hz の減衰は −0.98 dB に収まります。次の図の左が、この 2 つの検証です。入力は DC 2.0、ゆっくりしたドリフト、振幅 0.1 の 1 Hz の拍動を足したものです。
左:fc = 5 Hz では拍動が 1/5 に痩せ、fc = 0.5 Hz では元の高さに近い。右:段差に近い入力に対し、0.96 まで跳ねて −0.21 まで振れてから収まる
左の図で、fc = 5 Hz の線は破線の 1/5 ほどの高さしかありません。fc = 0.5 Hz の線は破線に近い高さで、位相が少し進んでいます。下駄は消え、ドリフトの分だけ 0 からわずかにずれています。1 次の HPF は、一定の傾きで動くドリフトに対して小さな一定値を残すからです。
設計例 2:振動・ノック検出
要件は、fs = 200 Hz、重力の DC と手の傾き(0.5 Hz 程度)を消す、叩いたときの速い振れ(5 Hz 以上)を残す、です。0.5 Hz と 5 Hz の差は 10 倍しかないので、1 次では分けきれません。2 次バターワースで fc = 2 Hz とします。
K = \tan(\pi \cdot 2/200) = 0.0314263 、D = 1.0454311 なので、b0 = b2 = 0.9565432、b1 = −1.9130865、a1 = −1.9111971、a2 = 0.9149758 です。
検証します。0.5 Hz は −24.1 dB、5 Hz は −0.11 dB で、要件どおりです。ステップ応答は上の図の右で、最初の値は b0 の 0.96、90 ms で 0 を横切り、175 ms で −0.21 の底に達し、410 ms で 1 % 以内に収まります。
C のコードに入れる値
基礎編の C の biquad の biquad_t に、次の係数を入れます。biquad_step() と biquad_reset() はそのまま使えます。
// 2 次バターワース HPF(fc = 20 Hz, fs = 1 kHz)
static biquad_t hpf = {
0.9149691441f, -1.8299382882f, 0.9149691441f,
-1.8226949252f, 0.8371816513f, 0.0f, 0.0f
};
// 2 次バターワース HPF(fc = 2 Hz, fs = 200 Hz)ノック検出の例
static biquad_t knock_hpf = {
0.9565432256f, -1.9130864511f, 0.9565432256f,
-1.9111970674f, 0.9149758348f, 0.0f, 0.0f
};
1 次の簡易形は 2 行です。
lpf += alpha * (x - lpf); // 脈波は alpha = 0.00314f
hpf = x - lpf; // lpf の初期値は最初の x にする
lpf を最初の入力で初期化するのは、起動直後の過渡を消すためです。biquad の biquad_reset() と同じ考え方で、これを忘れると DC の下駄が大きいほど長い時間振れます。
| 条件 | 種類 | 値 |
|---|---|---|
| fs 1 kHz、fc 5 Hz | 1 次の簡易形 | α = 0.0309(−3 dB は 5.16 Hz) |
| fs 1 kHz、fc 0.5 Hz | 1 次の簡易形 | α = 0.00314 |
| fs 200 Hz、fc 2 Hz | 2 次バターワース | 上のコードの係数 |
設計の手順がそろいました。次の章で、よくある場面に当てはめます。
💡 活用事例
ハイパスが効く場面と、効かない場面を並べます。まず表で全体を見て、そのあと 1 つずつ理由を説明します。
| 用途 | フィルタ | 設計値の例 |
|---|---|---|
| 温度センサ | LPF(HPF は不適) | 変化の速さに合わせる |
| マイク・脈波の DC 除去 | HPF(1 次簡易形) | fs 1 kHz、fc 5 Hz、α = 0.0309(脈波は fc 0.5 Hz、α = 0.00314) |
| 振動・ノック検出 | HPF(2 次) | fs 200 Hz、fc 2 Hz |
温度センサのドリフト。 「ドリフトだから HPF」は誤用の代表です。温度そのものがゆっくり動く信号なので、HPF を通すと測りたい値が消えてしまいます。ここで使うのは LPF(ローパスフィルタ)で、ADC のノイズを落とします。電源投入時のドリフトは、フィルタではなく校正や待ち時間の問題です。
マイクや脈波の DC 除去。 こちらは HPF が素直に効く場面です。マイクやピエゾ(圧電素子)の出力は 0 V を中心に上下に揺れるので、回路で VCC/2(電源電圧の半分)付近へ持ち上げて ADC に入れます。脈波センサの出力も、大きな一定の成分の上に小さな拍動が乗った形です。読んだ値から一定のオフセットとゆっくりしたドリフトを HPF で取り除き、0 を中心にそろえてから RMS や FFT にかけます。音や脈波のように、信号そのものが揺れていて平均値に意味が無いものが HPF の対象で、平均値こそが測りたい値である温度センサとはここが違います。切れ味は要らないので、lpf += alpha * (x - lpf); hpf = x - lpf; の 1 次の簡易形で十分です。fs = 1 kHz、fc = 5 Hz なら
\alpha = 1 - e^{-2\pi \cdot 5/1000} = 0.0309
、
a = 1 - \alpha = 0.9691
で、実際の −3 dB は 5.16 Hz です。ただし脈波は心拍の 1 Hz 前後が主役で、fc = 5 Hz のままでは 1 Hz が −14.3 dB 落ちて信号ごと消えます。脈波では fc を 0.5 Hz まで下げて α = 0.00314 とし、1 Hz の減衰を −0.98 dB に抑えます。
振動・ノック検出。 #2 の傾き表示と同じ加速度センサで、残す側を逆にした使い方です。傾き表示は重力のゆっくりした成分を残して振動を捨てましたが、ここでは重力(DC)と手の傾きのゆっくりした変化を HPF で消し、残った速い振れをしきい値や RMS で判定します。fs = 200 Hz、fc = 2 Hz の 2 次バターワース HPF なら、DC は 0 になり、0.5 Hz のゆっくりした傾きは −24.1 dB 落ち、5 Hz の減衰は −0.11 dB で済みます。叩いたときの急な立ち上がりは段差に近い入力です。段差の高さを 1 とすると、設計例のとおり 2 次の出力は 0.96 まで跳ね、90 ms で 0 を横切り、175 ms で −0.21 まで負に振れて、410 ms で 1 % 以内に収まります。波形を残したい信号では欠点だったこの動きが、検出ではむしろ都合のよいパルス状の山になります。跳ねて 0 へ戻るので、しきい値を越えた山を 1 回の叩きとして数えられます。絶対値で判定するときは、負の振れを数えないように、しきい値を最初の山の 2 割より上に置きます。
どの場面でも、まず「平均値に意味があるか」を問うています。次の章で、うまくいかないときの原因をまとめます。
🛠️ トラブルシューティング
実装で起きやすい症状を、原因と対策の組にしてまとめます。
| 症状 | よくある原因 | 対策 |
|---|---|---|
| DC が消えた | HPF や BPF の誤用 | 必要な帯域に DC があるなら LPF を使う |
| 測りたい信号まで小さくなった | fc が信号の下限より高い | 信号の下限の 1/2 以下に fc を下げる |
| 段差のあとに偽の山と谷が出る | 2 次 HPF のアンダーシュート | 1 次にする。段差に意味があるなら HPF を使わない |
| 起動直後に大きく振れる | 状態変数が 0 から始まる | 最初の値で biquad_reset() する |
| 出力が発散する | 極が単位円の外。係数の符号違いや丸め | 係数から極を計算する。a1, a2 の符号を見直す |
✅ まとめ
- 2 次バターワース HPF は、z = 1 に零点を 2 個置いて DC を完全に消す。分母は同じ fc の LPF と同じ
- 2 次の LPF と HPF はパワーを分け合い、1 次の組は振幅を分け合う
- 2 次はステップで負に振れる。極が複素数の組で、振動の成分を持つから。−0.21 という値は fc によらない
- DC と fc の 1/10 のうねりを足した信号を通すと、DC は 0、うねりは約 1/100 になり、高い側の振動だけが 0 のまわりに残る
- HPF はノイズを減らすフィルタではなく、低い側を捨てるフィルタ。白色ノイズの分散比は 0.956
- fc の基準は残したい信号の下限。平均値が測りたい値なら HPF は使わない
ここまでで、DC とドリフトを消す係数を要件から計算し、使ってはいけない場面も説明できるようになりました。次の一手は、手元のセンサの信号で「平均値に意味があるか」を問い、残したい下限の周波数を書き出すことです。#4 バンドパスとノッチの設計 では、低い側や高い側を丸ごと落とすのではなく、狙った周波数の近くだけを扱います。
よくある質問(FAQ)
Q. 1 次と 2 次の HPF はどう選べばいいですか?
A. 切れ味が要らないなら 1 次で十分です。DC を消すだけなら 2 行で書けて、段差でも負に振れません。消したい成分と残したい成分の周波数が 10 倍程度しか離れていないときは 2 次を選びます。fc = 20 Hz の例で 2 Hz の減衰を比べると、1 次の −20.6 dB に対して 2 次は −40.0 dB です。
Q. カットオフ周波数はどう決めればいいですか?
A. 残したい信号のいちばん低い周波数が基準です。脈波の例では 1 Hz を残すために fc を 0.5 Hz に置き、減衰を −0.98 dB に抑えました。消したい成分の周波数は、そのあとで減衰が足りるかを確かめる値です。
Q. 温度センサのドリフトに HPF を使っていいですか?
A. 使えません。温度そのものがゆっくり動く信号なので、測りたい値が消えます。ノイズを落とすなら LPF で、電源投入時のドリフトは校正や待ち時間で対応します。
Q. 入力の平均値を引くのとどう違いますか?
A. 全データの平均を引く方法は、記録したあとの処理なら単純で確実です。ただし未来のデータを使うので、マイコンで逐次処理するときは使えません。1 次の簡易形の lpf は、その平均を過去のデータだけで追いかける値です。
関連記事
この連載:
- #1 基礎編 — 伝達関数、ボード線図、白色ノイズの測り方、C の骨格
- #2 ローパスフィルタの設計 — 同じ手順でローパスを設計する
- #4 バンドパスとノッチの設計 — 狙った周波数を扱う
センサの入口:
- 加速度センサ KXR94-2050 で作るデジタル水平器 — ノック検出の土台になる加速度センサの読み方
- 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、オンライン版)
- Julius O. Smith III「Physical Audio Signal Processing」Bilinear Transformation
- Alan V. Oppenheim, Ronald W. Schafer「Discrete-Time Signal Processing」3rd ed.(Pearson)
係数の式
設計ツール(SciPy v1.18.0 Manual)