デジタルフィルタ入門 #1｜基礎編：サンプリング定理から伝達関数・ボード線図・C の biquad まで

はじめに

ADC の値をシリアルプロッタに流すと、たいていガタガタしています。そこで「とりあえず 10 回平均」を入れると線はきれいになります。ところが、同じ値をモータの速度制御に戻すと、今度はループが揺れ始める。電子工作でよくある流れです。

きれいになった分だけ、どこかで何かを払っています。この連載のゴールは、それを勘ではなく数字で言えるようになることです。

センサの値がガタガタしているとき、どのフィルタをどう設計すれば、どこまで落ち着き、何を失うのか。 それを数式と図で説明でき、自分のマイコンのコードに落とせるところまで進みます。

第 1 回の基礎編では、フィルタを 1 つずつ作る前に、そのための道具をそろえます。ADC が数列を作るときの決まり、式から周波数ごとの倍率と遅れを取り出す手順、ノイズで効き目を測る方法、そしてマイコンで動かす C の骨格です。#2 からはこの道具で、基本の 4 種類を 1 本に 1 つずつ設計します。ローパスは低い周波数を通し、ハイパスは高い周波数を通し、バンドパスは狙った範囲の周波数だけを通し、ノッチは狙った周波数だけを消すフィルタです。実機は使わず、すべて計算と数値シミュレーションで示します。図はサンプリング周波数 1 kHz にそろえたので、フィルタどうしを同じ土俵で比べられます。

この記事で分かること： サンプリングが信号の何を決めてしまうのか、差分方程式がなぜ 1 つの分数 H(z) にまとまるのか、その分数から振幅、位相、遅れの時間をどう取り出すのか、白色ノイズで効き目を測るときに何を計算しているのか、が分かります。

読み終えたときにできること： コードに書かれた差分方程式から H(z) を書き下し、狙いの周波数での倍率と遅れを手で計算できます。fs を見てエイリアシングの行き先を計算でき、フィルタの出力の PSD を Welch 法で求めて理論と比べられます。設計した係数を C の biquad に入れて動かせます。

🧭 まず全体像

デジタルフィルタとは、サンプリング（アナログの信号を一定の間隔で読み取ること）された数列を別の数列に変える計算のことです。抵抗やコンデンサの代わりに、足し算・掛け算と「1 サンプル前の値を覚えておくメモリ」でできています。信号の通り道の中では、次の位置に入ります。

図のアンチエイリアスは、ADC の前で高い周波数を落とすアナログのフィルタです。必要な理由は、この後の 2 つの節で説明します。

この先で繰り返し使う波の言葉を、先にそろえておきます。

言葉がそろったので、信号の通り道の入口である ADC から見ていきます。

サンプリング定理とナイキスト周波数

ADC が値を読む間隔をどこまで細かくすれば、元の信号の情報を失わずに済むのか。この小節で、その答えと、デジタルフィルタが扱える周波数の上限を確かめます。

信号に含まれる最高周波数を f_{\max} とします。サンプリング周波数 f_s （1 秒間に値を取り込む回数）が f_{\max} の 2 倍より高ければ、とびとびのサンプル列から元の連続信号（切れ目のないアナログの信号）を一意に復元できます。これがサンプリング定理（標本化定理・Nyquist–Shannon の定理）です。

裏返すと、 f_s で取り込んだ信号が持てる周波数の情報は f_s/2 までです。この上限をナイキスト周波数と呼びます。

次の小節では、この定理の条件を破ると何が起きるかを見ます。

エイリアシング（折り返し）

では、ナイキスト周波数より高い成分が ADC に入るとどうなるか。消えるのではなく、低い周波数のふりをして混ざります。これがエイリアシング（折り返し）です。この小節では、化ける仕組みと、なぜ ADC の前で防ぐしかないのかを確かめます。

なぜ区別がつかなくなるのか。ADC が見ているのは、サンプリング周期 T = 1/f_s （サンプルの時間間隔）ごとの点 t = nT の値だけです。周波数 f の余弦波に f_s の整数倍を足しても、点の上の値は 1 つも変わりません。

k は整数です。 f_s T = 1 なので k f_s \cdot nT = kn は整数になり、角度は 2\pi の整数倍だけ増えるので余弦の値は元に戻ります。余弦は偶関数（ f(-x) = f(x) を満たす関数）なので、 -f の波も同じ数列です。つまり f 、 f + k f_s 、 k f_s - f の波は、サンプルの上では同じ数列になります。ADC は点しか見ていないので、区別する手がかりがありません。

「折り返し」という名前は、周波数軸を f_s/2 のところで紙のように折る操作から来ています。折ると、 f_s/2 より上の成分 f は f_s - f の位置に重なります。fs = 1 kHz なら、970 Hz は 30 Hz に、630 Hz は 370 Hz に重なります（正弦波では符号が反転します）。どの周波数でも、0〜 f_s/2 に写る先は、 f を f_s で割った余りを取り、余りが f_s/2 を超えていれば f_s から引くと求まります。

次の図の 970 Hz の正弦波は、定理の条件を破った例です。復元には 1,940 Hz より高い f_s が要るのに、1 kHz でサンプリングしています。

同じ信号を、定理の条件を満たす fs = 2 kHz でサンプリングすると次のようになります。図の破線は sinc 補間の結果です。sinc 補間は、各サンプルの位置に sinc 関数（ \sin x / x の形の波）を置いて足し合わせ、点の間の値を求める計算です。

