はじめに

ADC の値に 50 Hz のハム(商用電源から配線に乗るノイズ)が乗っている。ローパスで落とそうとすると、信号の帯域まで削れて遅れも増えます。消したいのは 50 Hz だけで、それ以外には触りたくない。逆に、電力の 50 Hz だけを取り出したいこともあります。

この記事のゴールは、狙った周波数だけを扱うフィルタを、Q という 1 つの値から設計できるようになることです。バンドパス(BPF)は狙った帯域だけを通すフィルタで、ノッチは狙った周波数だけを消すフィルタです。どちらも山や谷の鋭さを Q で指定します。Q を上げれば鋭くなりますが、落ち着くまでの時間が長くなり、周波数のずれに弱くなります。その代償を式で追います。実機は使わず、すべて計算と数値シミュレーションで示します。

🗓️ 連載「デジタルフィルタ入門」
  • #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・コム・オールパス)

この記事で分かること: Audio EQ Cookbook の係数の式がどこから来るのか、Q と帯域幅、群遅延、収まる時間がどんな式で結ばれているのか、バンドパスとノッチがなぜナイキスト線図で同じ円を描くのかが分かります。

読み終えたときにできること: 中心周波数と許せる過渡の長さから Q を決め、係数を計算して C のコードに入れられます。周波数のずれに対してノッチがどこまで効くかを計算し、ボード線図、ステップ応答、白色ノイズで検証できます。

💡 対象読者と前提

マイコンで ADC の値を読んだことがある方を想定しています。伝達関数やボード線図は次の章で手短に説明するので、#1 基礎編 を読んでいなくても進められます。詳しくは基礎編の z 変換と伝達関数、周波数応答、群遅延、極と零点、白色ノイズと PSD で導出しています。


📘 最初に道具を 4 つだけ

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

伝達関数 H(z)。 マイコンのフィルタは、今の出力を過去の入力と出力の重み付き和で作る差分方程式で書けます。1 サンプルの遅れを z^{-1} という掛け算に置き換えると、この式は出力と入力の比 H(z) = Y(z)/X(z) という 1 つの分数になります。分母を 0 にする点が極、分子を 0 にする点が零点です。周波数を 0 から fs/2 へ上げることは、複素平面(複素数を横軸が実部、縦軸が虚部の点で表す平面)の半径 1 の円(単位円)の上を歩くことに対応し、極に近づくと山が、零点に近づくと谷ができます。分子も分母も 2 次の形を biquad と呼びます。

ボード線図。 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 倍です。

群遅延。 位相を周波数で微分して符号を変えた値 \tau_g = -d\angle H/d\omega で、その周波数付近の信号の包絡線(細かい波の山をなぞった外形の線)が何秒遅れて出てくるかを表します。

PSD。 パワースペクトル密度(PSD)は、周波数 1 Hz あたりのパワーの分布です。白色ノイズはすべての周波数に同じパワーが乗った信号で、フィルタを通すと PSD は周波数ごとに |H|^2 倍になります。出力の分散(ばらつきの 2 乗)を入力の分散で割った値が分散比で、白色ノイズではインパルス応答 h[n] の 2 乗和に等しくなります。

時間波形を読む言葉も先に置きます。ステップ応答は入力を 0 から 1 へ一度に切り替えたときの出力の動き、インパルス応答は 1 サンプルだけ 1 の入力を入れたときの出力の動きです。リンギングは、出力が振動しながら収まっていく様子です。過渡は入力を変えた直後の出力が動いている間で、落ち着いたあとが定常状態です。

道具がそろったので、Q と係数の式に進みます。


🎯 Q と Cookbook の式

バンドパスとノッチは、どちらも山や谷の鋭さを Q という 1 つの値で指定します。この章では Q を定義し、係数の式を導きます。

💡 ワード解説:Q(品質係数)

Q(クオリティファクタ)は中心周波数を帯域幅(振幅が −3 dB になる 2 点の間の幅)で割った値で、大きいほど鋭くなります。2 次系の減衰比 \zeta (振動の収まりやすさを表す値)とは Q = 1/(2\zeta) の関係にあります。

係数は、W3C の Audio EQ Cookbook(W3C が公開している biquad の係数式集)の式で求めます。f0 と Q を入れるだけで決まり、バンドパスとノッチで分母を共通にできるからです。 \omega_0 = 2\pi f_0/f_s として、まず帯域の広さを決める次の値を求めます。Q が大きいほど小さくなる量です。

