デジタルフィルタ入門 #11｜レート変換のフィルタの設計：CIC・コム・オールパスで間引きと補間、ΔΣ ADC の sinc フィルタと分数遅延

はじめに

#1 から #10 までは、サンプリング周波数 fs を 1 つに決めたうえで、伝達関数 H(z) を設計してきました。ところが手元の信号処理では、fs そのものを変えたい場面がよくあります。

• ΔΣ 型の ADC の中では、変調器が 1 bit の値を MHz の速さで吐き出しています。欲しいのは 24 bit の値を 1 秒に 10 回や 1000 回です。これを「平均すればいい」で済ませると、どれだけの雑音と妨害が残るのか
• 48 kHz の音声を、DAC の都合で 4 倍の 192 kHz にして出したい。0 を挟んで速さだけ合わせると何が聞こえるのか
• 2 つのセンサのサンプル時刻が 0.3 サンプルだけずれている。整数でないずれを、どうやって合わせるのか

どれも、サンプルの数を減らす（間引き）、増やす（補間）、または整数でない量だけずらす処理です。この記事では、これらをまとめてレート変換と呼びます。そして、レート変換には変える前か後に必ずフィルタが要ります。次の図が、その理由を 4 枚で示しています。

(a) は、振幅 1 の 50 Hz と振幅 0.5 の 900 Hz を足した 8 kHz の列を、そのまま 8 個に 1 個だけ残して 1 kHz にしたものです。1 kHz では 500 Hz より上の周波数を表せないので、900 Hz は 100 Hz に化けて振幅 0.5（−6.0 dB）のまま残りました。一度化けると、元からあった 100 Hz と区別できません。(b) は、先に 400 Hz のローパス（255 タップの FIR）を通してから残したもので、100 Hz の成分は −71.5 dB まで下がりました。(c) は逆向きで、1 kHz の列の間に 0 を 3 個ずつ挟んで 4 kHz にしたものです。元の 50 Hz と 200 Hz の振幅は 4 分の 1 に下がり、950 Hz や 1050 Hz に鏡に映したような成分（イメージ）が −12.0 dB で現れました。(d) は、そこへゲイン 4 のローパスを通したもので、イメージは −69.9 dB に下がり、50 Hz の振幅は 1.00 に戻っています。

この記事のゴールは、手元のレート変換について、比が整数か有理数か、どれだけ大きいか、掛け算をどれだけ使えるかを見分けて方式を選び、段数や語幅や係数を式と数値から決め、折り返しの量と通過域の平らさと遅れで設計を判定できるようになることです。間引きと補間のスペクトルを式で導き、ノーブル恒等式とポリフェーズで計算を減らす理屈を押さえてから、CIC、コム、オールパスの順に 1 段ずつ導き、手計算で確かめます。そのうえで、ΔΣ ADC の 1 bit の列の間引き、ロードセル用の ADC の 50/60 Hz 除去、分数遅延とハーフバンドを設計します。実機は使わず、すべて計算と数値シミュレーションで示します。

この記事で分かること： 間引きでなぜ折り返しが起き、補間でなぜイメージが出るのか（スペクトルの式）、R 個おきの係数のフィルタをなぜ間引きの後ろへ移せるのか（ノーブル恒等式）、ポリフェーズで計算がなぜ減るのか、CIC が掛け算なしで何を計算しているのか、そのゼロ点がなぜ折り返す帯域の真ん中に来るのか、積分器があふれても出力がなぜ正しいのか（ビット成長と巻き戻し）、10 SPS の sinc フィルタがなぜ 50 Hz と 60 Hz を同時に消すのか、振幅を変えないオールパスで何ができるのか（分数遅延、ハーフバンド、90° の位相差）が分かります。

読み終えたときにできること： 間引き比と通過域の幅から CIC の段数と語幅を決め、ドループ（通過域の端での振幅の落ち込み）を補償する短い FIR を作り、sinc フィルタの段数と出力レートから 50/60 Hz の減衰と整定時間を見積もり、Thiran の係数で整数でない遅れを作り、2 倍の間引きをオールパス 2 本のハーフバンドで済ませられます。Python と C で動かし、直接の計算との一致、折り返しの量、群遅延で設計を判定できます。

🔁 レート変換のフィルタとは

レート変換という言葉は、ふだんは「サンプリングレートを変える」くらいの意味で使われます。この記事では、それに伴うフィルタまで含めて、次のように定義します。

順序が定義の一部です。間引きでは、捨てる前にフィルタを通さないと折り返しが起きます。補間では、0 を挟んだ後でフィルタを通さないとイメージが残ります。どちらも、フィルタの役目はローパスで、そのカットオフは低いほうのサンプリング周波数の半分（ナイキスト周波数）で決まります。

語をほどきます。

• 間引き（デシメーション）：サンプルを R 個に 1 個だけ残して、サンプリング周波数を 1/R にする処理。ブロック図では ↓R と書く
• 補間（インターポレーション）：サンプルの間に値を足して、サンプリング周波数を L 倍にする処理。ここでは 0 を挟む操作を ↑L と書き、その後のローパスまで含めて補間と呼ぶ
• 折り返し（エイリアシング）：新しいナイキスト周波数より上の成分が、低い周波数に化けて重なる現象。#1 のサンプリングで起きるものと同じ
• イメージ：0 を挟んだときに、元のスペクトルが新しい帯域の中に鏡に映したように繰り返し現れる成分
• ナイキスト周波数：サンプリング周波数の半分。これより上の周波数は、そのサンプリング周波数では表せない

この記事では、レート変換に使う 3 つの道具を扱います。それぞれ 1 文で定義します。

• CIC フィルタ：積分器とコムだけを直列にした、掛け算を 1 つも使わない間引き・補間用のローパス（Cascaded Integrator–Comb、Hogenauer 1981）
• コムフィルタ： 1 - z^{-M} のように、遅らせた自分を引くフィルタ。周波数軸に等間隔の零点が並ぶ。移動平均の分子であり、ΔΣ ADC の sinc フィルタの中身
• オールパスフィルタ：すべての周波数で振幅が 1 のまま、位相だけを動かす IIR。整数でない遅れ（分数遅延、Thiran）と、2 本の組で作るハーフバンドに使う

何を渡すと何が返るのか

これまでのフィルタとの差を並べます。

時不変とは、入力を 1 サンプル遅らせると、出力も同じ形のまま 1 サンプル遅れる性質です。間引きはこれを満たしません。入力を 1 サンプルずらすと、残るサンプルの組が変わるからです。この性質がどこで効いてくるかは、制御工学の道具で見る節で確かめます。

名前の由来

CIC（Cascaded Integrator–Comb）は、E. B. Hogenauer の 1981 年の論文「An economical class of digital filters for decimation and interpolation」（IEEE Transactions on Acoustics, Speech, and Signal Processing）に由来します。要旨には、高いレートで動く積分器と、同じ数の低いレートで動くコムを直列にし、掛け算器を使わず記憶も少なく済む、と書かれています。積分器とコムの組の数は、折り返しやイメージの誤差の要求から決める、とも書かれています。

コムは櫛（comb）の意味で、振幅特性に等間隔の零点が並ぶ形が櫛の歯に似ていることから来ています。オールパスは「すべてを通す」の意味で、振幅を変えずに全周波数を通すことを表します。P. A. Regalia、S. K. Mitra、P. P. Vaidyanathan の 1988 年の解説「The digital all-pass filter: A versatile signal processing building block」は、オールパスの応用としてノッチ、相補フィルタ、多レートのフィルタ、群遅延の等化、ヒルベルト変換を挙げています。

ポリフェーズ（多相）の名は、M. G. Bellanger らの 1976 年の論文「Digital filtering by polyphase network: Application to sample-rate alteration and filter banks」で、サンプリングレートの変換に使う「ポリフェーズ回路網」として定義されています。ノーブル恒等式とポリフェーズ表現は、Vaidyanathan の 1990 年の解説論文（Proceedings of the IEEE）が多レート（複数のサンプリングレートが混在する）信号処理の基本の部品として整理しています。間引きと補間の理論全体は、R. E. Crochiere と L. R. Rabiner の 1981 年の解説と、1983 年の教科書『Multirate Digital Signal Processing』がまとめています。

前提

この記事の方法がうまく働くのは、次の条件がそろうときです。

1. 比が整数か有理数：R、L、L/R が整数で書ける。任意の比や、時間とともに変わる比は扱わない（Farrow 構造の名前だけ落とし穴で触れる）
2. 欲しい帯域が新しいナイキスト周波数より十分低い：間引きの後の fs,out/2 の手前に、フィルタの遷移帯（通過域から阻止域へ移る区間）を置ける
3. 固定小数点か浮動小数点かが決まっている：CIC は整数のレジスタで作るのが前提で、語幅（レジスタのビット数）を式で決める
4. 分数遅延とハーフバンドは帯域を限る：fs/2 の近くまで完全には作れない。どこまで使うかを帯域で決める

1 から 4 が破れたときに何が起きるかは、それぞれの節と設計例で「破る例」として見ます。

📐 間引きと補間の理論

フィルタを設計する前に、間引きと補間でスペクトルに何が起きるかを式で押さえます。この節の式が、後の CIC やハーフバンドで「どこに零点が要るか」「どこまで平らなら足りるか」を決める物差しになります。

サンプリング定理の復習

#1 のサンプリング定理は、帯域が fs/2 より下に収まった連続時間の信号なら、サンプル列から元の信号を完全に復元できる、というものでした。逆に、fs/2 より上の成分は、サンプリングで低い周波数に折り返します。間引きは、すでにサンプルになった列をもう一度サンプリングし直す処理なので、同じことが離散時間の中で起きます。

間引きのスペクトル：R 枚が重なる

R 個に 1 個を残す操作を式で書きます。

この y のスペクトルは、x のスペクトルを使って次のように書けます（導出は Crochiere–Rabiner、Vaidyanathan の解説にあります）。

式は 2 つのことを言っています。1 つ目は、周波数軸が R 倍に伸びることです。入力レートで θ/R だった成分が、出力レートでは θ に来ます。2 つ目は、それを 2π/R ずつずらした R 枚が重なることです。元の帯域が π/R より広いと、ずらした隣の 1 枚と重なり、重なった部分が折り返しです。重なりを防ぐには、間引く前に π/R（出力レートのナイキスト周波数）より上をローパスで消しておくしかありません。

補間のスペクトル：L 個のイメージ

L 倍に補間するときは、まずサンプルの間に L−1 個の 0 を挟みます。

この u のスペクトルは、次の 1 行で書けます。

周波数軸が 1/L に縮み、0 から 2π の間に元のスペクトルが L 回現れます。そのうち、0 を中心とする 1 つだけが元の信号で、残りの L−1 個がイメージです。カットオフ π/L のローパスでイメージを消すと、元の信号だけが残ります。

ここで振幅が 1/L になる点に注意します。0 を挟むと、サンプルの平均の値が 1/L に下がるからです。はじめにの図 (c) で 50 Hz の振幅が 0.25 になっていたのはこのためです。補間のローパスのゲインを L にすると、(d) のとおり 1 に戻ります。

次の図は、帯域 ±0.6π の三角形のスペクトルを例に、2 つの式を描いたものです。

(b) の和の高さが 0.5 になっているのが、式の 1/R です。帯域 0.6π の信号を R = 2 で間引くと、伸びた帯域は 1.2π になり π を超えるので、0.8π から π の区間で隣の 1 枚と重なりました。(c) の点線の内側（±π/3）だけが残すべき帯域で、外側の 2 つの三角形をローパスで消すのが補間のフィルタの仕事です。

有理数比 L/R：ローパスは 1 本で済む