同じ 970 Hz でも、fs = 2 kHz なら f_s > 2 f_{\max} を満たすので、点は元の波をなぞり、点の列から 970 Hz の波を復元できます。定理の条件を満たすか破るかの違いが、この 2 枚の差です。

いくつかの周波数について、どこに見えるかを上の求め方で計算すると次のようになります。970 Hz と 630 Hz は、上の 2 枚の図と同じ例です。

この表の範囲では、計算は 2 通りです。 f が f_s/2 を超えていれば f_s - f 、 f_s を超えていれば f - f_s です。

1 kHz の列では、30 Hz、970 Hz、1030 Hz が同じ 30 Hz に重なります。ADC のあとでは、この 3 つを区別できません。1030 Hz は 2 kHz でも f_s/2 を超えるので、970 Hz に化けます。定理が効くのは、入ってくる成分が本当に f_s/2 より下にあるときだけです。

周波数が f_s の整数倍に一致すると、行き先は 0 Hz です。50 Hz のハム（商用電源から配線に乗るノイズ）を 50 Hz でサンプリングすると毎回同じ位相を読むので、一定の値（DC）に見えます。

同じ現象は映像でも起きています。カメラのコマ数がサンプリング周波数に当たります。動画の車輪が逆回転して見えたり、回っているプロペラが止まって見えたりするのは、回転に対してコマ数が足りず、実際とは違う動きがコマの列に現れるためです。点滅する光の下で回転物が止まって見えるストロボ効果も、同じ原理です。

測定では、この化け方が害になります。fs = 1 kHz で 900 Hz のノイズが入ると、ADC のあとでは 100 Hz の成分として残り、信号と見分けがつきません。デジタルフィルタで 100 Hz を落とせば、そこにあった本物の 100 Hz まで一緒に失います。ノイズは、折り返す前に落とすしかありません。

ここまでが ADC の入口の話です。次の章からは、ADC を通ったあとの数列を扱う道具をそろえます。

📐 道具立て（制御工学の言葉で）

フィルタの効き目と遅れを数字で言うには、コードの式から周波数ごとの倍率と遅れの時間を取り出す道具が要ります。ここからの節は、すべて同じ 5 段の階段を上ります。差分方程式、z 変換、伝達関数、周波数応答、そして振幅・位相・群遅延です。

連続時間から離散時間へ

アナログ回路の信号は、時刻 t が切れ目なく流れる連続時間の関数 x(t) で、ふるまいは微分方程式（時間あたりの変化率を含む式）で書きます。マイコンが持てるのは、一定の間隔で読んだ値の列だけです。時刻を番号 n で指すこの世界が離散時間で、ふるまいは差分方程式（今の値を過去のサンプルから計算する式）で書きます。階段の 1 段目の前に、この 2 つをつなぎます。

サンプルとサンプルの時間間隔をサンプリング周期と呼び、 T = 1/f_s と書きます。連続時間の微小な刻み \Delta t に当たる量です。時刻を t \to nT と置き換えると、信号は数列になります。

角括弧の [n] は「 n 番目のサンプル」の意味で、この記事の式はすべてこの数列の式です。

次は微分です。微分は接線の傾きですが、数列には接線がありません。そこで、隣のサンプルとの差を T で割って近似します。これが差分で、1 サンプル前との差をとる形を後退差分と呼びます。

積分も同じように、幅 T の長方形の面積の和に置き換わります。

微分が差分に、積分が和に置き換わるので、微分方程式は差分方程式になります。

ここまでの置き換えを、電子工作でいちばん身近なアナログフィルタで試します。抵抗 R とコンデンサ C を 1 個ずつ使う RC ローパスです。狙いは、アナログ回路の式がマイコンの差分方程式にそのまま変わることを確かめることです。ここで出てくる式が、#2 で扱う 1 次 IIR の正体です。

入力電圧を x 、コンデンサの電圧（出力電圧）を y とします。R の両端の電圧は x - y なので、オームの法則から R を流れる電流が決まります。

出力から電流を取らないとすると、この電流はすべて C に入ります。コンデンサに流れ込む電流は C\,dy/dt なので、次の式が成り立ちます。

両辺に R を掛けて、 y を左辺に移します。

この積 RC を時定数と呼び、 \tau = RC と書きます。時定数は、入力を急に変えたとき、出力が最終値の約 63 %（ 1 - e^{-1} ）まで追いつくのにかかる時間です。5τ でほぼ（99 %）追いつきます。τ が大きいほど出力はゆっくり動き、カットオフは低くなります。これを使うと、回路の微分方程式は次の形です。

カットオフ周波数は f_c = 1/(2\pi RC) = 1/(2\pi\tau) で、RC フィルタ計算ツール と同じ式です。

微分を後退差分に置き換えます。

両辺に T を掛けて、 y[n] についてまとめます。

両辺を \tau + T で割って整理すると、次の形になります。

微分方程式が、過去の値から今の値を計算する 1 行の式になりました。これは #2 で扱う 1 次 IIR（1 つ前の出力を戻して計算するフィルタ。指数移動平均とも呼ぶ）と同じ形です。