\alpha_Q = \frac{\sin\omega_0}{2Q}

分母 A(z) = A_0 + A_1 z^{-1} + A_2 z^{-2} の係数は、バンドパスとノッチで共通です。

\begin{aligned} & A_0 = 1+\alpha_Q,\ A_2 = 1-\alpha_Q \\ & A_1 = -2\cos\omega_0 \end{aligned}

C のコードへ渡すときは、分子も分母も A_0 で割って先頭を 1 にそろえます。なお \alpha_Q は、#2 の 1 次 IIR の α とは別物です。

式の出どころ:アナログのプロトタイプから

Cookbook の式は、アナログの 2 次のフィルタ(プロトタイプ)を双一次変換したものです。中心の角周波数を 1 に正規化すると、プロトタイプは次の 2 つです。

G_{BP}(s) = \frac{s/Q}{s^2 + s/Q + 1} G_{N}(s) = \frac{s^2 + 1}{s^2 + s/Q + 1}

双一次変換は、s を z の式に置き換えてデジタルへ移す変換です。#2 ローパス と同じように、 K = \tan(\omega_0/2) と置いて中心周波数を合わせると、s は次の式に置き換わります。

s \leftarrow \frac{1}{K}\cdot\frac{1 - z^{-1}}{1 + z^{-1}}

分母 s^2 + s/Q + 1 に入れて、 K^2(1 + z^{-1})^2/(1 + K^2) を掛けます。

\begin{aligned} A(z) = {} & \frac{(1 - z^{-1})^2 + (K/Q)(1 - z^{-2})}{1 + K^2} \\ & + \frac{K^2 (1 + z^{-1})^2}{1 + K^2} \end{aligned}

z^{-1} のべきでまとめるときに、tan の半角の公式を使います。 \sin\omega_0 = 2K/(1 + K^2) 、 \cos\omega_0 = (1 - K^2)/(1 + K^2) です。定数項は 1 + K/(Q(1 + K^2)) = 1 + \alpha_Q になり、 z^{-1} の項は -2\cos\omega_0 、 z^{-2} の項は 1 - \alpha_Q です。Cookbook の分母が出ました。

分子も同じ操作です。BPF の s/Q は \alpha_Q(1 - z^{-2}) に、ノッチの s^2 + 1 は 1 - 2\cos\omega_0\, z^{-1} + z^{-2} になります。

H を 1 行に整理する

式が 2 つあると見通しが悪いので、周波数応答を 1 つの実数の変数で書き直します。 z = e^{j\theta} ( \theta = \omega T )を入れ、分子と分母に z を掛けて z + z^{-1} = 2\cos\theta 、 z - z^{-1} = 2j\sin\theta を使います。分子は 2j\alpha_Q\sin\theta 、分母は 2(\cos\theta - \cos\omega_0) + 2j\alpha_Q\sin\theta です。分子で割ると、次の形になります。

H_{BP}(e^{j\theta}) = \frac{1}{1 - jx} x(\theta) = \frac{\cos\theta - \cos\omega_0}{\alpha_Q \sin\theta}

x は実数で、 \theta = \omega_0 で 0 です。この 1 行から、振幅と位相がすぐ読めます。

|H_{BP}|^2 = \frac{1}{1 + x^2},\quad \angle H_{BP} = \arctan x

中心では x = 0 でゲインは 1、位相は 0 です。−3 dB の点は |H|^2 = 1/2 なので x = \pm 1 です。この記事の式の多くは、この 1 行から出ます。

📌 この 1 行から出る 3 つの式
  • 帯域幅:−3 dB の 2 点の間隔を角度で \Delta\theta とすると \tan(\Delta\theta/2) = \alpha_Q

  • 中心の群遅延: \tau_g = 1/\alpha_Q サンプル

  • 白色ノイズの分散比: \alpha_Q/(1 + \alpha_Q)

帯域幅の式は次のように出します。 x = \pm 1 の 2 点を t = \tan(\theta/2) で書くと、 \cos\theta = (1 - t^2)/(1 + t^2) 、 \sin\theta = 2t/(1 + t^2) なので、 c = \cos\omega_0 と置いて次の 2 次方程式になります。

(1 + c)\,t^2 \pm 2\alpha_Q t - (1 - c) = 0