48 kHz を 44.1 kHz にするような変換は、比が 147/160 の有理数になります。有理数比 L/R の変換は、補間 L → ローパス → 間引き R の順に並べます。

補間のイメージ消しと間引きの折り返し止めは、どちらも同じ高いレートのローパスなので、カットオフの低いほうに合わせた 1 本で両方を兼ねます。順序を逆（先に間引き）にすると、間引きの時点で情報が失われるので、必ず補間を先にします。SciPy の resample_poly は、この 1 本のローパスをポリフェーズで計算する関数です。

ノーブル恒等式：フィルタを低いレートへ移す

間引きの前のローパスは高いレートで回るので、計算が重くなります。そこで、フィルタを間引きの後ろへ移せないかを考えます。移せる条件を述べたのがノーブル恒等式です。

1 行目は、「R サンプルおきの値しか使わないフィルタ」は、間引いた後に低いレートで回しても同じ出力になる、という意味です。H(z^R) の出力のうち、残るサンプルが参照するのは、入力のうち R の倍数の位置だけだからです。2 行目は補間の側で、0 を挟む前に低いレートで回せる、という意味です。CIC のコムが間引きの後ろに置けるのは、この恒等式のおかげです。

ポリフェーズ分解

ふつうの FIR の H(z) は、R サンプルおきの係数だけではありません。そこで、係数を R 組に分けて、それぞれをノーブル恒等式に乗せます。

ポリフェーズ分解は、FIR の係数を「R で割った余り」で R 組に振り分けることです。各組は z^R の関数なので、ノーブル恒等式で間引きの後ろへ移せます。結果として、入力を R 本の列（相）に振り分け、それぞれを短い FIR に通して足す、という形になります。全部の FIR が出力レートで回るので、捨てる出力を計算しません。

次の図は、ノーブル恒等式と、長さ 8 の FIR を R = 2 でポリフェーズにした構造です。

(b) の 2 本の枝は、それぞれ 4 タップの FIR です。1 出力あたりの積和は 4 + 4 = 8 回で、入力 1 サンプルあたりに直すと 4 回です。素直に「長さ 8 の FIR を入力レートで回してから 2 つに 1 つ捨てる」と、入力 1 サンプルあたり 8 回（1 出力あたり 16 回）なので、半分になります。間引きの FIR なら「残す出力だけを計算する」工夫でも同じ数になりますが、ポリフェーズの形にしておくと、補間（0 との掛け算を飛ばす）や、後で見るオールパスのハーフバンドにそのまま使えます。

手計算：長さ 8 の FIR を R = 2 で 1 周

ポリフェーズが「全部計算して 2 つに 1 つ残す」と同じ値になるかを、手で確かめます。係数と入力は次のとおりです。

n = 8 の出力を、まず素直に計算します。 y[8] = \sum_{i=0}^{7} h[i]\, x[8-i] です。

同じ値をポリフェーズで計算します。E_0 は偶数番目の係数 (1, 3, 4, 2) で、偶数番目の入力 x[8], x[6], x[4], x[2] = (5, 2, 5, 4) に掛かります。E_1 は奇数番目の係数 (2, 4, 3, 1) で、1 つ遅らせた奇数番目の入力 x[7], x[5], x[3], x[1] = (6, 9, 1, 1) に掛かります。

一致しました。それぞれの枝は 4 回の積和で、n = 8 の 1 出力に 8 回です。素直な方法では、n = 7 の出力（捨てる側）にも 8 回かかっていました。全部の出力 (3, 15, 35, 64, 91) を Python で確かめるコードは、実装の節にあります。

🧱 CIC フィルタ

R が 64 や 256 のように大きいと、間引きの前のローパスは、入力レートで見て非常に狭い帯域を通す必要があり、FIR なら数百から数千タップになります。CIC は、この重いローパスを掛け算なしで済ませるための構造です。この節では、移動平均から CIC の伝達関数を導き、ノーブル恒等式で構造を軽くし、周波数応答、ドループ、語幅の順に見ます。

移動平均を積分器とコムに分ける

#2 の移動平均は、直近の P 個の和を P で割ったものでした。和の部分だけを伝達関数で書くと、等比数列の和の公式で 2 つに分かれます。

分母の 1/(1 - z^{-1}) は、入力を足し続ける積分器（y[n] = y[n−1] + x[n]）です。分子の 1 - z^{-P} は、P サンプル前の値を引くコムです。積分器で「これまでの総和」を作り、コムで「P サンプル前までの総和」を引けば、直近 P 個の和が残ります。掛け算は 1 つもありません。

CIC は、P = RM としたこの和を N 回重ねたものです。

つまり CIC は、長さ RM の移動和（移動平均の RM 倍）を N 回かけたローパスです。DC のゲインは (RM)^N で、出力をこの値で割れば移動平均の N 乗になります。インパルス応答は長さ RM の矩形を N 回畳み込んだ形（N = 2 で三角形）で、係数は左右対称なので、群遅延は一定の N(RM-1)/2 サンプル（入力レート）です。

ノーブル恒等式でコムを低いレートへ

このままでは、コムが RM サンプル分の遅延を入力レートで持つので、記憶が大きくなります。そこで、コムの式を書き直します。

右辺は「z を z^R に置き換えた 1 − z^{−M}」です。ノーブル恒等式から、間引きの前の 1 - z^{-RM} は、間引きの後の 1 - z^{-M} と同じです。積分器は z^R の形ではないので、間引きの前に残ります。こうして、N 個の積分器を入力レートで回し、R 個に 1 個を残し、N 個のコム（遅延 M）を出力レートで回す、という Hogenauer の構造が出ます。コムの記憶は 1 段あたり M 語で済みます。

次の図の上半分がその構造で、下半分は周波数応答です。

(a) の積分器 3 段とコム 3 段は、どれも足し算か引き算だけです。(b) の零点は fs,out の整数倍に並び、N を増やすほど零点のまわりの谷が深く広くなります。この零点の位置が、CIC がローパスとして都合のよい理由です。

周波数応答：sinc^N と零点の位置

H(z) に z = e^{j\theta} を入れ、 1 - e^{-jx} = 2j e^{-jx/2}\sin(x/2) を使うと、振幅は次の形になります。

R が大きいとき、分母の \sin(\theta/2) は θ/2 で近似でき、DC で 1 に正規化した振幅は sinc 関数（ \mathrm{sinc}(u) = \sin(\pi u)/(\pi u) ）の N 乗になります。f は周波数（Hz）です。

零点は \theta RM/2 = \pi k 、つまり周波数で f = k f_{s,out}/M （k = 1, 2, …）です。M = 1 なら fs,out の整数倍に零点が来ます。

ここが CIC の要です。間引きの式で見たとおり、出力の 0 Hz 付近へ折り返してくるのは、入力の fs,out、2fs,out、3fs,out … の付近の成分でした。CIC の零点は、ちょうどその真ん中に並んでいます。欲しい帯域が 0 \le f \le f_b なら、折り返して帯域に入るのは k f_{s,out} \pm f_b の帯で、その中でいちばん減衰が浅いのは、最初の零点の手前の f_{s,out} - f_b です。

図 (b) の R = 16、b = 0.1 では、fs,out − 0.1fs,out での減衰は N = 1 で −19.2 dB、N = 2 で −38.4 dB、N = 3 で −57.5 dB、N = 4 で −76.7 dB でした。段数を 1 つ増やすごとに同じ dB だけ深くなるのは、振幅が N 乗だからです。

通過域のドループと、帯域の幅で決まる折り返し

sinc は 0 から離れるにつれてなだらかに下がるので、通過域の中でも振幅が落ちます。この落ち込みをドループと呼びます。もう 1 つ、上の式から、折り返しの減衰は帯域の幅 b で大きく変わります。次の図は R = 256 の CIC で、この 2 つを描いたものです。

(a) と (b) は、同じ「N を増やす」が、通過域の平らさには損、折り返しの減衰には得という逆向きに効くことを示しています。M = 2 にすると零点が fs,out/2 に来るので、(a) の破線のとおり通過域が狭くなります。表にまとめます。

2 つ目の表から、CIC が向いているのは「出力レートに比べて欲しい帯域がずっと狭い」場合だと分かります。b = 0.4（fs,out の 40 % まで使いたい）では、N = 4 でも 24 dB しか折り返しを止められません。実際の設計では、CIC で最終のレートの 2 倍や 4 倍まで落とし、残りを FIR で落とすのが定石です。この形は設計例 1 で作ります。

手計算：R = 4、M = 1、N = 1 で 1 周

CIC が本当に移動和を間引いたものになるかを、手で確かめます。入力は整数 8 個です。

積分器で総和を取り、4 個に 1 個（n = 3 と 7）を残し、コムで 1 つ前の値を引きます。

直近 4 個の和を直接計算すると、n = 3 では 3 + 1 + 4 + 1 = 9、n = 7 では 5 + 9 + 2 + 6 = 22 で、一致しました。積分器が大きな値（31）まで育っても、コムの引き算で小さな値に戻ります。次の小節で見る語幅の話は、この「育ってから引く」性質から来ています。

ビット成長と、積分器があふれても正しい理由

CIC の出力は、DC のゲイン (RM)^N の分だけ入力より大きくなります。入力が B_in ビットの 2 の補数（負の数を、ビットを反転して 1 を足した形で表す、ふつうの符号付き整数の表し方）なら、出力を表すのに必要な語幅は次のとおりです（Hogenauer）。

この増え方をビット成長と呼びます。R = 256、M = 1、N = 3 なら 3 \times 8 = 24 ビット増えます。

ところが、積分器の中身はこれでは収まりません。入力に直流が乗っていれば、積分器は際限なく増え続けるからです。それでも CIC が動くのは、次の理由によります。

• 2 の補数（または符号なし）の足し算と引き算は、B ビットのレジスタの中で 2^B を法とする計算（あふれた分を捨てて巻き戻る計算）になる
• 積分器もコムも足し算と引き算だけなので、CIC 全体の計算は「本当の値を 2^B で割った余り」を正しく求めている
• 本当の出力が B ビットの範囲に収まっていれば、余りは本当の値そのものである

つまり、途中の積分器が何周回っても、最後の出力が B_out ビットに収まれば正しい値が出ます。条件は 2 つで、全段を B_out ビット以上の同じ語幅で作ること、あふれたときに飽和（最大値で止める処理）させず巻き戻すことです。Hogenauer の論文は、後段で下位ビットを捨てて語幅を減らす方法（レジスタの刈り込み）も示していますが、この記事では全段を同じ語幅にします。

次の図は、1 bit の ΔΣ 変調器の出力（0 と 1 の列、B_in = 1）を R = 256、M = 1、N = 3 の CIC に通したものです。式の B_out は 1 + 24 = 25 ビットです。

(a) の積分器は、入力から 0.7 ms ほどで 25 ビットを使い切り、その後はのこぎり波のように何度も巻き戻っています。それでも (b) の 25 ビットの出力は、64 ビットで計算した基準と 400 個すべてで一致しました（差 0）。23 ビットにすると、出力の本当の値が 2^{23} を超えるところで余りが本当の値と食い違い、400 個中 199 個が壊れました。

入力が 0/1 の符号なしなので、ここではレジスタも符号なしで読んでいます。出力の値は 0 から 2^{24} の整数で、25 ビットの符号なしなら表せます。入力を ±1 の 2 の補数（2 ビット）で入れるなら、B_in = 2 で 26 ビットが要ります。この正弦波の振幅 0.5 の例では出力が 0.25 \times 2^{24} から 0.75 \times 2^{24} の間にあるので 24 ビットでも壊れませんが、入力がすべて 1 になると出力は 2^{24} になり、24 ビットでは 0 に化けます。語幅は、最大の入力のときの出力で決めます。

ドループを打ち消す補償 FIR