数字を入れてみます。fc = 200 Hz の RC は \tau = 1/(2\pi \cdot 200) = 0.796 ms で、fs = 1 kHz（T = 1 ms）なら \alpha = 1/(0.796 + 1) = 0.557 です。T が τ より長いので、差分の近似は粗くなります。

そこで、微分を近似する代わりに、1 サンプルでの減り方そのものを RC に合わせる決め方もあります。入力が 0 のとき、RC の出力は e^{-t/\tau} の指数関数で減るので、1 サンプル T のあいだに e^{-T/\tau} 倍になります。一方、上の 1 行の式は、入力が 0 なら y[n] = (1-\alpha)\,y[n-1] で、1 サンプルごとに 1-\alpha 倍です。この 2 つの倍率を等しいと置くと \alpha = 1 - e^{-T/\tau} = 0.715 で、同じ回路でも値が変わります。 \alpha の決め方は 1 通りではないので、#2 で比べます。

インパルスも置き換えます。インパルスは、ごく短い時間だけ立ち上がる鋭いパルスです。連続時間では \delta(t) （デルタ関数）と書き、これに対応するのが単位サンプル \delta[n] で、 n = 0 のときだけ 1、それ以外は 0 の数列です。 \delta[n] を入れたときの出力をインパルス応答と呼びます。この後の FIR、IIR の名前も、この出力の話です。

時間領域での対応をまとめると次のようになります。

アナログの式をマイコンの式へ移す対応がそろいました。次は、マイコンで書くフィルタの式を 1 つの形にまとめます。

差分方程式：マイコンで書くそのままの式

マイコンの中のフィルタは、次の形に必ず書けます。今の出力を、今と過去の入力、過去の出力の重み付き和（それぞれに係数を掛けて足したもの）で作る式です。

b_k が入力側（フィードフォワード）、 a_k が出力の戻し（フィードバック）の係数です。設計とは、この係数を決める作業にほかなりません。ただ、この形のままでは、係数と周波数の特性の関係が見えません。

z 変換と伝達関数

差分方程式には遅れの項が混ざっているので、係数を眺めても周波数の特性は読み取れません。そこで、遅れを掛け算で書ける形へ式を移します。その道具が z 変換で、この小節の到達点は、フィルタを 1 つの分数 H(z) で表すことです。

z 変換は、数列を多項式のように扱う道具です。数列の各要素に z^{-n} の重みを掛けて足し合わせます。信号は n = 0 から始まり、それより前は 0 とします。

書き下すと x[0] + x[1]z^{-1} + x[2]z^{-2} + \cdots です。z は「1 サンプル後ろ」を数える目印で、 z^{-n} が付いた項が n 番目のサンプルを表します。

この定義から、フィルタの計算でいちばん使う性質が出ます。1 サンプル遅らせた数列 x[n-1] を z 変換してみます。

m = n-1 と置き換えると、次のようになります。 x[-1] = 0 なので、和は m = 0 から始まります。

1 サンプル遅らせることは、 z^{-1} を掛けることと同じです。 2 サンプル遅らせるなら z^{-2} を掛けます。

これを使うと、差分方程式を掛け算の式に変えられます。出発点は次の 2 次の差分方程式です。「差分方程式」の節の一般式を M = L = 2 にした形です。

z 変換は線形です。線形とは、定数倍は定数倍のまま、和は和のまま移る性質のことで、両辺を項ごとに変換できます。遅れている項には先ほどの性質を使います。

Y(z) の項を左辺に集めて、両辺をくくります。

両辺を割ると、出力と入力の比が取れます。この比が伝達関数です。

差分方程式の係数が、そのまま分子（分数の上の式）と分母（分数の下の式）に並びます。分子が入力側の b_k 、分母が出力側の a_k です。設計が係数を決める作業であることが、式の形からそのまま見えます。式に出てくる z^{-1} の最大のべきを次数と呼び、この式は 2 次です。

この 2 次の形を biquad（バイクアッド）と呼びます。分子も分母も 2 次なので bi-quadratic です。この連載の #2〜#4 の IIR フィルタは、1 次 IIR を除いてすべて biquad 1 段で書けます。

同じ手順を、RC ローパスから導いた式で確かめます。

括弧を開いて、今の入力の項と 1 つ前の出力の項にまとめます。

両辺を z 変換します。

Y(z) を左辺に集めて、比を取ります。

分母を 0 にする z を極と呼びます（「極と零点」の節で説明します）。この式の極は z = 1-\alpha で、α が小さいほど極は 1 に近づきます。

アナログの回路や制御では、微分方程式をラプラス変換（時間の関数を複素数 s の式に変える、z 変換のアナログ版）して伝達関数 G(s) を書きます。これを知っている方は、次の対応で読み替えられます。先ほどの表が時間領域の対応で、こちらは変換したあとの領域の対応です。

表の言葉を先に確認します。極は伝達関数の分母を 0 にする点で、「極と零点」の節で説明します。安定は、入力を止めれば出力が収まり、発散（際限なく大きくなること）しない状態です。左半面は複素数 s を描いた平面の左半分（実部が負の領域）、単位円は半径 1 の円です。1 次遅れは、RC ローパスのように時定数 τ で指数関数的に追いつく要素です。

フィルタを表す H(z) が手に入りました。次の小節で、この分数から周波数ごとの倍率を取り出します。

周波数応答：z を単位円の上に置く

知りたいのは、ある周波数の成分が何倍になり、どれだけずれるかです。H(z) からこの 2 つを取り出すのが周波数応答です。式に複素数が出てくるので、言葉を先にそろえます。