下側の点を t_1 、上側の点を t_2 とします。解と係数の関係(2 次方程式の解の和と積を係数で表す公式)から、 t_2 - t_1 = 2\alpha_Q/(1 + c) 、 t_1 t_2 = (1 - c)/(1 + c) です。tan の差の公式に入れると、1 + c が消えます。

\tan\frac{\theta_2 - \theta_1}{2} = \frac{t_2 - t_1}{1 + t_1 t_2} = \alpha_Q

群遅延は位相の微分です。 \theta = \omega_0 で x は 0 なので d(\arctan x)/d\theta = dx/d\theta で、x の分子を微分すると -\sin\omega_0 、分母はそのまま \alpha_Q\sin\omega_0 です。符号を変えて 1/\alpha_Q サンプルになります。分散比は、ノイズの章で導きます。

式がそろったので、バンドパスから順に形を見ます。


📶 バンドパス

バンドパスは、狙った周波数の近くだけを取り出します。この章では、山の形が Q でどう変わるかを見ます。Cookbook の「constant 0 dB peak gain」の形は、ピークのゲインが Q によらず 1 です。

H_{BP}(z) = \frac{\alpha_Q\,(1 - z^{-2})}{A(z)}

零点は z = 1 と z = −1 にあり、DC と fs/2 を消します。極は角度 \omega_0 の付近にあり、そこに山を作ります。f0 = 100 Hz、Q = 5 のとき、−3 dB の点は 91.0 Hz と 109.7 Hz で、帯域幅は 18.7 Hz です。

f0 = 100 Hz、fs = 1 kHz で確かめます。 \omega_0 = 0.62832 、 \sin\omega_0 = 0.58779 なので、Q = 5 なら \alpha_Q = 0.05878 です。先ほどの式から \Delta\theta = 2\arctan(0.05878) = 0.11743 rad で、周波数に直すと 0.11743 \times 1000/(2\pi) = 18.7 Hz です。中心の群遅延は 1/\alpha_Q = 17.0 サンプル、つまり 17.0 ms です。

Q = 2、5、15 のバンドパスのボード線図、極零点図、ステップ応答

Q を上げると山が鋭くなり、極が単位円へ近づく。中心では位相が 0 を通るので、f0 ちょうどの成分は遅れずに出てくる

Q を 2、5、15 と上げると、−3 dB の帯域幅は 46.4 Hz、18.7 Hz、6.2 Hz と狭くなります(左上)。どの Q でも頂点は 0 dB のままで、裾だけが引き締まります。位相(左下)は f0 より下で進み、f0 で 0 を横切り、f0 より上で遅れます。+90° から −90° への変化は Q が高いほど急です。位相の傾きは群遅延なので、位相が 0 になるのは定常状態の正弦波の話で、包絡線は遅れます。f0 の群遅延は 6.8 ms、17.0 ms、51.0 ms です。

右上の極零点図では、極の半径が 0.862、0.943、0.981 と単位円へ寄っていきます。極が単位円に近いほど応答が長く尾を引く様子は、右下のステップ応答に出ています。入力は 1 のままなのに出力が 0 へ戻るのは、z = 1 の零点が DC を消すためです。最初の振れは Q を上げるほど 0.34、0.17、0.06 と小さくなり、振動が最大値の 1 % 以下に収まるまでの時間は 33 ms、78 ms、232 ms と長くなります。インパルス応答でも同じで、31 ms、80 ms、235 ms です。狭い帯域を選ぶほど、長い時間の信号が要ります。

極の半径も式で書けます。2 つの極は共役(虚部の符号だけが違う組)なので、積は半径の 2 乗で、分母の定数項に等しくなります。

r = \sqrt{\frac{1 - \alpha_Q}{1 + \alpha_Q}}

インパルス応答の振幅は 1 サンプルごとに r 倍になるので、1 % になるサンプル数 n は r^n = 0.01 から n = \ln 0.01/\ln r です。 \alpha_Q が小さければ \ln r \approx -\alpha_Q なので、次の目安になります。

n \approx \frac{4.6}{\alpha_Q}

Q = 2、5、15 では、 r^n = 0.01 の式で 31、78、235 サンプルです。上のインパルス応答の 31 ms、80 ms、235 ms とほぼ合います。Q が決めているのは \alpha_Q で、帯域幅も遅れも収まる時間も、全部この 1 つの数から出てきます。