CIC の後ろには、ドループを打ち消す短い FIR を置くのが定石です。欲しい特性は、通過域で 1/\lvert H_{CIC}\rvert 、つまり sinc の逆数で、周波数とともにゆるやかに持ち上がる形です。

この記事では、出力レートで動く対称な FIR の係数を、0 から f_b までの 1/\lvert H_{CIC}\rvert に最小 2 乗で合わせて決めます（対称なので位相は直線で、群遅延は一定）。R = 256、N = 3、fs,out = 4 kHz、f_b = 1.6 kHz で 15 タップにすると、CIC との積は 0〜1.6 kHz で ±0.0084 dB に収まりました。係数は中央が 1.634、両隣が −0.428 で、外へ行くほど小さく符号が交互に変わります。周波数応答は、制御工学の道具の節のボード線図で見ます。

ただし、補償 FIR は高い周波数を持ち上げるので、1.6 kHz で +7.3 dB、fs,out/2 の 2 kHz では +9.6 dB のゲインを持ちます。この帯域に折り返したノイズがあれば、それも持ち上げます。設計例 1 で、その代償を数で確かめます。

補間用の CIC

補間では、構造を逆に並べます。出力レートの側（低いレート）でコムを N 段回し、0 を R−1 個挟み、高いレートで積分器を N 段回します。コムを先に置くのは、ノーブル恒等式の 2 行目で、0 を挟む前に移したからです。

0 を挟むと平均の値が 1/R になるので、DC のゲインは次のとおりです。

R = 4、M = 1、N = 3 なら 64/4 = 16 で、入力に 1 を並べると出力は 16 に落ち着きました。補間用の CIC は、積分器が最後に来るので、出力の直流が積分器の中に残り続けます。入力に直流があっても発散しないのは、前段のコムが直流をいったん 0 にしているからです。

〰️ コムフィルタと ΔΣ ADC の sinc フィルタ

CIC の中にあったコムは、単独でも使えるフィルタです。この節では、コムの零点を見たうえで、ΔΣ ADC の中の sinc フィルタが CIC そのものであること、そして出力レートを選ぶだけで電源の 50 Hz と 60 Hz を同時に消せることを導きます。この節の後半は、ロードセルや熱電対を ADC で読むときに、データシートの数字を読み解くのにそのまま使えます。

コムの零点

いちばん簡単なコムは、M サンプル前の自分を引く FIR です。

振幅は \lvert H_c(e^{j\theta})\rvert = 2\lvert\sin(\theta M/2)\rvert で、零点は \theta = 2\pi k/M 、周波数では f = k f_s/M （k = 0, 1, 2, …）です。

k = 0 も零点なので、コムだけでは直流も消えます。ここに積分器 1/(1 - z^{-1}) を掛けると、DC の零点が極と打ち消し合い、M で割れば長さ M の移動平均になります。

#2 の移動平均は、コムと積分器の組だったと読み替えられます。#4 のノッチが 1 つの周波数だけを消したのに対し、コムは等間隔の多数のノッチを一度に置きます。

分母に遅れを置くと、逆の形の IIR コムになります。

こちらは同じ間隔で鋭いピークが並びます。音響では、短い遅れの帰還を重ねて残響（リバーブ）の響きを作るのに使われてきました。この記事では名前だけにとどめます。

次の図の (a) が、コム・移動平均・IIR コムの振幅です。残りの 3 枚は、すぐ後の sinc フィルタの話で使います。

(a) の移動平均 8（橙）が、コム（青）の DC の零点だけを取り除いた形になっていることが見て取れます。IIR コム（紫）は、コムの零点の位置でなく、その間の山の位置にピークを立てています。

ΔΣ 変調と sinc フィルタ

ΔΣ（デルタシグマ）変調は、入力と出力の差を積分して 1 ビット（または少ないビット）で量子化し、それを帰還する回路で、入力より何十倍も速く動かすことで、量子化の誤差（量子化ノイズ）を高い周波数へ押しやる方式です。押しやる働きをノイズシェーピングと呼びます。変調器の出力を低い周波数だけ残してから間引けば、少ないビットの高速な列から、多いビットの低速な値が得られます。ΔΣ 型の ADC は、この変調器と、間引きのデジタルフィルタを 1 チップに入れたものです。この記事では変調器の回路には踏み込まず、出てくる 1 bit の列をどう間引くかだけを扱います。

この間引きのフィルタに、CIC がそのまま使われています。データシートでは「sinc3 フィルタ」「sinc4 フィルタ」と呼ばれ、数字は段数 N を表します。周波数応答が sinc の N 乗だからです。たとえば Texas Instruments のブリッジセンサ用 24 ビット ADC である ADS1232 のデータシートは、sinc4 のデジタルフィルタを使い、変調器のサンプリング周波数は 76.8 kHz（クロック 4.9152 MHz のとき）、出力は 10 SPS か 80 SPS（SPS：1 秒あたりのサンプル数）と記しています。10 SPS なら R = 76800/10 = 7680 です。

10 SPS は 50 Hz と 60 Hz を同時に消す

出力レートが fs,out の sinc^N は、fs,out の整数倍に零点を持ちます。fs,out = 10 SPS なら、零点は 10、20、30、40、50、60 Hz、…です。50 Hz にも 60 Hz にも零点があるので、商用電源から乗る誘導ノイズを、東日本（50 Hz）でも西日本（60 Hz）でも同時に消せます。ADS1232 のデータシートも、10 Hz の整数倍が除去されるので 50 Hz と 60 Hz を同時に除去できる、と書いています。

ただし、電源の周波数は少し揺れ、ADC のクロックにも誤差があります。零点からずれたときにどれだけ消えるかが、段数 N で変わります。図の (c) が、その近くの拡大です。

表の値は、それぞれの区間の端で、いちばん浅い値です。N = 1（単純な平均）では、零点から 2 % ずれただけで 34 dB しか消えません。N = 3 や 4 では、零点の谷が N 倍の dB で深くなるので、2 % ずれても 100 dB 以上が残ります。データシートが sinc4 の利点として「零点ごとに 4 つの零が同じ場所にある」と書き、内部発振器の ±3 % のばらつきでも 100 dB 以上の除去が得られる、としているのはこのためです（ADS1232 の電気的特性の表では、内部発振器・10 SPS・50 Hz または 60 Hz ±1 Hz で最小 100 dB、標準 110 dB）。

比べるために、出力レートを 50 SPS にした場合を計算すると、零点は 50 Hz の整数倍で、60 Hz は零点に乗りません。60 Hz での減衰は N = 1 で −16.1 dB、N = 3 で −48.4 dB、N = 4 で −64.6 dB にとどまります。50 と 60 の両方を消したいなら、出力レートは 10 SPS（またはその整数分の 1）を選びます。

段数の代償：遅れと整定時間

段数を増やすと、通過域も狭くなり、応答も遅くなります。理想化した sinc^N（R が大きい極限）で −3 dB の周波数を計算すると、10 SPS で N = 1 は 4.42 Hz、N = 3 は 2.62 Hz、N = 4 は 2.27 Hz でした。ADS1232 のデータシートは 10 SPS の −3 dB を 2.4 Hz としていて、この記事の理想化した計算の 2.27 Hz とは少しずれます。

整定時間（入力が変わってから出力が新しい値に落ち着くまでの時間）は、インパルス応答の長さで決まります。sinc^N のインパルス応答は、出力レートで数えて N サンプル分の長さです。図の (d) は、出力周期の途中（2.37 周期の時刻）にステップが入ったときの出力です。横軸は出力の番号です。N = 1 は番号 2 で 0.63、番号 3 で 1 に届きました。N = 3 は番号 2 で 0.04、番号 3 で 0.60、番号 4 で 0.99、番号 5 で 1 です。N = 4 は番号 6 で 1 に届きます。ステップが出力周期の途中で入ると、ステップを含む最初の出力（ここでは番号 2）は一部だけを含むので、それを 1 個目と数えると、完全に新しい値になるのは N + 1 個目です。ADS1232 のデータシートも、入力が急に変わったときは 4 回分のデータ周期で整定し、ステップの後に完全に整定したデータを得るには 5 回の読み出しが要る（最初の 4 回は捨てる）と書いています。10 SPS の sinc4 なら、それは 0.5 秒です。

🌀 オールパスフィルタ

ここまでのフィルタは、振幅を削ることで仕事をしてきました。オールパスは逆に、振幅には一切触れず、位相（と群遅延）だけを動かします。それだけで、整数でない遅れ、2 倍の間引き、90° の位相差という 3 つの仕事ができます。この節では 1 次のオールパスの性質を導いてから、その 3 つを順に作ります。

1 次のオールパス：振幅がちょうど 1 になる理由

1 次のオールパスは、係数 a を 1 つだけ持つ IIR です。

差分方程式は y[n] = a\,x[n] + x[n-1] - a\,y[n-1] で、 y[n] = a\,(x[n] - y[n-1]) + x[n-1] と書き直せば掛け算は 1 回です。

分子と分母が「係数の並びを逆にした」鏡像の関係にあることが、振幅 1 の理由です。 z = e^{j\theta} を入れて、分子から e^{-j\theta} をくくり出します。

右辺のかっこの中は、分母 1 + a\,e^{-j\theta} の複素共役です。共役どうしの絶対値は等しく、 \lvert e^{-j\theta}\rvert = 1 なので、すべての θ で \lvert A(e^{j\theta})\rvert = 1 です。数値で確かめると、a = −0.9、−0.5、0.3、0.9 の 4 通りで、1 からのずれは最大 1.6 × 10⁻¹⁵ でした（倍精度の丸めの大きさ）。

位相と群遅延

同じ分解から位相も出ます。分子のかっこの偏角を ψ とおくと、分母の偏角は −ψ です。

#1 の群遅延は \tau(\theta) = -d\phi/d\theta でした。ψ を微分して整理すると、群遅延は次の形になります。

途中を 1 段ずつ書くと、 d\psi/d\theta = (a\cos\theta + a^2)/(1 + 2a\cos\theta + a^2) で、 \tau = 1 - 2\,d\psi/d\theta の分子は 1 + 2a\cos\theta + a^2 - 2a\cos\theta - 2a^2 = 1 - a^2 になります。

安定の条件は極 −a が単位円の内側にあること、つまり \lvert a\rvert \lt 1 です。このとき分子の 1 - a^2 は正なので、群遅延はどの周波数でも正です。代表の値を並べます。

a = 0 は単なる 1 サンプルの遅れです。a を負にすると低域が遅れ、正にすると高域が遅れます。DC の群遅延は (1-a)/(1+a) で、ここから「DC 付近で D サンプル遅らせたい」なら a = (1-D)/(1+D) と決められます。これが 1 次の Thiran で、すぐ後で一般の次数へ広げます。

次の図は、この位相と群遅延と、すぐ後の手計算の結果です。

(b) のとおり、群遅延は周波数で大きく変わります。オールパスは振幅を変えないかわりに、波形の各成分を周波数ごとに違う時間だけ遅らせます。この「周波数で変わる遅れ」を、遅れを作りたい帯域で平らにするのが分数遅延で、逆に積極的に使うのがハーフバンドとヒルベルトの組です。

手計算：a = 0.5 に fs/4 の余弦を通す

振幅が変わらず位相だけが遅れることを、4 サンプルで確かめます。入力は fs/4 の余弦 x[n] = \cos(\pi n/2) 、つまり 1, 0, −1, 0 の繰り返しです。十分時間がたった後の出力（定常状態）を y = (0.8, 0.6, −0.8, −0.6) と見込み、差分方程式 y[n] = 0.5\,x[n] + x[n-1] - 0.5\,y[n-1] に代入して、自分自身と矛盾しないかを見ます。