入力に、複素数で書いた正弦波 x[n] = e^{j\omega nT} を入れて考えます。この数列は 1 サンプル遅らせると全体が e^{-j\omega T} 倍になります。つまり、遅れを表す z^{-1} が e^{-j\omega T} という数に置き換わり、H(z) がそのまま出力と入力の比になります。

伝達関数の z に e^{j\omega T} を入れると、周波数ごとの応答が得られます。その周波数の正弦波を入れたときに、振幅と位相がどう変わるかを表す量です。 T = 1/f_s はサンプリング周期、 \omega = 2\pi f です。複素数は絶対値と偏角に分けて書けるので、次の形に整理します。

絶対値 |H| がその周波数を何倍にするか（振幅特性）、偏角 \angle H がどれだけ位相をずらすか（位相特性）を表します。

振幅はこれで読めます。位相は角度のままでは扱いにくいので、次の小節で時間に直します。

群遅延：「何ミリ秒遅れるか」の正体

位相は角度なので、そのままでは「何ミリ秒遅れるか」の感覚がつかめません。そこで周波数で微分して、時間の次元に直します。

これが群遅延で、その周波数付近の信号の包絡線（細かい波の山をなぞった外形の線）が何秒遅れて出てくるかを表します。位相が周波数に比例して直線なら、群遅延は一定です。全周波数が同じ時間だけ遅れるので、波形の形が崩れません。これを直線位相と呼びます。

振幅と位相の読み方がそろいました。次は、1 点ずつ計算しなくても特性の形が分かる図を導入します。

極と零点：z 平面で周波数特性を読む

周波数応答を 1 点ずつ計算しなくても、H(z) の分母と分子を 0 にする点の位置から、グラフの形を読めます。この小節で、その読み方を身につけます。

H(z) の分子を 0 にする点が零点、分母を 0 にする点が極です。複素数 z を描いた平面（z 平面）で単位円の上を歩くことが、周波数を 0 から f_s/2 へ上げることに対応します。 |H| の大きさは、歩いている点から各零点までの距離の積を、各極までの距離の積で割ったものに比例します。z 平面に極と零点を打った図を、極零点図と呼びます。

極が単位円の内側なら安定、外に出ると発散します。 これは s 平面の左半面に対応する条件です。

極があるかどうかは、フィルタを 2 種類に分ける基準にもなります。

FIR と IIR

過去の出力を計算に戻すかどうかで、フィルタの性質は大きく変わります。過去の出力を戻さない（ a_k = 0 ）のが FIR、戻すのが IIR です。インパルスを入れたときの応答が、有限の長さで終わるか、理論上は無限に続くかが名前の由来です。FIR の係数の個数をタップ数と呼びます。次の表の急峻さ（切れ味）は、通す帯域から落とす帯域へ移るときの振幅の落ち方の急さです。量子化は、値を限られたビット数に丸めることで、ここでは係数や途中の計算値の丸めを指します。

5 段の階段を 1 本で通す：2 点平均

道具が全部そろったので、いちばん小さなフィルタで階段を上から下まで通します。題材は 2 点平均、つまり今の値と 1 つ前の値を足して 2 で割る計算です。手で全部計算できる大きさなので、式の 1 つ 1 つが何をしているかを確かめられます。fs = 1 kHz、 T = 1 ms とします。

1 段目は差分方程式です。 b_0 = b_1 = 1/2 で、ほかの係数は 0 です。

2 段目と 3 段目で z 変換して比を取ると、伝達関数が出ます。分子と分母に z を掛けた右の形にすると、零点と極が読めます。

零点は分子を 0 にする z = -1 、極は分母を 0 にする z = 0 です。単位円の上の z = -1 は fs/2 に当たるので、500 Hz が完全に消えることがここで分かります。

4 段目は z = e^{j\theta} （ \theta = \omega T ）を入れる周波数応答です。指数の半分の e^{-j\theta/2} をくくり出すと、オイラーの式から出る e^{j\theta/2} + e^{-j\theta/2} = 2\cos(\theta/2) が使えます。

5 段目で振幅、位相、群遅延に分けます。 e^{-j\theta/2} は絶対値が 1 なので振幅には効かず、位相だけを決めます。

位相が ω に比例する直線なので、群遅延は周波数によらず一定です。

数字を入れます。f = 250 Hz は \theta = \pi/2 なので、 |H| = \cos(\pi/4) = 0.707 、つまり −3.0 dB で、位相は −45° です。2 点平均のカットオフは fs/4 だと分かります。

次の図の左がこの計算のボード線図です。右の極零点図でも同じ値を読みます。250 Hz は z = j の点で、零点までの距離は \sqrt{2} 、極までの距離は 1 です。距離の比に分子の定数 1/2 を掛けると、 \sqrt{2}/2 = 0.707 で同じ値になります。

もう 1 つ、分母のある biquad でも同じ計算をします。#2 で設計する 20 Hz の 2 次バターワース（b0 = 0.0036217、b1 = 0.0072434、b2 = 0.0036217、a1 = −1.8226949、a2 = 0.8371817）が、100 Hz をどこまで落とすかを求めます。 \theta = 2\pi \cdot 100/1000 なので、 z = e^{j\theta} = 0.80902 + 0.58779j です。これを H(z) の分子と分母に入れて絶対値を取ると、分子は 0.013103、分母は 0.34973 になります。