次のノッチは、この山を裏返した谷を作るフィルタです。


🕳️ ノッチ

ノッチは、バンドパスとは逆に 1 つの周波数だけを消します。ハムのように周波数が分かっている妨害に使います。分母は先ほどの A(z) のままで、分子だけが変わります。

H_{N}(z) = \frac{1 - 2\cos\omega_0\, z^{-1} + z^{-2}}{A(z)}

分子の零点は単位円の上、角度 \pm\omega_0 にあります。したがって f0 ちょうどの正弦波は、計算上は完全に 0 になります。そのすぐ内側に極を置くことで、少し離れた周波数では極と零点が打ち消し合い、ゲインが 1 に戻ります。

Q = 2、10、30 のノッチのボード線図、単位円の近くに極と零点が並ぶ図、ステップ応答

Q を上げると谷は細くなり、拡大図の極が零点へ寄っていく。極が単位円に近づくほど、ステップ応答のリンギングは長く続く

Q は「谷の幅」と「落ち着くまでの時間」のトレードオフ(一方を良くすると他方が悪くなる関係)です。f0 = 50 Hz、fs = 1 kHz で計算すると次のようになります。

Q −3 dB の幅 50.2 Hz の減衰 1 % へ収まる時間
2 24.5 Hz −35.8 dB 58 ms
10 4.9 Hz −21.8 dB 296 ms
30 1.6 Hz −12.5 dB 886 ms

3 列目は、消したい成分が狙いから 0.2 Hz ずれたときの効き方です。Q を上げるほど周辺の信号を傷つけませんが、周波数がずれたときの消え方は急に悪くなります。

表の値も、先ほどの式で読めます。f0 = 50 Hz では \sin\omega_0 = 0.30902 なので、Q = 10 の \alpha_Q は 0.01545 です。帯域幅は 2\arctan(0.01545) \times 1000/(2\pi) = 4.9 Hz、1 % へ収まる時間の目安は 4.6/0.01545 = 298 サンプルで、表の 296 ms に近い値です。

💡 Cookbook と scipy.signal.iirnotch の Q は同じではない

SciPy の iirnotch は Q を「f0 ÷ デジタルの −3 dB 幅」と定義しています。Cookbook の Q はアナログのプロトタイプの値です。低い周波数ではほぼ同じで、f0 = 50 Hz、Q = 10 の係数の差は最大 0.0005 です。

高い周波数では別物になります。f0 = 400 Hz、Q = 5 では、SciPy の −3 dB 幅が 80 Hz なのに対し、Cookbook の式では 18.7 Hz です。f0 が fs/4 を超えるなら、どちらの定義の Q なのかを確かめてください。

差の出どころは \tan(\Delta\theta/2) = \alpha_Q = \sin\omega_0/(2Q) です。 \sin\omega_0 は 100 Hz と 400 Hz で同じ値なので、Cookbook の幅は f0 に比例せず、どちらも 18.7 Hz になります。

次の章では、ここまでのフィルタを別の図で見直します。


📊 ナイキスト線図で見る

ボード線図は振幅と位相を別々のグラフに分けて描きます。ナイキスト線図は、同じ情報を複素平面の 1 本の曲線にまとめた図です。 周波数を 0 から f_s/2 まで動かしながら、 H(e^{j\omega T}) の点を描いていきます。原点からの距離が振幅、実軸からの角度が位相です。

この図で見るのは、周波数を上げたときに曲線がどこを通るかです。灰色の円は半径 1 の単位円で、内側へ入るほど振幅は小さくなります。黒丸に周波数を添え、矢印で周波数が上がる向きを示しました。

バンドパスとノッチの周波数応答を複素平面に描いた 2 枚の図。どちらも中心 0.5、半径 0.5 の円を描く

バンドパスとノッチは、中心 0.5・半径 0.5 のきれいな円を描く

バンドパスとノッチが同じ円になるのは偶然ではありません。分子を足すと分母と等しくなるので、同じ f0 と Q なら全周波数で次の式が成り立ちます。

H_{BP}(z) + H_{N}(z) = 1

ノッチの出力は「入力からバンドパスの出力を引いたもの」と同じです。

円になる理由も、先ほどの 1 行で分かります。

\frac{1}{1 - jx} - \frac{1}{2} = \frac{1}{2}\cdot\frac{1 + jx}{1 - jx}