4 行とも見込みどおりの値に戻り、この列が定常状態です。4 点の DFT の 1 番目のビン（fs/4 の成分）を、入力と出力で比べます。

\lvert Y[1]\rvert = \sqrt{1.6^2 + 1.2^2} = 2 で、振幅は入力と同じです。偏角は \arctan(-1.2/1.6) = -36.87^\circ で、位相の式 -90^\circ + 2\arctan(0.5) = -90^\circ + 53.13^\circ と一致します。3:4:5 の直角三角形がそのまま出てくる、きれいな例です。

この遅れ 36.87° を時間に直すと、fs/4 の 1 周期 4 サンプルの 36.87/360 で 0.41 サンプルです。表の群遅延 0.60 とは違います。0.41 は位相遅延（その周波数の正弦波そのものが遅れる時間 −φ/θ）、0.60 は群遅延（その周波数の付近の包絡線が遅れる時間）で、オールパスでは 2 つが一致しません。同じ計算を Python で確かめるコードは、実装の節にあります。

分数遅延：Thiran の最大平坦群遅延

2 つのセンサのサンプル時刻が 0.3 サンプルずれているとき、片方を 0.3 サンプルだけ遅らせられれば時刻がそろいます。整数でない遅れを作るフィルタを分数遅延フィルタと呼びます。理想の伝達関数は z^{-D} （D は整数でない）で、振幅は 1、群遅延はどこでも D です。振幅がもともと 1 のオールパスなら、残る課題は「群遅延をどれだけの帯域で D に保てるか」だけです。

J.-P. Thiran は 1971 年の論文で、群遅延を DC のまわりで最大限に平ら（DC での群遅延の微分ができるだけ高い次数まで 0）にする再帰形のフィルタを、z 平面で直接導きました。オールパスの形で書くと、次のとおりです（Laakso らの 1996 年の解説の形）。

K = 1 では a_1 = (1-D)/(1+D) となり、前の小節で 1 次オールパスの DC の群遅延から決めた値と一致します。

K = 3、D = 3.3 の係数は、 a_1 = -0.209302 、 a_2 = 0.051338 、 a_3 = -0.006248 です。極の絶対値は 0.206、0.206、0.148 で、どれも単位円の十分内側にあります。

D の選び方には 2 つの条件があります。1 つは安定で、極が単位円の内側に入るのは D \gt K - 1 のときです。もう 1 つは帯域で、D が次数 K に近いほど群遅延が平らな帯域が広く、 K - 0.5 \le D \le K + 0.5 がよく使われる目安です。0.3 サンプルずらしたいなら、整数部を K に合わせて D = K + 0.3 とし、もう一方の信号を K サンプル（整数）遅らせてそろえます。

次の図は、K = 3 で D が範囲の中（3.3）と外（2.3、1.3）のとき、それに FIR のラグランジュ補間（多項式で間の値を推定する FIR。1 次なら線形補間）を並べたものです。

(a) の Thiran D = 3.3 は、0.1 fs で 3.298、0.2 fs で 3.235 と目標に張り付きます。範囲外の D = 2.3 は、0.2 fs までは 2.31 と悪くありませんが、0.3 fs で 2.44、0.4 fs で 3.29 と大きくそれます。(c) の D = 1.3 は、 D \gt K - 1 = 2 を破るので極が 1.98 になり、インパルス応答が増え続けます。

ラグランジュ補間の FIR は、係数が式で決まり、どんな D でも安定という利点があります。そのかわり、表のかっこのとおり振幅が高域で落ちます。Laakso らの解説は、FIR（ラグランジュなど）とオールパス（Thiran など）の分数遅延を、係数の更新のしやすさと計算量も含めて整理しています。振幅を落としたくないならオールパス、D を毎サンプル変えたいなら FIR（その代表が落とし穴の節で名前を挙げる Farrow 構造）という分け方になります。

2 本のオールパスで作るハーフバンド

ハーフバンドフィルタは、通過域と阻止域が fs/4 を中心に対称なローパスで、2 倍の間引きと補間に使います。FIR で作ると、中央以外の係数が 1 つおきに 0 になる性質があります。ここでは、Vaidyanathan や Regalia らの解説で取り上げられている、2 本のオールパスを足し引きして作る IIR のハーフバンドを導きます。

2 本の枝の位相が低域でそろい、高域で 180° ずれるように c_0、c_1 を選べば、足した H_LP は低域を通して高域を打ち消し、引いた H_HP はその逆になります。z² の 1 次が 2 本と z^{−1} で、全体の次数は 5 です。

この形には、係数をどう選んでも成り立つ性質があります。 p = A_0/2 、 q = z^{-1}A_1/2 とおくと、平行四辺形の恒等式 \lvert p+q\rvert^2 + \lvert p-q\rvert^2 = 2\lvert p\rvert^2 + 2\lvert q\rvert^2 と、オールパスの振幅 1 から、次が出ます。

この性質をパワー相補と呼びます。ローパスで削った電力は、ちょうどハイパスに出てきます。fs/4 では z^{-2} = -1 なので、2 本の枝の位相差は c_0、c_1 によらず 90° になり、 \lvert H_{LP}\rvert^2 = 1/2 、つまり −3.01 dB です。

この記事では、通過域 0〜0.2 fs、阻止域 0.3〜0.5 fs として、阻止域の最大を最小にするように c_0、c_1 を数値最適化で決めました。 c_0 = 0.236471 、 c_1 = 0.714542 で、阻止域の減衰は −36.2 dB、通過域のリプル（振幅の波打ち）は 0.001 dB です。

(a) の破線（和）は、数値で 1 からのずれが最大 8.9 × 10⁻¹⁶ でした。(b) の位相差（破線）が低域で 0°、高域で 180° になっているのが、足して通し、引いて通す仕組みそのものです。

ポリフェーズにすると、このフィルタは非常に軽くなります。A_0(z²) と A_1(z²) は z² の関数なので、ノーブル恒等式で 2 倍の間引きの後ろへ移せます。偶数番目の入力を A_0 に、1 つ遅らせた奇数番目の入力を A_1 に通し、出力レートで足して 1/2 にするだけです。1 次オールパスは掛け算 1 回なので、1 出力あたりの掛け算は 2 回です。同じ阻止域の FIR との比較は、設計例 3 で行います。

ヒルベルト変換をオールパスの組で作る

#10 のヒルベルト変換は、すべての周波数の位相を −90° 回し、振幅を変えないフィルタでした。#10 では、それを窓をかけた FIR で近似し、位相はちょうど −90° で、誤差は振幅に出ました。逆に、振幅を 1 に固定して誤差を位相に逃がす作り方があります。2 本のオールパスを並べて、出力どうしの位相差が帯域の中でほぼ 90° になるように係数を選ぶ方法です。

ヒルベルト変換そのものとの違いは、2 本の枝の両方が位相を動かすことです。どちらの出力も元の信号とは位相が違いますが、2 本の間の位相差だけが 90° になります。解析信号（#10 の、実部と虚部が 90° ずれた複素数の列）を作る目的なら、H_I の出力を実部（同相成分、in-phase）、H_Q の出力を虚部（直交成分、quadrature）にすれば足ります。

0.02〜0.48 fs で位相差の 90° からのずれの最大を最小にするように、4 つの係数を数値最適化で決めました。 c_{I} = (-0.132005,\ -0.694810) 、 c_{Q} = (-0.420224,\ -0.904797) で、ずれは最大 0.40° です。#10 の FIR（L = 63、ハミング窓）と並べます。

(a) と (b) は、誤差の置き場所が FIR とオールパス対で正反対であることを示しています。(c) のとおり、オールパス対の遅れは 0.25 fs で約 1.9 サンプルと短い一方、帯域の端では 10 サンプルを超え、周波数で変わります。2 本の枝の群遅延はほぼ重なっているので、2 本の間の時刻はそろっています。

性能を 1 つの物差しで比べるため、解析信号に残る負の周波数の成分（本来 0 になるべき成分）の大きさを計算しました。振幅の誤差も位相の誤差も、ここに漏れとして出ます。

オールパス対の漏れは、0.02〜0.48 fs の中で最大 −49.2 dB でした。掛け算は 4 分の 1、遅れも短く済みます。そのかわり、波形の形（成分ごとの時刻のそろい方）は変わります。包絡線や瞬時周波数だけが欲しい用途ならオールパス対、波形の形を保ちたいなら FIR、という使い分けになります。

群遅延の等化

オールパスのもう 1 つの使い道を、名前だけ挙げておきます。#2 のバターワースのような IIR は、通過域の中で群遅延がうねります。その逆向きのうねりを持つオールパスを直列につなぐと、振幅を変えずに群遅延を平らに近づけられます。これを群遅延の等化と呼び、Regalia らの解説が応用の 1 つに挙げています。遅れの合計は増えるので、制御ループの中ではなく、記録や通信の波形の形をそろえたい場面で使います。

📊 周波数と制御工学の道具で見る

ここまでの道具を、連載で使ってきた周波数と制御工学の物差しで読み直します。レート変換は、制御ループの中ではセンサの前処理として現れることが多いので、位相余裕と、折り返した成分がループにどう入るかを中心に見ます。

CIC をループの中に入れると：群遅延と位相余裕

#2 の「制御ループに入れると何が起きるか」と同じループを使います。モータはゲイン 1、時定数 50 ms の 1 次遅れ、1 ms 周期（1 kHz）の PI 制御（ K_p = 3 、 K_i = 60 ）で、フィルタなしのゲイン交差（開ループのゲインがちょうど 1 になる周波数）は約 9.6 Hz、位相余裕（ゲイン交差での位相が −180° までどれだけ残っているか）は 88.4° でした。

ここで、速度センサを 8 kHz でサンプルし、CIC で 8 分の 1 に間引いてから 1 kHz のループに渡す構成を考えます。CIC の群遅延は、入力レートで数えて N(RM - 1)/2 サンプルです。遅れ τ_d の分だけ、ゲイン交差での位相が 360^\circ \times f_{gc} \times \tau_d 減ります。

まず、#2 と同じ −3 dB 20 Hz にそろえて比べます。R = 8 のまま M を変えて −3 dB を 20 Hz に合わせると、N = 1 は M = 22（長さ 176 の移動平均を 8 kHz で回すのと同じで、#2 の移動平均 22 を 1 kHz で回すのと同じ時間幅）、N = 2 は M = 16、N = 3 は M = 13 になりました。

同じ −3 dB にそろえると、段数が多いほど遅れが長く、位相余裕が減ります。sinc^N は N が大きいほど通過域の肩がなだらかなので、同じ −3 dB にするには長い窓が要るからです。N = 3 では 26.0° まで下がりました。CIC を「ループの帯域を決めるローパス」として使うのは、遅れの面で得になりません。

一方、表の最後の行は、M = 1 の CIC を「折り返し止め」としてだけ使う構成です。−3 dB は 263 Hz と高く、遅れは 1.3 ms、位相余裕は 83.9° で、なしとの差は 4.5° です。ループの帯域を決めるのは後段のフィルタに任せ、CIC は 8 kHz から 1 kHz に落とすときの折り返しを止めるだけ、という役割分担なら、ほとんど位相を失いません。

もう 1 つ、間引きのタイミングにも注意します。CIC の出力のタイミングが、ループのサンプル時刻とそろっていないと、最大で出力の 1 サンプル（ここでは 1 ms）の余分な遅れが入ります。9.6 Hz のゲイン交差では、 360^\circ \times 9.6 \times 0.001 = 3.5^\circ です。表の位相余裕は、出力の時刻が最新の入力とそろっている場合の値です。

折り返した成分は「見かけの低周波」の外乱になる