dB にすると 20\log_{10} 0.03747 = -28.5 dB です。この値は、実装の章の Python が周波数特性から出力する値と同じです。どんなに複雑なフィルタでも、計算の手順はこの 1 本道です。

一周したところで、次の章では、1 つの周波数ではなくノイズ全体に対する効き目を測る道具を用意します。

📏 効き目を測る道具

ここまでは 1 つの周波数ごとの話でした。実際のセンサノイズは、あらゆる周波数が混ざった信号です。#2〜#4 では、その代表である白色ノイズをフィルタに通して効き目を測ります。この章では、そのとき使う言葉と計算の手順をまとめます。結果の図は、フィルタごとに各記事で見ます。

まず、ばらつきを数字で言うための言葉をそろえます。

白色ノイズと PSD

白色ノイズがフィルタを通るとパワーの分布がどう変わるかは、式 1 本で書けます。この連載の計算の土台になる式です。シミュレーションの条件は、全記事で次のものを使います。

条件は fs = 1 kHz、N = 2^16 点、σ = 1 のガウス白色ノイズです（乱数列を決める初期値のシードは固定）。

理論は単純です。片側 PSD は、パワーを正の周波数だけにまとめた PSD です。 \sigma^2 のパワーが 0〜 f_s/2 に均等に広がるので、入力の片側 PSD は S_x = 2\sigma^2/f_s です。フィルタは周波数ごとに振幅を |H| 倍するので、パワーは |H|^2 倍になります。

出力の分散は、この PSD を積分した値です。パーセバルの定理は、信号のエネルギーを時間の側で足しても周波数の側で足しても同じ値になる、という定理です。これを使うと、 |H|^2 を全周波数で平均した値はインパルス応答 h[n] の 2 乗和に等しくなるので、分散は次の式でも求まります。

最後の式の出どころを確かめます。全周波数の平均を取るには \theta = \omega T を -\pi から \pi まで積分して 2\pi で割ります。周波数応答は H(e^{j\theta}) = \sum_n h[n] e^{-j\theta n} なので、絶対値の 2 乗は 2 重の和です。

e^{-j\theta(n-m)} を 1 周ぶん平均すると、 n = m の項だけが 1 で残り、それ以外は 0 になります。正と負の部分が打ち消し合うからです。残るのは \sum_n h[n]^2 で、これが |H|^2 の全周波数の平均です。

式の意味は単純です。 |H|^2 の平均が大きいフィルタほど、ノイズがたくさん通ります。その平均は、インパルス応答を 2 乗して足すだけで求まります。

分散比と ENBW

ノイズがどれだけ減るかを 1 つの数字で言うために、出力の分散を入力の分散で割った分散比を使います。白色ノイズなら、前の小節の式から \sum h[n]^2 がそのまま分散比です。

帯域の広さを 1 つの数字で比べるのが等価雑音帯域幅（ENBW）です。同じ出力パワーになる理想の長方形フィルタの幅を表します。全帯域の幅 f_s/2 に分散比を掛けて、次の式で求まります（ピークのゲインが 1 のとき）。

2 点平均で試します。インパルス応答は 1/2 が 2 個なので、 \sum h[n]^2 = 1/4 + 1/4 = 1/2 です。周波数の側から計算しても、 |H|^2 = \cos^2(\theta/2) の平均は 1/2 で一致します。ENBW は 500 Hz × 1/2 = 250 Hz です。同じ白色ノイズを 2 点平均に通すシミュレーションでは、分散比は 0.501 でした。

分散比が小さいほどノイズは減ります。ただし分散比だけを見ると遅れの代償が見えないので、#2 からは必ず群遅延と並べて読みます。

ピリオドグラムと Welch 法

分散比は 1 つの数字でした。どの周波数がどれだけ残ったかを見るには、シミュレーションの出力から PSD を推定します。この小節では、その推定の手順と、理論と重ねるときの条件を決めます。

いちばん素朴な方法は、信号全体を 1 回 FFT（高速フーリエ変換：信号を周波数成分に分ける計算）して、周波数ごとのパワーを取ることです。この結果をピリオドグラムと呼びます。ところが、ノイズのピリオドグラムはギザギザが激しく、点数を増やしても滑らかになりません。点数を増やしても周波数の刻みが細かくなるだけで、1 点ごとのばらつきは減らないからです。σ = 1 の白色ノイズで計算すると、真の PSD に対する 1 点ごとの標準偏差は、約 3,800 点でも約 61,000 点でも 100 % 前後のままでした。

そこで Welch 法を使います。1 回の長い FFT の代わりに、短い FFT を何度も行って平均する方法です。手順は次の 3 つです。

1. 信号を 4096 点の区間に分ける
2. 各区間に窓を掛けて FFT し、パワーを取る
3. 全区間で平均する

独立な K 区間を平均すれば、ばらつきの標準偏差は 1/\sqrt{K} 倍になります。代わりに周波数の刻みは区間の長さで決まり、ここでは 1 kHz ÷ 4096 = 0.24 Hz です。この連載の用途には十分な細かさです。