右辺の (1 + jx)/(1 - jx) は、分子と分母が共役なので絶対値が常に 1 です。つまり H_{BP} は、x がどんな値でも 1/2 から距離 1/2 の位置にあります。これが中心 0.5、半径 0.5 の円です。

バンドパスは 0 Hz で原点を出発し、中心の周波数で 1 に達して、fs/2 で原点へ戻ります。ノッチは 1 - H_{BP} なので、同じ円の上を、1 から出発して f0 で原点を通り、1 へ戻ります。原点を通る、つまり f0 で完全に 0 になることが、図からも読めます。

形の読み方が分かったところで、次の章では、ノイズに対する効き目を計算します。


🎲 白色ノイズを通す

実際のセンサノイズは、あらゆる周波数が混ざった信号です。その代表である白色ノイズを BPF とノッチに通します。条件は、fs = 1 kHz、N = 2^16 点、σ = 1 のガウス白色ノイズです(乱数列を決める初期値のシードは固定)。導出と測り方は基礎編の 白色ノイズと PSD と Welch 法 にまとめました。

時間波形

次の図は、同じ 1 本の白色ノイズを BPF(100 Hz、Q = 5)に通し、過渡が収まったあとの 0.5 秒ぶんを切り出したものです(8000 サンプル目から)。縦軸は 2 段共通なので、高さの差がそのまま振幅の差です。

生の白色ノイズと BPF を通した時間波形の 2 段の図

縦軸は 2 段共通。BPF の出力には周期 10 ms ほどの波が現れ、その振幅がゆっくりふらつく

周期がおよそ 10 ms の波が現れます。入力は乱数なのに正弦波のように見えるのは、フィルタが 91.0〜109.7 Hz の成分だけを残すためです。ただし本物の正弦波とは違い、振幅がゆっくりふらつきます。帯域幅が中心周波数に比べて狭い信号(狭帯域信号)では、近い周波数の成分どうしが強め合ったり打ち消し合ったりするからです。図の中でも、山の高さが 0.5 を超える箇所と 0.1 を下回る箇所があります。RMS(2 乗平均の平方根)は 0.237 倍、分散比(出力の分散を入力の分散で割った値)は 0.0563 です。

Welch 法の PSD と理論値

どの周波数がどれだけ残ったかを、周波数ごとの数字にします。推定には、短い区間の FFT(信号を周波数成分に分ける計算)を平均する Welch 法を使います。条件はハン窓(区間の端を 0 に絞る重み)、4096 点、50 % オーバーラップ(区間を半分ずつ重ねて取ること)で、理由は基礎編の Welch 法 で説明しています。

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

生の白色ノイズと BPF の出力の PSD に理論曲線を重ねた図

100 Hz の山の頂点は生のノイズの水準に届き、裾は破線に沿って落ちる

シミュレーションと理論の差は平均で約 0.7 dB でした。

分散は何分の 1 になるか

ノイズがどれだけ減るかを 1 つの数字で言うために、分散比を比べます。ENBW(等価雑音帯域幅)は、同じ出力パワーになる理想の長方形フィルタの幅で、分散比に f_s/2 を掛けた値です。

フィルタ 分散比(理論) 分散比(シミュレーション) ENBW
BPF 100 Hz、Q = 5 0.0555 0.0563 27.8 Hz
ノッチ 50 Hz、Q = 10 0.9848 0.9843 492.4 Hz

理論の列は閉じた式で書けます。先ほどの H_{BP} = 1/(1 - jx) の実部は 1/(1 + x^2) で、 |H_{BP}|^2 と同じです。分散比は |H|^2 の全周波数の平均なので、 H_{BP} そのものの平均の実部に等しくなります。

H の全周波数の平均はインパルス応答の先頭 h[0] です。 H(e^{j\theta}) = \sum h[n]e^{-j\theta n} を 1 周平均すると、 n = 0 の項だけが残るからです。 h[0] は分子の b_0 なので、次の式になります。

\sum_n h[n]^2 = b_0 = \frac{\alpha_Q}{1 + \alpha_Q}

f0 = 100 Hz、Q = 5 なら 0.0555 で、表の値です。ノッチは |H_{N}|^2 = 1 - |H_{BP}|^2 なので、分散比は 1/(1 + \alpha_Q) です。f0 = 50 Hz、Q = 10 で 0.9848 になります。バンドパスとノッチも、パワーを分け合っています。