次に、折り返しがループに何をするかを見ます。センサに 1010 Hz の機械振動（振幅 0.2）が乗っているとします。1 kHz で直接サンプルすると、1010 Hz は #1 の折り返しの式のとおり \lvert 1010 - 1000\rvert = 10 Hz に化けます。10 Hz はループのゲイン交差（9.6 Hz）のすぐそばで、ループはこれを本物の速度の変動だと思って追いかけます。

(b) の赤は、1 kHz で直接サンプルした場合の速度です。センサの値ではなく、モータの速度そのものが 300 ms 以降も幅 0.28（0.86〜1.14）で揺れ続けました。化けた 10 Hz を打ち消そうとして、ループが本当にモータを揺らしたからです。緑は 8 kHz でサンプルして CIC（R = 8、M = 1、N = 3）で落とした場合で、揺れの幅は 1.8 × 10⁻⁵ でした。1010 Hz は CIC の最初の零点（1000 Hz）のすぐそばにあり、N = 1 で −39.9 dB、N = 2 で −79.7 dB、N = 3 で −119.6 dB 消えるからです。

一度折り返した成分は、どんなフィルタでも本物の 10 Hz と区別できません。折り返しは間引く前にしか止められない、というのが、この小節の結論です。

オールパスの群遅延は周波数で変わる

オールパスの節で見たとおり、1 次オールパスの群遅延は (1-a^2)/(1 + 2a\cos\theta + a^2) で、周波数とともに変わります。分数遅延として使えるのは、群遅延が D に張り付いている帯域だけです。Thiran 3 次 D = 3.3 では、群遅延のずれが 0.1 fs で 0.002 サンプル、0.2 fs で 0.065 サンプル、0.3 fs で 0.34 サンプルでした。「0.3 サンプルずらす」という目的に対して、0.3 fs ではずれが目的と同じ大きさになるので、使える帯域はおよそ 0.2 fs までです。設計例 3 では、これを「ずらした後の誤差」の帯域として数で決めます。

ボード線図で見る CIC と補償 FIR

#1 のボード線図（振幅と位相を周波数に対して描いた図）で、CIC のドループと補償 FIR を重ねます。次の図の (b) がそれで、(a) は次の小節で使います。

(b) の黒い線（積）は、0〜1.6 kHz で ±0.0084 dB に収まっています。デシベルでは掛け算が足し算になるので、補償 FIR の曲線は CIC の曲線を 0 dB の線で折り返した形になっています。1.6 kHz より上では補償が追いつかず、積は 2 kHz で −2.2 dB まで下がります。補償 FIR は対称な係数なので位相は直線で、群遅延は 7 サンプル（出力レート）で一定です。

多レート系は時不変でない

最後に、間引きが時不変でないことを確かめます。図の (a) は、係数 (1, 2, 3, 2, 1)/9 の FIR を通してから 2 つに 1 つ残す系に、インパルス δ[n] と、それを 1 サンプル遅らせた δ[n−1] を入れたものです。出力は (0.11, 0.33, 0.11, 0, …) と (0, 0.22, 0.22, 0, …) で、後者は前者をずらした形になっていません。残るサンプルの組が入れ替わり、FIR の偶数番目の係数を見るか奇数番目を見るかが変わったからです。

時不変でない系には、厳密には 1 つの伝達関数が定まりません。入力をずらすたびに応答が変わるからです。それでもこの記事でこれまで伝達関数で議論してきたのは、間引きの前のフィルタが出力レートのナイキスト周波数より上を十分に消していれば、出力に効くのは通過域の成分だけで、それは「低いレートの線形時不変なフィルタ」として扱って差し支えないからです。折り返しが残っているあいだは、この近似が成り立ちません。逆に言えば、折り返しの減衰を設計の物差しにすることは、多レートの系を伝達関数で扱ってよい範囲を決めることでもあります。

🛠️ 設計例

ここからは、要件 → 方式 → パラメータ → 検証の順に、3 つのレート変換を設計します。どれも計算と数値シミュレーションで、値はすべて記事の Python で再現できます。

設計例 1：ΔΣ 変調器の 1 bit の列を間引く

要件： 1.024 MHz で出てくる 1 bit の列（0 と 1）から、4 kHz の多ビットの値を作る。欲しい帯域は 0〜1.6 kHz で、その中の振幅の平らさは ±0.1 dB。

入力の作り方： 実物の変調器の代わりに、2 次の ΔΣ 変調器を数式どおりに回して 1 bit の列を作りました。入力 x[n] は振幅 0.5 の 125 Hz の正弦波（±1 を全振幅とする）に、標準偏差 10⁻⁶ のごく小さなノイズを足したものです。積分器を 2 段持ち、出力 v[n] を ±1 として帰還します。

方式と段数： R = 256 の CIC で 4 kHz に落とします。段数は、量子化ノイズの形から決めます。次の図の (b) のとおり、2 次の変調器のノイズは周波数とともに 1 桁あたり 40 dB で増えます。一方、CIC の零点の付近では、振幅は零点からの距離の N 乗に比例して小さくなります（1 桁あたり 20N dB）。折り返してくるのは零点のまわりの帯なので、N が 2 だとノイズの増え方と減衰が釣り合い、N が 3 で初めて減衰が勝ちます。

(a) では、0 と 1 が細かく入れ替わり、その 1 の割合（密度）が入力の値をなぞっています。(c) の SNR（0〜2 kHz の信号とノイズの電力の比）は、N = 1 で 39.3 dB、N = 2 で 83.2 dB、N = 3 で 102.1 dB でした。単純な平均（N = 1）では、高域にあったノイズが折り返して床を 60 dB 以上押し上げています。N = 3 を選びます。

語幅： 入力は 0/1 の B_in = 1、 N\log_2(RM) = 3 \times 8 = 24 で、B_out = 25 ビットです。CIC の節の図のとおり、25 ビットの符号なしレジスタで 64 ビットの計算と完全に一致しました。

ドループと補償： 1.6 kHz（0.4 fs,out）のドループは −7.26 dB です。CIC の節で作った 15 タップの補償 FIR を通すと、0〜1.6 kHz で ±0.0084 dB に収まり、要件の ±0.1 dB を満たします。ただし、補償 FIR が 2 kHz 付近を +9.6 dB 持ち上げるので、SNR は 102.1 dB から 97.2 dB に下がりました。0〜1.6 kHz だけで数えた SNR は 103.3 dB で、下がった分は帯域の外（1.6〜2 kHz）のノイズです。

破る例：通過域が広いと折り返しが止まらない： ここまでの数字は、入力が 125 Hz の正弦波だけで、ほかの周波数の信号がない場合です。CIC の節の表のとおり、b = 0.4 では N = 3 でも折り返しの減衰は −17.8 dB しかありません。入力に 2.4 kHz の信号（妨害）があると、4 kHz に落とした時点で 1.6 kHz に化け、元の 17.8 dB 下の大きさで帯域の中に残ります。量子化ノイズは 2.4 kHz 付近では小さいので SNR には出てきませんが、入力の信号そのものの折り返しは止められていません。

直し方：2 段に分ける： CIC で 2 倍のレート（8 kHz）まで落とし、そこで補償と折り返し止めを兼ねた FIR で 2 分の 1 にします。

FIR の係数は、0〜1.6 kHz で CIC の逆数に、2.4〜4 kHz で 0 に合わせる重み付きの最小 2 乗を、誤差の大きいところの重みを上げながら繰り返して求めました（振幅の誤差の最大を小さくする工夫）。2 段構成の SNR は 0〜2 kHz で 99.4 dB、0〜1.6 kHz で 102.0 dB で、2.4 kHz の妨害は 89.8 dB 以上消えます。CIC の 1 段目の折り返しは、b = 1.6/8 = 0.2 なので、8 kHz の倍数 ± 1.6 kHz の帯で −37.9 dB です。この帯の信号が問題になるなら、変調器の前のアナログのローパスで落とします。ADS1232 のデータシートも、デジタルフィルタの応答は変調器のサンプリング周波数（76.8 kHz）の倍数で繰り返すと書き、PGA（プログラマブルゲインアンプ。ADC に内蔵された、利得を切り替えられる増幅器）の出力に 0.1 µF の外付けコンデンサを付けると内部の 2 kΩ と合わせて 720 Hz のローパスになり、変調器に入る前に帯域を制限して折り返しを防ぐ、と説明しています。

設計例 2：ロードセル用 ADC の 10 SPS と 50/60 Hz 除去

要件： ロードセル（ひずみゲージのブリッジ）の出力を 24 ビットの ΔΣ ADC で読む。配線に商用電源の 50 Hz または 60 Hz が乗り、電源の周波数とクロックのずれを合わせて ±2 % を見込む。そのずれの範囲で 100 dB 以上消したい。重さが変わってから 1 秒以内に値が落ち着いてほしい。

方式： sinc^N の零点を 50 Hz と 60 Hz の両方に置くため、出力レートは 10 SPS にします。コムの節の図のとおり、零点は 10 Hz の整数倍に並びます。

段数： コムの節の表から、±2 % の端での減衰は N = 3 で −101.8 dB、N = 4 で −135.8 dB です。N = 3 でも要件の 100 dB にわずかに届きますが、±3 % では −91.6 dB まで下がります。余裕を見て N = 4（sinc4）を選びます。ADS1232 はこの構成（sinc4、10 SPS）で、データシートの値は内部発振器で最小 100 dB、外部クロックで最小 120 dB（50 Hz または 60 Hz の ±1 Hz）です。

整定時間： sinc4 では、入力が変わってから 5 個目の出力で完全に整定します（最初の 4 個を捨てる）。10 SPS なら 0.5 秒で、要件の 1 秒に収まります。N = 3 なら 4 個目で 0.4 秒です。

破る例：出力レートの選び方： 出力レートを 50 SPS にすると、60 Hz は零点に乗らず、sinc4 でも 60 Hz の減衰は −64.6 dB です。出力レートが 50 Hz と 60 Hz の公約数（10、5、2、1 Hz と、その整数分の 1）でないと、どちらか、または両方が零点から外れます。たとえば 15 SPS では 60 Hz は零点に乗りますが、50 Hz は乗りません（零点は 45 と 60 Hz）。50/60 Hz の両方を零点に乗せられる出力レートは、10 Hz を整数で割った値（10、5、3.33、2.5 SPS …）です。

マイコンの ADC で同じことをするなら： マイコン内蔵の 12 ビット ADC を 1 kHz で回し、100 個ずつ平均して 10 SPS にするのは、R = 100 の N = 1 の CIC と同じです。±2 % では 34 dB しか消えません。同じ 1 kHz のサンプルを N = 3 の CIC（R = 100、M = 1）に通すと、零点の位置は同じまま、±2 % で約 102 dB になります。語幅は 12 + \lceil 3\log_2 100 \rceil = 12 + 20 = 32 ビットで、uint32_t にちょうど収まります。代償は、整定が 4 個目（0.4 秒）になることです。なお、1 kHz のサンプリングそのものの折り返し（500 Hz より上）は、ADC の前のアナログのローパスで止めておきます。

設計例 3：分数遅延とハーフバンド

要件 1： fs = 8 kHz の 2 つのセンサのうち、B が A より 0.3 サンプル遅れてサンプルされている。A を 0.3 サンプル遅らせて B と時刻をそろえたい。ずらした後の誤差を 1 %（−40 dB）以下に保てる帯域を、できるだけ広く取りたい。

方式： A に Thiran 3 次（D = 3.3）を通し、B を 3 サンプル遅らせます。比較のために、ラグランジュ補間 3 次（D = 1.3、B は 1 サンプル遅らせる）と線形補間（D = 0.3）も計算しました。誤差は、理想の遅れ e^{-j\theta D} との差の大きさ \lvert H(e^{j\theta}) - e^{-j\theta D}\rvert で測ります。この値は、その周波数の正弦波をずらしたときに残る誤差の振幅（入力の振幅に対する比）です。