手順 2 の窓は、区間の切り口の問題を抑えるためのものです。信号を区間に切り出すと、端で信号が急に切れます。この不連続が、元の信号にない周波数成分を作ります。これがスペクトル漏れ（漏れ）です。そこで、区間の端で 0 になる滑らかな重みを信号に掛けます。この重みが窓関数（窓）です。何も掛けずに切り出すのは、区間の中を 1 倍する四角い窓を掛けることと同じなので、矩形窓と呼びます。ハン窓は cos 形の山の窓で、両端が 0、中央が 1 になります。よく使われる標準的な窓で、区間の点数を N とすると次の式で表されます。

次の図の左で 2 つの窓の形を、右で切り出したあとのスペクトルを比べます。右の入力は 100.22 Hz の正弦波 1 本だけなので、本来なら線は 1 本で、裾は出ないはずです。

正弦波から 10 Hz 離れた位置の値は、矩形窓の約 −42 dB に対してハン窓は約 −106 dB です。この裾の差が、PSD の深い谷を測るときに効きます。矩形窓のままだと、近くの大きな成分から漏れた裾が谷を埋めてしまいます。

窓にも代償があります。区間の端を弱めるので、そのまま区間を並べると、境目付近のデータがほとんど使われません。そこで区間を半分ずつ重ねて取ります。これが 50 % オーバーラップです。先頭の 4096 点を除いた 61,440 点で比べると（#2 以降でフィルタの出力が落ち着くまでの区間を捨てるので、同じ範囲にそろえています）、重ねなければ 15 区間で標準偏差は 25 %、半分ずつ重ねると 29 区間で 19 % になりました。

ここまでの条件（ハン窓、4096 点、50 % オーバーラップ）で、生の白色ノイズの PSD を求めます。σ = 1 の水準の理論値は 10\log_{10}(2/f_s) = -27.0 dB です。次の図は、灰色が 1 回の FFT、青が Welch 法、破線が理論の水準です。見るのは、青の線が破線のまわりに集まるかどうかです。

#2〜#4 では、同じ条件でフィルタの出力の PSD を求め、 S_y(f) の式で計算した理論値と重ねます。ここまでで測り方がそろいました。次の章で、フィルタをコードにする骨格を作ります。

🛠️ 実装の骨格

#2〜#4 で設計する係数を動かすコードを、ここでまとめて書きます。各記事では、このコードに入れる係数だけを示します。浮動小数点の biquad、固定小数点の 1 次 IIR、Python からの係数の持ち出し、既存のライブラリ、FIR の順に見ていきます。

C の biquad（Direct Form II Transposed）

まず、この連載の IIR フィルタをまとめて動かせる biquad を C で書きます。biquad の実装にはいくつかの形があります。伝達関数は同じでも、計算の順番と状態変数の持ち方が違います。Direct Form は差分方程式をそのまま信号の流れ図にした形で、II はそれを組み替えて遅延メモリを共有し、状態変数を 2 個に減らしたものです。Transposed（転置）は、その流れ図の矢印を逆向きにした形を指します。浮動小数点ならこの Direct Form II Transposed（DF2T）が定番で、状態変数は 2 個で済みます。零点の計算が極より先に来るので、極が作る大きなゲインで途中の値が膨らみにくい、という利点もあります。係数の符号は SciPy と同じです。

typedef struct {
    float b0, b1, b2, a1, a2;  // a0 = 1 に正規化した係数
    float z1, z2;              // 状態変数
} biquad_t;

static inline float biquad_step(biquad_t *f, float x)
{
    float y = f->b0 * x + f->z1;
    f->z1 = f->b1 * x - f->a1 * y + f->z2;
    f->z2 = f->b2 * x - f->a2 * y;
    return y;
}

// 入力 x0 が十分長く続いた状態から始める（起動直後の過渡を消す）
static inline void biquad_reset(biquad_t *f, float x0)
{
    float y = x0 * (f->b0 + f->b1 + f->b2) / (1.0f + f->a1 + f->a2);
    f->z2 = f->b2 * x0 - f->a2 * y;
    f->z1 = f->b1 * x0 - f->a1 * y + f->z2;
}

// 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
};

例の係数は、#2 で設計する 2 次バターワース LPF です。タイマ割り込みなどで、一定の周期で biquad_step() を呼びます。係数は fs とセットなので、呼ぶ周期がぶれるとカットオフも一緒にぶれます。delay() で作った周期は処理時間の分だけ伸びることに注意してください。

このコードが本当に H(z) を計算しているかを、式で確かめます。z1 と z2 は次のサンプルへ持ち越す値で、z 変換すると z^{-1} が付きます。コードの 3 行は次の 3 本の式です。

下から順に代入して Z_1 と Z_2 を消すと、次の式になります。

Y を左辺に集めれば、「z 変換と伝達関数」で導いた H(z) そのものです。

biquad_reset() は、入力 x_0 が十分長く続いたあとの状態を計算しています。過渡（入力を変えた直後の、出力が動いている間）が終わると出力は y = H(1)\,x_0 で一定になるので、これを z1 と z2 の式に入れれば、最初の 1 サンプルから落ち着いた状態で始められます。 H(1) は z = 1 、つまり DC での倍率です。

浮動小数点の計算が遅いマイコンでは、次の小節の整数の書き方が役に立ちます。

固定小数点の 1 次 IIR

固定小数点は、小数点の位置を決めておき、整数で小数を表す方法です。FPU（浮動小数点の計算を行う専用回路）のないマイコンでは、RC から導いた 1 次 IIR の α を 2 のべき乗にしてシフト（ビットをずらす演算。右へ 1 つで 1/2）で計算する書き方がよく使われます。ただし、整数のまま書くと止まります。