分散比 0.0555 の BPF は、ノイズの大部分を捨てます。ノッチの 0.985 は、ほとんど捨てません。

ノッチは 50 Hz だけを消す

実際の信号では、白色ノイズに 50 Hz のような決まった周波数の妨害が重なることがあります。ノッチがその妨害だけを消せるかを確かめます。50 Hz、振幅 1 の正弦波に σ = 1 の白色ノイズを足し、Q = 10 のノッチに通します。

左は 50 Hz の線がノッチで消える PSD、右は 50 Hz 入力の開始直後に出力が減衰していく時間波形

左で 50 Hz の鋭い線だけが消え、周囲のノイズの水準はそのままであることを見る。右は 50 Hz を入れた直後の過渡で、Q が高いほど消えるまで時間がかかる

過渡が収まったあとの 50 Hz 成分の振幅は、入力の 1.003 に対して出力で 0.00015 でした。全体の分散比は 0.984 なので、ノイズはほとんど減りません。ノッチは「ノイズを減らすフィルタ」ではなく、「決まった周波数の妨害を消すフィルタ」です。

ここまでで、効き目が数字になりました。次の章では、設計の手順にまとめます。


🧪 設計手順と実装

バンドパスとノッチの設計は、#2 ローパス と同じ 4 段で進めます。

  1. 要件:fs、中心周波数 f0、許せる過渡の長さ、周波数のずれの大きさを書き出す
  2. 式: \alpha_Q = \sin\omega_0/(2Q) 、幅は \tan(\Delta\theta/2) = \alpha_Q 、収まる時間は約 4.6/\alpha_Q サンプル
  3. 係数:Cookbook の式を A_0 で割って C の初期化子にする
  4. 検証:ボード線図で幅とずれたときの減衰、ステップ応答で過渡、白色ノイズで分散比を確かめる

設計例 1:50 Hz ハムのノッチ

要件は、fs = 1 kHz、f0 = 50 Hz、落ち着くまで 300 ms 程度まで許す、電源周波数のずれ 0.2 Hz でも −20 dB 以上落とす、とします。

収まる時間の目安 4.6/\alpha_Q = 300 から \alpha_Q = 0.0153 、 Q = \sin\omega_0/(2\alpha_Q) = 0.30902/0.0307 \approx 10 です。Q = 10 なら \alpha_Q = 0.01545 、A0 = 1.01545 で、Cookbook の式を A0 で割ると b0 = 0.9847842、b1 = a1 = −1.8731709、a2 = 0.9695685 です。

検証します。表のとおり 50.2 Hz の減衰は −21.8 dB で要件を満たし、1 % へ収まるのは 296 ms です。

設計例 2:電力計の基本波

電力の基本波(高調波を除いた 50 Hz そのものの成分)を取り出します。高調波は、基本波の整数倍の周波数の成分です。要件は、fs = 1 kHz、第 3 高調波の 150 Hz と第 5 高調波の 250 Hz を −20 dB 以下に落とす、200 ms で落ち着く、とします。Q = 5 を試します。

\alpha_Q = 0.30902/10 = 0.030902 、A0 = 1.030902 なので、b0 = 0.0299754、b1 = 0、b2 = −0.0299754、a1 = −1.8450964、a2 = 0.9400492 です。収まる時間は r = 0.96956 から 149 ms で、中心の群遅延は 1/\alpha_Q = 32.4 ms です。

次の図で検証します。入力は、振幅 1 の基本波に振幅 0.3 の第 3 高調波と振幅 0.15 の第 5 高調波を足した波形です。

左は f0 50 Hz、Q 5 の BPF の振幅特性で 150 Hz が −23.3 dB、250 Hz が −29.8 dB。右は高調波を含む入力と BPF の出力、基本波の時間波形

左の黒丸が高調波の減衰。右では出力が 100 ms ほどで破線の基本波に重なる

150 Hz は −23.3 dB、250 Hz は −29.8 dB で、要件を満たします。出力に残る高調波の振幅は計算で 0.021 と 0.005 です。シミュレーションでは、基本波との差は 107 ms で 0.05 以内に入り、200 ms より後の差は最大 0.025 でした。

周波数のずれにも注意が要ります。50 Hz から 1 Hz ずれると、振幅は −0.18 dB でほとんど変わりませんが、位相は 49 Hz で +11.6°、51 Hz で −11.4° 回ります。電圧と電流に同じフィルタを通せば位相の関係は崩れませんが、フィルタの出力の位相を時刻の基準には使えません。