(a) の −40 dB を横切る周波数（誤差 1 % の帯域）は、Thiran 3 次で 1525 Hz（0.19 fs）、ラグランジュ 3 次で 1097 Hz、線形補間で 394 Hz でした。1 kHz での誤差は、それぞれ −63.6 dB、−43.1 dB、−23.9 dB です。

Thiran は掛け算が多く遅れも 3.3 サンプル（0.41 ms）と長いかわりに、帯域が最も広く、振幅も落ちません。センサの帯域が 400 Hz 以下なら線形補間で足ります。

要件 2： その 8 kHz の列を 4 kHz に落としたい。0〜1.6 kHz を残し、4 kHz に落としたときに 0〜1.6 kHz へ折り返す 2.4〜4 kHz を消したい。

方式： 2 分の 1 なので、ハーフバンドを使います。オールパスの節で求めた 5 次の IIR（c_0 = 0.236471、c_1 = 0.714542）と、等リプル（阻止域の波打ちの高さがそろうように係数を決める設計。SciPy の remez）で設計した 31 タップの FIR ハーフバンドを比べます。FIR は中央が 0.5、それ以外の偶数番目の係数が 0 になり、0 でない係数は 17 個です。

同じ −36 dB の阻止域を FIR ハーフバンドで作ると、19 タップ（掛け算 5 回）が要りました。掛け算の数と遅れでは IIR が有利で、波形の形（位相の直線性）を保つ点では FIR が有利です。

破る例：遷移帯は fs/4 をまたぐ： 表の最後の行のとおり、ハーフバンドは係数をどう選んでも fs/4 で −3 dB（IIR）か −6 dB（FIR）になります。FIR のハーフバンドは振幅が fs/4 を中心に対称（ H(\theta) + H(\pi - \theta) = 1 ）、IIR はパワー相補で fs/4 が半分の電力だからです。要件が「0〜1.8 kHz を残し、2.0 kHz から消したい」のように遷移帯が fs/4（2 kHz）の片側に寄っていると、ハーフバンドでは満たせません。2.0〜2.4 kHz の成分は 4 kHz に落とすと 1.6〜2.0 kHz に折り返します。ハーフバンドが使えるのは、この区間を捨ててよいとき、つまり欲しい帯域が fs,out/2 より遷移帯の幅の分だけ低いときです。

✅ 検証

記事の Python の関数と C のコードを、独立な計算と突き合わせました。乱数のシードは 20261105 で、図ごとに分けています。

CIC：移動平均の N 乗との一致とビット成長

−128〜127 の乱数 4096 個を入力に、CIC の出力と、長さ RM の矩形を N 回畳み込んだ FIR を通して R 個に 1 個を残したものを比べました。R = 4・M = 1・N = 1、R = 8・M = 1・N = 3、R = 16・M = 2・N = 4、R = 256・M = 1・N = 3 の 4 通りすべてで、差は整数で 0 でした。

ビット成長の式も、入力を最小値（B_in = 8 ビットの −128）に固定して確かめました。R = 8・M = 1・N = 3 では出力が −65536 = -2^{16} で、式の B_out = 17 ビットの最小値に一致します。R = 16・M = 2・N = 4 では -2^{27} で、B_out = 28 ビットに一致しました。補間用の CIC の DC のゲインは、R = 4・M = 1・N = 3 で 16（式は 64/4）、R = 8・M = 1・N = 2 で 8（式は 64/8）でした。

オールパス、Thiran、ハーフバンド

K = 1 の Thiran は、0.05 fs ですでに群遅延が 0.7088 と目標からずれ始めています。次数を上げるほど平らな帯域が広がることが、ここからも分かります。

SciPy との比較

SciPy 1.17.1 の decimate は、FIR を選ぶとハミング窓の firwin で長さ 20q+1 の FIR を作り、zero_phase=False なら lfilter を通して q 個に 1 個を残します。resample_poly の既定の窓は Kaiser 窓（β で形を調整できる窓関数。既定は β = 5.0）です。どちらも、記事の自前の関数と倍精度の丸めの範囲で一致しました。

シードを変えても

設計例 1（2 次の ΔΣ、R = 256、N = 3、4096 出力）の SNR は、シード 20261105・1・2 で 102.08 dB、102.01 dB、101.69 dB でした。乱数は入力の初期位相と小さなノイズだけなので、ばらつきは 0.4 dB 以内です。

C と Python

C（float と uint32_t）を MSVC でビルドし、Python の出力と比べました。CIC（R = 256、N = 3、25 ビット）は 300 個の出力がすべて整数で一致しました。float の 1 次オールパス、Thiran 3 次、ハーフバンドは、標準正規の乱数 4000 個に対して、差の最大がそれぞれ 3.6 × 10⁻⁷、4.3 × 10⁻⁷、2.6 × 10⁻⁷ でした（出力の最大は 2.8〜4.2 で、float の丸めの大きさ）。

チェックリスト

1. 間引きの前に、fs,out/2 より上を消したか：折り返しの減衰を、欲しい帯域の幅 b から計算する
2. 補間のローパスのゲインは L か：0 を挟むと振幅が 1/L になる
3. CIC の語幅は B_in + ⌈N log2(RM)⌉ 以上か：途中のあふれは巻き戻しでよい。飽和させない。C では符号なしの型
4. CIC の出力を (RM)^N で割ったか：割らないと、後段のゲインが数万倍になる
5. ドループは要件の中か：はみ出すなら補償 FIR。補償 FIR が高域のノイズを持ち上げることも数える
6. ループに入れるなら、遅れを位相余裕へ：群遅延 N(RM−1)/2 と、出力のタイミングのずれ（最大 1 出力サンプル）

🕳️ 落とし穴

フィルタなしで間引く

いちばん多いのは、サンプルを単に間引くことです。はじめにの図 (a) のとおり、900 Hz が 100 Hz に化け、後から消す方法はありません。マイコンで「重いので 10 回に 1 回だけ使う」も同じ間引きです。間引くなら、少なくとも R 個の平均（N = 1 の CIC）を取り、帯域が広いなら段数を増やします。

補間のゲインを忘れる

0 を挟んだだけでは振幅が 1/L になります。補間のローパスのゲインを L にするのを忘れると、出力が小さくなります。逆に、0 ではなく同じ値を L 回くり返す（ゼロ次ホールド）と、振幅は保たれますが、sinc 状の特性が掛かり、イメージも残ります。

CIC の語幅が足りない、飽和させる

語幅を B_out より小さくすると、検証の図の 23 ビットのように出力が壊れます。また、DSP（信号処理用のプロセッサ）やマイコンの飽和演算（あふれたら最大値で止める）を積分器に使うと、巻き戻しが成り立たず、やはり壊れます。CIC のレジスタは、巻き戻る足し算で作ります。

CIC の DC ゲインのスケーリング

CIC の出力は (RM)^N 倍です。R = 256、N = 3 なら約 1677 万倍です。2 のべき乗なら右シフトで割れますが、R = 100 のように 2 のべき乗でないと、掛け算 1 回（逆数を掛ける）が要ります。割る前に後段へ渡すと、補償 FIR やループのゲインが桁違いになります。

ドループを放置する

sinc^N のドループは、fs,out の 25 % で N = 3 なら −2.7 dB です。センサの振幅を測る用途では、この分だけ高い周波数の成分が小さく出ます。帯域の端まで振幅を使うなら、補償 FIR を入れるか、ドループを読み取り値の補正に入れます。

非整数の比を「近い整数」で済ませる

44.1 kHz を 48 kHz にするとき、比 160/147 を「だいたい 1.09 倍」として整数の操作で近似すると、ずれが周期的にたまり、ときどきサンプルを飛ばしたり重ねたりすることになります。これは周期的な位相のぶれ（ジッタ）になり、スペクトルに細い線が立ちます。有理数比なら L/R の補間・間引きで正確に作ります。

ハーフバンドの遷移帯と fs/4

設計例 3 のとおり、ハーフバンドは fs/4 で必ず −3 dB か −6 dB になります。欲しい帯域が fs/4 に近いと、通過域の端が削れるか、折り返しが残ります。fs/4 の近くまで使いたいなら、ハーフバンドではない一般のローパスを使います。

Thiran の D が範囲の外

D が K - 1 以下だとフィルタは発散し、 K - 0.5 \le D \le K + 0.5 から外れると平らな帯域が狭くなります。遅れを時間とともに変える用途では、D が範囲から出ないように整数部を切り替える処理が要り、切り替えの瞬間に過渡が出ます。

制御ループの中の間引きの遅れ

制御工学の道具の節のとおり、同じ −3 dB でも段数の多い CIC ほど遅れが長く、位相余裕が減ります。ループの帯域を決めるのは別のフィルタに任せ、CIC は M = 1 で折り返し止めに使うと、遅れは小さく済みます。

任意の比、変わり続ける比

2 つの独立したクロックの間のように、比が有理数で書けない、あるいは時間とともにゆっくり変わる変換（非同期のレート変換）は、この記事の範囲の外です。分数遅延の D を毎サンプル変える必要があり、そのために係数を D の多項式で表した FIR（Farrow 構造、Farrow 1988）がよく使われます。ここでは名前だけにとどめます。

💻 実装

C では、CIC 間引き（uint32_t、N = 3）、1 次オールパス、Thiran 3 次、ポリフェーズの IIR ハーフバンドを、#1 の biquad と同じ「構造体に状態を持たせ、1 サンプルずつ回す関数」の流儀で書きます。Python では、記事の数値を出した関数そのものを載せます。

C：CIC・オールパス・Thiran・ハーフバンド

ヘッダです。CIC の段数と差分遅延はコンパイル時の定数にしています。

/* rate.h : CIC 間引き・1 次オールパス・Thiran 3 次・ポリフェーズ IIR ハーフバンド */
#ifndef RATE_H
#define RATE_H
#include <stdint.h>

#define CIC_N 3          /* 段数 */
#define CIC_M 1          /* 差分遅延 */

/* CIC 間引き。レジスタは uint32_t（符号なしの巻き戻しは C で定義された動作） */
typedef struct {
    uint32_t integ[CIC_N];         /* 積分器（入力レート） */
    uint32_t comb_z[CIC_N][CIC_M]; /* コムの遅延（出力レート） */
    int      R;                    /* 間引き比 */
    int      count;                /* 0 .. R-1 */
} cic_t;

void cic_init(cic_t *c, int R);
/* 入力 1 サンプルを積む。出力の番（R 個目）なら 1 を返し *out に入れる */
int  cic_push(cic_t *c, uint32_t x, uint32_t *out);

/* 1 次オールパス A(z) = (a + z^-1)/(1 + a z^-1) */
typedef struct { float a, x1, y1; } ap1_t;
void  ap1_init(ap1_t *s, float a);
float ap1_step(ap1_t *s, float x);

/* Thiran 3 次（分母 1 + a1 z^-1 + a2 z^-2 + a3 z^-3、分子は逆順） */
typedef struct { float a[4]; float x[3], y[3]; } thiran3_t;
void  thiran3_init(thiran3_t *s, float D);
float thiran3_step(thiran3_t *s, float x);

/* 2 本のオールパスのハーフバンドで 2 分の 1 に間引く（入力 2 個で出力 1 個） */
typedef struct { ap1_t p0, p1; float xodd_prev; } hb_t;
void  hb_init(hb_t *s, float c0, float c1);
float hb_step2(hb_t *s, float x_even, float x_odd);
#endif

本体です。CIC の積分器とコムは、足し算と引き算だけです。

#include "rate.h"