// だめな例：差が 64 未満になると増分が 0 になる
y += (x - y) >> 6;

差が 64 未満になると、右シフトの結果が 0 になって値が動かなくなります。入力が 1000 なら、差が 63 になる 937 で止まる計算です。入力を 0 から 1000 へ上げ、途中で 0 へ戻して、3 通りの計算を比べたのが次の図です。

対策は、内部の状態に小数ビットを持たせることです。次のコードは下位 16 ビットを小数部に使う形式で、これを Q16 と呼びます。値を 65536 倍した整数として持つのと同じです（フィルタの Q とは無関係です）。

#include <stdint.h>

typedef struct {
    int32_t state;  // 下位 16 ビットが小数部（Q16）
} iir1_q16_t;

// shift = 6 で alpha = 1/64。x は ±32767 の範囲を想定
static inline int32_t iir1_q16_step(iir1_q16_t *f, int32_t x, int shift)
{
    f->state += (x * 65536 - f->state) >> shift;
    return f->state >> 16;
}

fs = 1 kHz で shift = 6 なら、α = 1/64 で −3 dB は約 2.5 Hz です。設計の注意点は次の 3 つです。

• 桁あふれ：値が変数の範囲を超えること。12 ビット ADC の値（〜4095）なら、状態は 2^28 程度に収まり int32 で足りる。16 ビットの値を入れるときは範囲を計算し直す
• α が小さいときの桁落ち：桁落ちは、小さい値が丸めで消えること。小数ビットが足りないと同じ止まり方が起きる。必要な小数ビットは、少なくとも shift より大きくとる
• 負の数の右シフト：C 言語では処理系定義（結果をコンパイラごとに決める）。GCC は算術シフト（符号ビットを保って右へずらす）だが、移植するコードでは処理系の仕様を見ておく

ここまでは、係数が手元にある前提でした。次の小節で、係数の作り方を見ます。

Python で設計して C へ持っていく

SciPy で係数を設計し、周波数特性を見てから C へ貼り付けるのが手早い流れです。butter() の output="sos" は、biquad を縦につないだ形（second-order sections）で係数を返します。SciPy のドキュメントも、一般のフィルタ処理には sos を勧めています。

import numpy as np
from scipy import signal

fs = 1000.0

# 2 次バターワース LPF（4 次なら biquad 2 段になる）
sos = signal.butter(2, 20, btype="low", fs=fs, output="sos")

# 50 Hz ノッチ
b, a = signal.iirnotch(50, 10, fs=fs)

# 周波数特性を確かめる
f, h = signal.sosfreqz(sos, worN=4096, fs=fs)
i = np.argmin(np.abs(f - 100))
print(f"100 Hz: {20 * np.log10(abs(h[i])):.1f} dB")

# 数列にかけてみる
x = np.random.default_rng(0).normal(size=1 << 16)
y = signal.sosfilt(sos, x)
print(f"variance ratio: {np.var(y) / np.var(x):.4f}")

# C の初期化子として出力する（sos の 1 行は b0 b1 b2 a0 a1 a2）
for s in sos:
    print("{ %.10ff, %.10ff, %.10ff, %.10ff, %.10ff, 0.0f, 0.0f }," % (
        s[0], s[1], s[2], s[4], s[5]))

sosfilt は各段を DF2T で計算し、高い次数での数値誤差を抑える設計です。iirnotch は #4 で扱うノッチの設計関数です。先ほどの C と同じ計算順を Python で単精度（32 ビットの float）に写し、lfilter の倍精度の結果と比べると、2 万サンプルの最大差は 2×10^−6 でした。

作った係数は、自分で書いたコードだけでなく、既存のライブラリにも渡せます。

CMSIS-DSP と ESP-DSP

biquad を自分で書かず、マイコン向けの信号処理ライブラリを使う方法もあります。CMSIS-DSP は Arm の Cortex-M 向け、ESP-DSP は Espressif の ESP32 向けのライブラリです。どちらのライブラリも biquad を部品として持っています。伝達関数は同じでも、内部の形と係数の符号が同じとは限りません。

CMSIS-DSP の関数は arm_biquad_cascade_df2T_f32() です。差分方程式は a_1 y[n] を足す形で書かれていて、ドキュメントにも「設計ツールによっては符号が逆なので、そのときは反転する」と注意があります。SciPy の値は -a1, -a2 にして渡します。係数の並びは 1 段につき {b0, b1, b2, a1, a2} で、段数ぶん続けて並べます。

ESP-DSP の関数は dsps_biquad_f32() です。ソースでは、a0 = 1 とした {b0, b1, b2, a1, a2} を渡し、d0 = x - a1*w0 - a2*w1 の形で計算します。SciPy の値をそのまま使えます。dsps_biquad_gen_lpf_f32() のような係数生成の関数もあり、こちらの周波数は正規化周波数、つまりサンプリング周波数で割った 0〜0.5 の値で渡します。

ここまでは biquad の話でした。波形の形を保ちたいときは、次の小節の FIR に切り替えます。

直線位相が要るなら FIR

波形の形そのものが大事なら FIR を使います。SciPy の firwin は窓関数法で係数を作ります。窓関数法は、理想のローパスのインパルス応答（無限に続く sinc 関数）を有限の長さで切り出し、窓を掛けて係数にする設計法です。