設計例 3:オーディオ EQ の 1 バンド

Cookbook の本来の用途です。peakingEQ は狙った帯域だけを持ち上げる、または下げる形で、バンドパスやノッチと分母の形が同じです。fs = 48 kHz、f0 = 1 kHz、Q = 1.4、−6 dB を設計します。

Cookbook の式は、ゲインから A = 10^{\mathrm{dBgain}/40} を求め、 \alpha_Q を A で分けて分子と分母に配ります。

\begin{aligned} & b_0 = 1 + \alpha_Q A,\ b_2 = 1 - \alpha_Q A \\ & b_1 = -2\cos\omega_0 \end{aligned} \begin{aligned} & A_0 = 1 + \frac{\alpha_Q}{A},\ A_2 = 1 - \frac{\alpha_Q}{A} \\ & A_1 = -2\cos\omega_0 \end{aligned}

\theta = \omega_0 では \cos\theta - \cos\omega_0 の項が消えるので、残るのは分子の \alpha_Q A と分母の \alpha_Q/A の比で、ゲインは A^2 です。dB にすると 40\log_{10}A で、狙いの値に戻ります。

数字を入れると、 \omega_0 = 0.13090 、 \alpha_Q = 0.046616 、A = 0.70795 で、A0 で割ると b0 = 0.9691836、b1 = a1 = −1.8603877、b2 = 0.9072574、a2 = 0.8764410 です。f0 のゲインは計算で −6.00 dB になりました。

これを帯域の数だけ直列につなぐと EQ になります。次の図は、100 Hz を +4 dB、1 kHz を −6 dB、8 kHz を +3 dB にした 3 段の例です。

100 Hz、1 kHz、8 kHz の peakingEQ の振幅特性と、3 段を直列にした合計の振幅特性

細い線が各段、太い線が合計。直列につなぐと dB は足し算になるので、帯域ごとに独立して調整できる

C のコードに入れる値

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

// BPF(f0 = 100 Hz, Q = 5, fs = 1 kHz)
static biquad_t bpf100 = {
    0.0555154112f, 0.0f, -0.0555154112f,
    -1.5282081665f, 0.8889691775f, 0.0f, 0.0f
};

// ノッチ(f0 = 50 Hz, Q = 10, fs = 1 kHz)
static biquad_t notch50 = {
    0.9847842466f, -1.8731709497f, 0.9847842466f,
    -1.8731709497f, 0.9695684932f, 0.0f, 0.0f
};

// BPF(f0 = 50 Hz, Q = 5, fs = 1 kHz)電力計の基本波
static biquad_t bpf50 = {
    0.0299754084f, 0.0f, -0.0299754084f,
    -1.8450964177f, 0.9400491832f, 0.0f, 0.0f
};

// peakingEQ(f0 = 1 kHz, Q = 1.4, -6 dB, fs = 48 kHz)
static biquad_t eq1k = {
    0.9691835856f, -1.8603877431f, 0.9072573700f,
    -1.8603877431f, 0.8764409556f, 0.0f, 0.0f
};

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


💡 活用事例

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

用途 フィルタ 設計値の例
50 Hz ハム ノッチ、移動平均 Q = 10、または N = 20
オーディオ EQ biquad の直列 1 バンドにつき biquad 1 段
電力計の基本波 BPF f0 = 50 Hz、Q = 5