void cic_init(cic_t *c, int R)
{
    for (int i = 0; i < CIC_N; i++) {
        c->integ[i] = 0;
        for (int k = 0; k < CIC_M; k++) c->comb_z[i][k] = 0;
    }
    c->R = R;
    c->count = 0;
}

int cic_push(cic_t *c, uint32_t x, uint32_t *out)
{
    uint32_t v = x;
    for (int i = 0; i < CIC_N; i++) {      /* 積分器：足すだけ。あふれても巻き戻る */
        c->integ[i] += v;
        v = c->integ[i];
    }
    if (++c->count < c->R) return 0;
    c->count = 0;
    for (int i = 0; i < CIC_N; i++) {      /* コム：M 個前を引く */
        uint32_t d = c->comb_z[i][CIC_M - 1];
        for (int k = CIC_M - 1; k > 0; k--) c->comb_z[i][k] = c->comb_z[i][k - 1];
        c->comb_z[i][0] = v;
        v = v - d;
    }
    *out = v;                              /* 正しいのは下位 B_out ビット */
    return 1;
}

void ap1_init(ap1_t *s, float a) { s->a = a; s->x1 = 0.0f; s->y1 = 0.0f; }

float ap1_step(ap1_t *s, float x)
{
    float y = s->a * x + s->x1 - s->a * s->y1;
    s->x1 = x;
    s->y1 = y;
    return y;
}

void thiran3_init(thiran3_t *s, float D)
{
    /* a_k = (-1)^k C(3,k) prod_{n=0}^{3} (D-3+n)/(D-3+k+n) */
    static const float C3[4] = {1.0f, 3.0f, 3.0f, 1.0f};
    s->a[0] = 1.0f;
    for (int k = 1; k <= 3; k++) {
        float p = 1.0f;
        for (int n = 0; n <= 3; n++) p *= (D - 3.0f + n) / (D - 3.0f + k + n);
        s->a[k] = ((k & 1) ? -1.0f : 1.0f) * C3[k] * p;
    }
    for (int i = 0; i < 3; i++) { s->x[i] = 0.0f; s->y[i] = 0.0f; }
}

float thiran3_step(thiran3_t *s, float x)
{
    const float *a = s->a;   /* 分子は a3 + a2 z^-1 + a1 z^-2 + z^-3 */
    float y = a[3] * x + a[2] * s->x[0] + a[1] * s->x[1] + s->x[2]
            - a[1] * s->y[0] - a[2] * s->y[1] - a[3] * s->y[2];
    s->x[2] = s->x[1]; s->x[1] = s->x[0]; s->x[0] = x;
    s->y[2] = s->y[1]; s->y[1] = s->y[0]; s->y[0] = y;
    return y;
}

void hb_init(hb_t *s, float c0, float c1)
{
    ap1_init(&s->p0, c0);
    ap1_init(&s->p1, c1);
    s->xodd_prev = 0.0f;
}

float hb_step2(hb_t *s, float x_even, float x_odd)
{
    /* 偶数番目は A0 へ、奇数番目は 1 つ遅らせて A1 へ（どちらも出力レートで回す） */
    float y = 0.5f * (ap1_step(&s->p0, x_even) + ap1_step(&s->p1, s->xodd_prev));
    s->xodd_prev = x_odd;
    return y;
}

使い方の例です。ΔΣ の 1 bit を 1 個ずつ渡し、出力の番が来たら下位 25 ビットを取り出して 2^{24} で割ります。

cic_t cic;
cic_init(&cic, 256);
uint32_t raw;
if (cic_push(&cic, bit, &raw)) {                  /* bit は 0 か 1 */
    uint32_t v = raw & 0x1FFFFFFu;                /* 下位 25 ビット（0 .. 2^24） */
    float density = (float)v / 16777216.0f;       /* 1 の割合（0 .. 1） */
    float x = 2.0f * density - 1.0f;              /* ±1 の目盛りに戻す */
}

hb_t hb;
hb_init(&hb, 0.23647102f, 0.71454215f);           /* 5 次ハーフバンド */
float y = hb_step2(&hb, x0, x1);                  /* 8 kHz の 2 個 → 4 kHz の 1 個 */

入力が ±1 の符号付きなら、レジスタは uint32_t のまま、入力を 2 の補数のビット列として uint32_t に入れ、最後に下位 B_out ビットを符号拡張（最上位ビットが 1 なら上位を 1 で埋める）して int32_t に戻します。B_out が 32 を超えるときは uint64_t にします。

Python：記事の関数

記事の図と数値は、次の関数で計算しています。NumPy だけで動きます。

import numpy as np
from math import comb as binom


def wrap(v, width, signed=True):
    """width ビットのレジスタに巻き戻す（None なら何もしない）。"""
    if width is None:
        return v
    m = 1 << width
    if not signed:
        return v % m
    return ((v + (m >> 1)) % m) - (m >> 1)


def cic_decimate(x, R, M=1, N=3, width=None, signed=True):
    """整数の CIC 間引き。出力は n = R-1, 2R-1, ... の位置、DC ゲインは (RM)^N。"""
    y = np.asarray(x, dtype=np.int64)
    for _ in range(N):                      # 積分器（入力レート）
        y = wrap(np.cumsum(y), width, signed)
    y = y[R - 1::R]                         # R 個に 1 個を残す
    for _ in range(N):                      # コム（出力レート、遅延 M）
        d = np.zeros_like(y)
        d[M:] = y[:-M]
        y = wrap(y - d, width, signed)
    return y


def comb(x, M):
    """コム y[n] = x[n] - x[n-M]。"""
    x = np.asarray(x, float)
    y = x.copy()
    y[M:] -= x[:-M]
    return y


def allpass1(x, a, x1=0.0, y1=0.0):
    """1 次オールパス (a + z^-1)/(1 + a z^-1)。"""
    y = np.zeros(len(x))
    for n, v in enumerate(x):
        y[n] = a * v + x1 - a * y1
        x1, y1 = v, y[n]
    return y


def thiran_coeffs(D, K):
    """Thiran の K 次オールパスの分母の係数 a_0..a_K（a_0 = 1）。"""
    a = np.ones(K + 1)
    for k in range(1, K + 1):
        p = 1.0
        for n in range(K + 1):
            p *= (D - K + n) / (D - K + k + n)
        a[k] = (-1) ** k * binom(K, k) * p
    return a


def halfband_iir(x, c0, c1):
    """(A0(z^2) + z^-1 A1(z^2))/2 で 2 分の 1 に間引く（ポリフェーズ）。"""
    x = np.asarray(x, float)
    xe = x[0::2]
    xo = np.concatenate(([0.0], x[1::2]))[:len(xe)]
    return 0.5 * (allpass1(xe, c0) + allpass1(xo, c1))