taps = signal.firwin(101, 20, fs=1000.0)   # 101 タップ、カットオフ 20 Hz
y = signal.lfilter(taps, 1.0, x)           # 遅れは (101 - 1) / 2 = 50 サンプル

ここでの 20 Hz は −3 dB ではなく、振幅が半分になる点（−6 dB）です。この例の −3 dB は 15.9 Hz、100 Hz の減衰は −71 dB で、#2 の 2 次バターワースの −28.5 dB より格段に深く落ちます。代わりに遅れは 50 ms で、計算量は 1 サンプルごとに 101 回の掛け算です。

✅ まとめ

• サンプリングで情報が持てるのは fs/2 まで。それより上の成分は低い周波数に化け、デジタルでは取り返せない
• デジタルフィルタは差分方程式で、z 変換すると伝達関数 H(z) になる。z を単位円に置くと周波数特性が出る
• 振幅は dB、位相の傾きは群遅延として読む。極は山、零点は谷を作り、極が単位円の内側にあることが安定の条件
• 白色ノイズの出力は |H|^2 の形の PSD になり、分散比はインパルス応答の 2 乗和で決まる
• PSD はピリオドグラムのままではばらつくので、Welch 法で平均する
• biquad の 3 行のコードは H(z) そのもので、係数を入れ替えれば同じコードで全部のフィルタが動く

ここまでで、コードの式から周波数ごとの倍率と遅れを計算し、ノイズで効き目を測り、係数をマイコンに入れる道具がそろいました。次の一手は、この道具で実際にフィルタを設計することです。#2 ローパスフィルタの設計 では、まず「移動平均 22 点」と「1 次 IIR」と「2 次バターワース」を同じカットオフで比べます。手元の ADC の値に 2 点平均を入れ、分散が半分になることを確かめるところから始めるのもおすすめです。

よくある質問（FAQ）

Q. サンプリング周波数はカットオフの何倍必要ですか？

A. 理論の下限は 2 倍（fc < fs/2）ですが、双一次変換の歪みとアンチエイリアスの余裕を考えると、fc の 10 倍以上あると扱いやすくなります。この連載の fc = 20 Hz、fs = 1 kHz は 50 倍です。逆に 1000 倍を超えるほど差が開くと、係数の丸めが効いてくるので fs を下げることを考えます。

Q. FIR と IIR はどちらが良いのですか？

A. 優劣ではなく目的で選びます。遅れとメモリを節約したいなら IIR、波形の形を保ちたいなら FIR です。同じ 20 Hz でも、2 次バターワースの遅れは約 11 ms、101 タップの FIR は 50 ms です。そのかわり FIR の 100 Hz の減衰は −71 dB まで届きます。

Q. 2 次より高い次数が必要なときはどうしますか？

A. biquad を直列につなぎます。SciPy の butter(N, ..., output="sos") は N 次のフィルタを biquad の列として返します。高い次数を 1 本の分子・分母の多項式で計算すると数値誤差が大きくなるので、段に分けるのが定石です。次数を上げると遅れも増えるので、ループの中では位相余裕も計算し直します。

Q. ADC のあとのデジタルフィルタで、エイリアシングを消せませんか？

A. 消せません。折り返した成分は ADC の時点で本物の信号と同じ数列になるので、見分ける手がかりが残っていません。ADC の前にアナログのアンチエイリアスフィルタを入れるか、fs を上げて折り返す成分を減らします。

Q. FPU のないマイコンでも IIR は使えますか？

A. 使えます。1 次 IIR なら、α を 2 のべき乗にしてシフトで計算できます。ただし整数のままだと差が小さくなったところで止まるので、内部の状態に小数ビットを持たせます。この記事の Q16 のコードがその形です。

関連記事

この連載：

• #2 ローパスフィルタの設計 — 移動平均、1 次 IIR、2 次バターワースを同じカットオフで比べる
• #3 ハイパスフィルタの設計 — DC とドリフトを消す
• #4 バンドパスとノッチの設計 — Q と Audio EQ Cookbook で狙った周波数を扱う

フィルタの計算：

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

センサの入口：

• 加速度センサ KXR94-2050 で作るデジタル水平器 — LPF を入れたくなる代表例。ADC 値から角度を計算する
• Arduino の AD 変換とシリアル通信 — ADC の値を読む基本
• I2C・SPI・UART の違い — デジタル出力のセンサをつなぐ入口
• RS-485・CAN・1-Wire の違い — 基板の外のセンサをつなぐ入口

参考

取得日はいずれも 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「Mathematics of the Discrete Fourier Transform (DFT)」Sampling Theorem
• Alan V. Oppenheim, Ronald W. Schafer「Discrete-Time Signal Processing」3rd ed.（Pearson）

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

• scipy.signal.butter
• scipy.signal.iirnotch
• scipy.signal.freqz
• scipy.signal.group_delay
• scipy.signal.welch
• scipy.signal.sosfilt
• scipy.signal.firwin

マイコン向けライブラリ

• Arm CMSIS-DSP「Biquad Cascade IIR Filters Using a Direct Form II Transposed Structure」
• Espressif ESP-DSP「API Reference」
• Espressif esp-dsp「dsps_biquad_f32_ansi.c」（GitHub）