50 Hz ハムの除去。 配線に乗った商用電源のノイズです。fs = 1 kHz、Q = 10 のノッチで、幅 4.9 Hz だけ削れます。同じ fs で N = 20 の移動平均(#2)でも 50 Hz と高調波(整数倍の周波数の成分)に谷が入るので、高調波も問題なら移動平均のほうが向いています。

オーディオの EQ。 EQ(イコライザ)は、帯域ごとに音量を上げ下げする機能です。Audio EQ Cookbook は名前のとおりこの用途の式集です。fs = 48 kHz で、peakingEQ(狙った帯域だけを持ち上げる、または下げる形)の biquad を帯域の数だけ直列につなぎます。1 バンドは biquad 1 段で、係数の計算式はこの記事のノッチやバンドパスと同じ形です(設計例 3)。

電力計の基本波の抽出。 基本波は、高調波を除いた 50 Hz そのものの成分です。その基本波だけの実効値を見たいとき、RMS の計算の前に BPF を入れます。f0 = 50 Hz、Q = 5 なら高調波を落とせます。電圧と電流の両方に同じフィルタを通せば、位相の関係も崩れません。全高調波を含む真の実効値を測るなら、逆にこの BPF は入れません(設計例 2)。

どの場面でも、狙う周波数の幅と、許せる時間の長さから Q を決めています。次の章で、うまくいかないときの原因をまとめます。


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

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

症状 よくある原因 対策
ノッチが浅い f0 のずれに対して Q が高い Q を下げる
消えるまで時間がかかる Q が高い Q を下げる。収まる時間は約 4.6/α_Q サンプル
帯域幅が計算と合わない Cookbook と iirnotch の Q を混同 tan(Δθ/2) = α_Q で計算し直す
中心の周波数がずれている 呼び出し周期が設計の fs と違う 実際の呼び出し周期で計算し直す
出力が発散する 極が単位円の外。係数の符号違いや丸め 係数から極を計算する。a1, a2 の符号を見直す
起動直後に大きく振れる 状態変数が 0 から始まる 最初の値で biquad_reset() する

✅ まとめ

  • Cookbook の係数は、アナログのプロトタイプの双一次変換から導ける。f0 と Q を入れるだけで決まり、分母は BPF とノッチで共通
  • 周波数応答は 1/(1 - jx) の 1 行に整理でき、帯域幅、群遅延、分散比が \alpha_Q の式で書ける
  • Q は谷の幅と過渡の長さのトレードオフ。収まる時間は約 4.6/\alpha_Q サンプル
  • BPF とノッチは足すと 1 で、ナイキスト線図で中心 0.5、半径 0.5 の同じ円を描く
  • ノッチはノイズを減らすフィルタではなく、決まった周波数の妨害を消すフィルタ

ここまでで、ローパス、ハイパス、バンドパス、ノッチの 4 種類を要件から設計できるようになりました。手元の ADC の値を FFT で眺めて、50 Hz の線が立っていたら Q = 10 のノッチを入れてみるのが、次の一手です。

ここまでのフィルタは、信号の中身を知らないまま周波数だけで切り分けます。ゆっくりした信号の周波数帯にノイズが重なっていると、線形フィルタ(この連載のフィルタのように、入力を 2 倍すると出力も 2 倍になるもの)では分けられません。次回は、対象の動きのモデルを持ち込んで推定するカルマンフィルタへ進みます。


よくある質問(FAQ)

Q. ノッチの Q はいくつにすればいいですか?

A. 消したい周波数がどれだけ揺れるかで決めます。この記事の例では、Q = 10 なら 0.2 Hz ずれても −21.8 dB 落ち、落ち着くまでは約 300 ms です。Q = 30 だと谷の幅は 1.6 Hz と狭いですが、同じずれで −12.5 dB にとどまり、落ち着くまで 900 ms ほどかかります。まず Q = 5〜10 あたりから試すのが無難です。

Q. 50 Hz を消すのはノッチと移動平均のどちらがいいですか?

A. 消したいのが 50 Hz だけならノッチです。fs = 1 kHz、Q = 10 で幅 4.9 Hz だけを削り、周辺の信号はほとんど残ります。高調波も問題なら、N = 20 の移動平均が 50 Hz とその整数倍に谷を置けます。代わりに信号の帯域も削れます。

Q. バンドパスの Q を上げるとノイズはどれだけ減りますか?

A. 分散比は \alpha_Q/(1 + \alpha_Q) で、Q が大きいところではほぼ Q に反比例します。f0 = 100 Hz では Q = 2、5、15 で 0.128、0.056、0.019 です。代わりに中心の群遅延は 6.8 ms、17.0 ms、51.0 ms と Q に比例して伸びます。

Q. ノッチの周波数がずれるときはどうしますか?

A. Q を下げて谷の幅を広げます。Q = 2 なら 0.2 Hz ずれても −35.8 dB 落ちます。代わりに周辺の信号も削れます。周波数の揺れが大きいなら、固定のノッチではなく、周波数を追いかける適応フィルタが候補で、連載の #8 で扱う予定です。


関連記事

この連載:

センサの入口:


参考

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

理論

係数の式

設計ツール(SciPy v1.18.0 Manual)