def polyphase_decimate(x, h, R):
    """FIR h で R 分の 1 に間引く。出力は「全部フィルタして n = mR を残す」と同じ。"""
    x = np.asarray(x, float)
    h = np.asarray(h, float)
    L = -(-len(h) // R) * R
    hp = np.zeros(L)
    hp[:len(h)] = h
    nout = -(-len(x) // R)
    y = np.zeros(nout)
    for k in range(R):
        Ek = hp[k::R]                       # E_k の係数 h[k], h[k+R], ...
        idx = np.arange(nout) * R - k       # 入力の相 k：x[mR - k]
        xk = np.zeros(nout)
        ok = idx >= 0
        xk[ok] = x[idx[ok]]
        y += np.convolve(xk, Ek)[:nout]
    return y

Thiran のフィルタそのものは、分母 a と、それを逆に並べた分子で scipy.signal.lfilter(a[::-1], a, x) と書けます。周波数応答と群遅延は scipy.signal.freqz と scipy.signal.group_delay で確かめられます。

手計算を確かめる

本文の 4 つの手計算を、上の関数で確かめるコードです。

import numpy as np

# 1) ポリフェーズ：長さ 8 の FIR を R = 2 で
h = np.array([1, 2, 3, 4, 4, 3, 2, 1.0])
x = np.array([3, 1, 4, 1, 5, 9, 2, 6, 5, 3.0])
print("全部計算して偶数番目:", np.convolve(x, h)[:len(x)][0::2])
print("ポリフェーズ        :", polyphase_decimate(x, h, 2))

# 2) CIC：R = 4、M = 1、N = 1
xi = np.array([3, 1, 4, 1, 5, 9, 2, 6])
print("積分器:", np.cumsum(xi), " CIC:", cic_decimate(xi, 4, 1, 1))

# 3) 1 次オールパス a = 0.5、入力 cos(pi n / 2) の定常状態
y = allpass1(np.tile([1, 0, -1, 0.0], 50), 0.5)[-4:]
Y = np.fft.fft(y)
print("y:", y, " Y[1]:", Y[1], " |Y[1]|:", abs(Y[1]),
      " 偏角[deg]:", np.degrees(np.angle(Y[1])))

# 4) Thiran 3 次 D = 3.3
print("Thiran:", thiran_coeffs(3.3, 3))

出力は次のとおりです。

全部計算して偶数番目: [ 3. 15. 35. 64. 91.]
ポリフェーズ        : [ 3. 15. 35. 64. 91.]
積分器: [ 3  4  8  9 14 23 25 31]  CIC: [ 9 22]
y: [ 0.8  0.6 -0.8 -0.6]  Y[1]: (1.6-1.2j)  |Y[1]|: 2.0  偏角[deg]: -36.86989764584402
Thiran: [ 1.         -0.20930233  0.05133831 -0.00624752]

4 つとも、本文の手計算（91、9 と 22、0.8・0.6・−0.8・−0.6 と −36.87°、Thiran の係数）と一致します。

ライブラリ

• SciPy：scipy.signal.decimate（間引き。既定は 8 次のチェビシェフ I 型の IIR（通過域に小さな波打ちを許して肩を急にした IIR）で、FIR も選べる）、scipy.signal.resample_poly（有理数比 up/down をポリフェーズで）、scipy.signal.upfirdn（0 を挟む → FIR → 間引きを 1 つの関数で）があります。CIC の専用の関数はないので、上の cic_decimate か、矩形の FIR を upfirdn に渡して確かめます
• Arm CMSIS-DSP：FIR の間引き arm_fir_decimate_f32 と、FIR の補間 arm_fir_interpolate_f32 があります（q15、q31 版もあります）。リポジトリに CIC の関数は見当たりません
• Espressif ESP-DSP：FIR の間引き dsps_fird_f32（16 ビット固定小数点の dsps_fird_s16 も）と、ポリフェーズの再サンプリング（dsps_resampler_ph_*）、多レートの再サンプリング（dsps_resampler_mr_*）があります。こちらも CIC の関数は見当たりません

💡 活用事例

レート変換のフィルタが使われている場面を、1〜2 文ずつ挙げます。

• ΔΣ 型の ADC と DAC：ADC は変調器の 1 bit（または数ビット）の列を sinc フィルタ（CIC）と FIR で間引いて多ビットにする。DAC は逆に、補間フィルタでレートを上げてから ΔΣ 変調器に渡す
• PDM 出力の MEMS マイク：MEMS（半導体の工程で作る微小な機械構造）のマイクが 1 bit の列（PDM：パルス密度変調。1 の密度で値を表す方式）を出し、受け側のマイコンや DSP が間引きのフィルタで PCM（各サンプルを多ビットの数で表す、ふつうのデジタル音声の形式）の音声にする。CIC と FIR を組み合わせた構成がよく使われる
• SDR（ソフトウェア無線。復調などの処理をデジタルで行う無線機）：高速な ADC の出力を、FPGA の中で CIC で大きく間引き、補償 FIR とハーフバンドで仕上げる。掛け算なしで R を大きく取れる CIC の利点がいちばん効く場面
• オーディオのサンプリング周波数の変換：44.1 kHz と 48 kHz の間の変換は、比が 147/160（または 160/147）の有理数比になる
• ロードセル、熱電対、ブリッジセンサ：10 SPS の sinc フィルタで 50 Hz と 60 Hz を同時に消す。設計例 2 の構成
• センサの時刻合わせ：複数のセンサやマイクのサンプル時刻のずれを、分数遅延で合わせる。センサ融合や、マイクを並べて音の方向を推定する処理で使う
• 制御の多レート構成：内側の電流ループは速く、外側の速度や位置のループは遅く回す。速いレートのセンサ値を遅いループに渡すところで、CIC などの折り返し止めを入れる

🧭 使い分けの勘所：CIC・sinc・ハーフバンド・オールパス

ここまでの結果を、判断の材料としてまとめます。方式を選ぶ前に、比と帯域と計算の制約を確かめます。

判断の軸

方式を決める前に答えておく 7 つの問いです。上から順に確かめます。

④が抜けがちです。CIC の折り返しの減衰は、段数よりも先に、帯域の幅 b で決まります。b = 0.4 では N = 3 でも 17.8 dB しか止まりません。

方式の比較

間引き・補間の側の 3 つです。

遅れと 50/60 Hz の側の 3 つです。

決定木

同じ判断を 1 本の流れにします。上から順に答えていくと、方式と最初に決める値に着きます。

ハーフバンドを選んだら、次に位相の直線性が要るかで FIR か IIR かを決めます。CIC を選んだら、帯域の幅 b から段数 N を、B_in と N log2(RM) から語幅を決めます。

用途別の早見表

よくある場面ごとに、最初に試す方式と値の目安を並べます。

勘所

• 間引きの前に必ずフィルタ：折り返した成分は、後からどんなフィルタでも本物と区別できない。最低でも R 個の平均を取る
• R が大きいなら CIC と短い補償 FIR：掛け算なしで R を大きく取り、ドループだけを短い FIR で直す。帯域が広いなら CIC で 2〜4 倍のレートまでにとどめ、残りを FIR で落とす
• CIC の語幅は式で決め、巻き戻しに任せる： B_{in} + \lceil N\log_2(RM)\rceil ビット。途中のあふれは気にしない。飽和演算と符号付きのあふれは使わない
• 段数は「ノイズの増え方」に勝つ数：2 次の ΔΣ なら N = 3。同じ −3 dB なら、段数が多いほど遅れが長い
• 10 SPS は 50 と 60 を同時に消す：零点が 10 Hz の整数倍。段数を増やすと、周波数のずれに強くなる代わりに整定が遅くなる
• 2 倍ならハーフバンド：fs/4 で −3 dB（IIR）か −6 dB（FIR）になることを前提に、遷移帯を fs/4 の両側に置く
• 位相だけ動かすならオールパス：振幅を変えずに、整数でない遅れ（Thiran、D は K の近く）や 90° の位相差を作る
• ループの中では遅れで選ぶ：CIC は M = 1 で折り返し止めに使い、帯域を決めるのは別のフィルタに任せる

✅ まとめ

• レート変換は、fs を変える前（間引き）か後（補間）に、必ずローパスが要る。間引きは R 枚のスペクトルが重なって折り返し、補間は L−1 個のイメージが出る
• ノーブル恒等式で、z^R の形のフィルタは低いレートへ移せる。ポリフェーズはそれを一般の FIR に広げ、捨てる出力を計算しない
• CIC は長さ RM の移動和の N 乗で、掛け算なしで作れる。零点が fs,out/M の整数倍、つまり折り返してくる帯の真ん中に来る。帯域の幅 b が小さいほど折り返しは深く消え、b = 0.4 では N = 3 でも 17.8 dB にとどまる
• CIC の語幅は B_in + ⌈N log2(RM)⌉ ビットで、途中の積分器があふれても巻き戻しで正しい出力になる。R = 256、N = 3 の 1 bit 入力なら 25 ビット
• ΔΣ ADC の sinc フィルタは CIC そのもので、10 SPS なら 50 Hz と 60 Hz に同時に零点が来る。±2 % のずれで N = 1 は 34 dB、N = 4 は 136 dB
• オールパスは振幅 1 のまま位相だけを動かす。Thiran の係数で整数でない遅れを作り（D = 3.3 で誤差 1 % の帯域 0.19 fs）、2 本の組でハーフバンド（掛け算 2 回で −36 dB）と 90° の位相差（掛け算 4 回で ±0.4°）を作れる
• 制御ループの中では、同じ −3 dB でも段数の多い CIC ほど遅れが長い。CIC は M = 1 で折り返し止めに使うと、位相余裕の損は 4.5° で済み、1010 Hz の振動が 10 Hz に化けてモータを揺らすのを止められる

ここまでで、比と帯域の幅と計算の制約から方式を選び、CIC の段数と語幅、補償 FIR、sinc フィルタの出力レートと段数、Thiran の次数と D、ハーフバンドの係数を、式と数値から決められるようになりました。次の一手は、手元の ΔΣ ADC（ロードセル用や温度計測用）のデータシートを開いて、デジタルフィルタの段数（sinc3 か sinc4 か）と出力レートを読み、この記事の表で 50/60 Hz の減衰と整定時間を見積もってみることです。

よくある質問（FAQ）

Q. 平均を取ってから間引くのでは、だめなのですか？

A. R 個の平均を取ってから間引くのは、N = 1 の CIC そのものなので、何もしないよりずっと良い方法です。ただし、平均の零点のまわりの谷は浅く、帯域の幅 b = 0.05 で 25.6 dB、b = 0.4 で 5.9 dB しか折り返しを止められません。2 次の ΔΣ の 1 bit の列を R = 256 で間引いた例では、平均（N = 1）の SNR は 39 dB、N = 3 の CIC は 102 dB でした。平均を 3 回重ねる（CIC の N = 3）だけで、計算はほとんど増えずに大きく改善します。

Q. CIC の積分器があふれても、本当に大丈夫ですか？

A. 全段を同じ語幅（B_in + ⌈N log2(RM)⌉ ビット以上）で作り、あふれたときに巻き戻る足し算（2 の補数か符号なし）を使えば、出力は正しくなります。足し算と引き算だけの計算は、語幅を法とする余りの計算として正しく、本当の出力がその語幅に収まるからです。だめなのは、飽和演算を使う場合と、C の符号付き整数のあふれ（未定義動作）に頼る場合です。C では uint32_t か uint64_t で作ります。

Q. sinc3 と sinc4 は、どちらを選べばいいですか？

A. 50/60 Hz を消すのが目的なら、電源とクロックのずれの見込みで決めます。10 SPS で ±2 % のずれなら、sinc3 で約 102 dB、sinc4 で約 136 dB です。代わりに、入力が変わってから整定するまでの出力の数が、sinc3 は 4 個、sinc4 は 5 個で、10 SPS ならそれぞれ 0.4 秒と 0.5 秒です。−3 dB の帯域も sinc4 のほうが少し狭くなります。

Q. 2 倍の間引きには、IIR と FIR のどちらのハーフバンドを使えばいいですか？

A. 軽さと遅れの短さを取るなら IIR、波形の形を保つなら FIR です。5 次の IIR ハーフバンドは 1 出力あたり掛け算 2 回で −36 dB、群遅延は 0.16〜0.41 ms と短いですが、周波数で変わります。同じ −36 dB の FIR は 19 タップ（掛け算 5 回）で、群遅延は一定です。音声や振動の解析で波形の形を比べるなら FIR、包絡線やレベルだけを見るなら IIR が向いています。

Q. 2 つのセンサのサンプル時刻のずれは、どこまで分数遅延で直せますか？

A. ずれ自体は何サンプルでも直せますが、正確に直せる帯域に限りがあります。0.3 サンプルを Thiran 3 次で直す例では、誤差 1 % 以下の帯域は fs の 0.19 倍（8 kHz で 1525 Hz）、線形補間では 394 Hz でした。信号の帯域が fs の 5 % 程度までなら線形補間で足り、それより広ければ Thiran やラグランジュの次数を上げます。ずれが時間とともに変わるなら、Farrow 構造のような可変の分数遅延が要ります。

関連記事

この連載：

• #1 基礎編 — サンプリング定理と折り返し、群遅延、C の biquad
• #2 ローパスフィルタの設計 — 移動平均（＝コムと積分器）と、制御ループの位相余裕
• #3 ハイパスフィルタの設計 — DC とドリフトを消す
• #4 バンドパスとノッチの設計 — 1 つの周波数を消すノッチ（コムは等間隔のノッチ）
• #5 カルマンフィルタの設計 — モデルと分散で切り分ける
• #6 EKF・UKF の設計 — 線形でない対象へ広げる
• #7 パーティクルフィルタの設計 — 分布をサンプルで持つ
• #8 適応フィルタの設計 — 入力と誤差から係数を学習する
• #9 平滑化 — 外れ値を落とし、形を保つ
• #10 検出 — 周波数・波形・包絡線を取り出す。FIR のヒルベルト変換

信号の入口：

• Arduino の AD 変換（ADC） — サンプルを取り込む入口。間引く前の平均の取り方
• 加速度センサ KXR94-2050 で作るデジタル水平器 — センサの値を読んで処理する題材

参考

取得日はいずれも 2026-11-10 です。図と数値はすべて筆者が Python（NumPy、SciPy、Matplotlib）で計算したものです。乱数のシードは固定です。

CIC と多レート信号処理

• E. B. Hogenauer「An economical class of digital filters for decimation and interpolation」（IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(2), 155–162, 1981）
• R. E. Crochiere, L. R. Rabiner「Interpolation and decimation of digital signals — A tutorial review」（Proceedings of the IEEE, 69(3), 300–331, 1981）
• M. G. Bellanger, G. Bonnerot, M. Coudreuse「Digital filtering by polyphase network: Application to sample-rate alteration and filter banks」（IEEE Transactions on Acoustics, Speech, and Signal Processing, 24(2), 109–114, 1976）
• P. P. Vaidyanathan「Multirate digital filters, filter banks, polyphase networks, and applications: A tutorial」（Proceedings of the IEEE, 78(1), 56–93, 1990）

オールパスと分数遅延

• P. A. Regalia, S. K. Mitra, P. P. Vaidyanathan「The digital all-pass filter: A versatile signal processing building block」（Proceedings of the IEEE, 76(1), 19–37, 1988）
• J.-P. Thiran「Recursive digital filters with maximally flat group delay」（IEEE Transactions on Circuit Theory, 18(6), 659–664, 1971）
• T. I. Laakso, V. Välimäki, M. Karjalainen, U. K. Laine「Splitting the unit delay: Tools for fractional delay filter design」（IEEE Signal Processing Magazine, 13(1), 30–60, 1996）
• C. W. Farrow「A continuously variable digital delay element」（IEEE International Symposium on Circuits and Systems, 2641–2645, 1988）

ΔΣ ADC

• Texas Instruments「ADS1232, ADS1234 24-Bit Analog-to-Digital Converters for Bridge Sensors」（SBAS350H, 2025 年 6 月改訂）

教科書

• R. E. Crochiere, L. R. Rabiner「Multirate Digital Signal Processing」（Prentice Hall, 1983）
• A. V. Oppenheim, R. W. Schafer「Discrete-Time Signal Processing」（Pearson）
• f. j. harris「Multirate Signal Processing for Communication Systems」（Prentice Hall）

ライブラリ

• SciPy「scipy.signal.decimate」
• SciPy「scipy.signal.resample_poly」
• SciPy「scipy.signal.upfirdn」
• Arm CMSIS-DSP「arm_fir_decimate_f32.c」
• Arm CMSIS-DSP「arm_fir_interpolate_f32.c」
• Espressif ESP-DSP「dsps_fir.h」
• Espressif ESP-DSP「dsps_resampler.h」