ソフトウェア無線（SDR）や $\Delta\Sigma$ 型ADCの出力は、しばしば数十MHz〜数百MHzという非常に高いサンプリングレートで生成されます。この超高速のストリームをそのままFPGAやマイコンで一般的なFIRフィルタにかけて間引こうとすると、1サンプルあたり何十回もの乗算が必要になり、回路規模も消費電力も現実的でなくなります。「乗算器をひとつも使わずに、加算と遅延だけで高効率に間引きできるフィルタはないのか？」という問いに対する、エレガントで実用的な答えが CICフィルタ（Cascaded Integrator-Comb filter、縦続積分器・くし形フィルタ） です。
CICフィルタは、Eugene Hogenauer が1981年に提案したことから Hogenauer フィルタとも呼ばれ、今日では SDR受信機のフロントエンド、$\Delta\Sigma$ ADC/DACのデシメーション・インターポレーション段、レーダーやデジタル無線機のサンプルレート変換など、高速ストリーム処理のいたるところで使われています。本記事を読むと、なぜ「積分器」と「くし形フィルタ」を縦続するだけで移動平均フィルタになるのか、その周波数特性がなぜ $|\mathrm{sinc}|^N$ の形になるのか、そして実用上避けて通れない「通過帯域のドループ（垂れ下がり）」をどう補償するのかが、数式の導出とPython実装の両面からすっきり理解できるようになります。
本記事の内容
- 乗算器なしで動くCICフィルタの直感的な仕組み（積分器＋くし形フィルタ＝移動平均）
- 積分器 $1/(1-z^{-1})$ とくし形 $1-z^{-M}$ の伝達関数の導出
- 縦続が矩形窓（移動平均）と等価になることの証明
- 周波数応答が $|\mathrm{sinc}|^N$ 型になることの導出とエイリアス帯の減衰設計
- デシメーションとの結合（Noble恒等式による積分器の高速側配置）
- 固定小数点実装でのビット成長（bit growth）の見積もり
- 通過帯域ドループとFIRコンペンセータ（CICコンペンセータ）による補償
- Pythonでのintegrator/comb構造による実装、周波数特性、デシメーション動作、補償前後の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、z変換と伝達関数 $H(z)$、単位円上での周波数応答 $H(e^{j\omega})$、そしてデシメーション（間引き）の概念は前提として用います。CICフィルタは積分器（IIR的な再帰構造）とくし形（FIR的な構造）を組み合わせた、両者の性質を併せ持つフィルタなので、FIR/IIRの基礎を押さえておくと理解が格段に速くなります。
CICフィルタとは
レートの高い信号を「間引いて」レートを下げるとき、いきなり1個おきに値を捨てると（単純デシメーション）、ナイキスト周波数を超える成分が低い周波数に折り返してくる「エイリアシング」が起きてしまいます。これを防ぐには、間引く前に高周波成分を落とす低域通過フィルタ（アンチエイリアシングフィルタ）が必要です。ところが入力レートが数百MHzともなると、普通のFIRフィルタでは乗算回数が膨大になり、ハードウェアに載りません。
ここで発想を変えてみましょう。最も単純な低域通過フィルタは何でしょうか？ それは 移動平均 です。直近 $RM$ 個のサンプルの平均を取るだけで、細かい変動（高周波）はならされ、ゆっくりした変動（低周波）だけが残ります。移動平均は「全タップ係数が等しい矩形窓FIRフィルタ」であり、係数がすべて $1$ なので、本来なら乗算は一切要りません。ただ素朴に実装すると、窓の長さ $L$ ごとに $L-1$ 回の加算が要ります。
CICフィルタの巧妙さは、この移動平均を 「積分器（累積和を取る）」と「くし形フィルタ（$M$ サンプル前との差分を取る）」の組み合わせ で実現する点にあります。積分器で全部の過去を足し上げ、くし形で「窓の外に出た古い分」を引き算する——こうすると、窓の長さがどれだけ長くても、1サンプルあたりの演算は「1回の加算」と「1回の減算」だけで済みます。窓長に依存しない一定コストで、しかも乗算ゼロ。これがCICフィルタが超高速ストリームで重宝される理由です。
イメージとしては、レジのお釣り計算に似ています。長い買い物リストの合計を毎回ゼロから足し直すのではなく、累積合計をひとつ覚えておき、新しい品物を足し、古くなった品物を引く——差分だけを更新していくのです。
CICフィルタの正式名称 Cascaded Integrator-Comb が示す通り、これを以下のように配置します。
- 積分器（Integrator） を $N$ 段、高いサンプルレート側に縦続する
- 間に レート変更（デシメーションなら $\downarrow R$、インターポレーションなら $\uparrow R$） を挟む
- くし形フィルタ（Comb） を $N$ 段、低いサンプルレート側に縦続する
段数 $N$ を増やすほど、移動平均が何重にもかかり、阻止帯域の減衰が深くなります。次節では、この積分器とくし形フィルタそれぞれの伝達関数を求め、それらの縦続が本当に移動平均になることを数式で確かめていきましょう。
積分器とくし形フィルタの伝達関数
積分器（累積器）
積分器は、入力のこれまでの値をすべて足し上げる累積和器です。差分方程式は次のように書けます。
$$ y_I[n] = y_I[n-1] + x[n] $$
「ひとつ前の出力に、今の入力を足す」という単純な再帰です。両辺をz変換しましょう。$y_I[n-1]$ のz変換は時間シフト性質より $z^{-1} Y_I(z)$ なので
$$ Y_I(z) = z^{-1} Y_I(z) + X(z) $$
$Y_I(z)$ について整理すると
$$ Y_I(z)\left(1 – z^{-1}\right) = X(z) $$
したがって積分器の伝達関数は
$$ H_I(z) = \frac{Y_I(z)}{X(z)} = \frac{1}{1 – z^{-1}} $$
これは極が $z = 1$（単位円上）にある、本質的に不安定すれすれの構造です。実際、積分器単体に有界な入力を入れると出力は際限なく成長しかねません。後で見るように、CICフィルタが安定に動作するのは、この積分器の「無限大に発散しようとする性質」をくし形フィルタが正確に打ち消すからです。
くし形フィルタ
くし形フィルタは、今の入力から「$M$ サンプル前の入力」を引き算するフィルタです。$M$ は 差分遅延（differential delay） と呼ばれ、通常 $M=1$ または $M=2$ を使います。差分方程式は
$$ y_C[n] = x[n] – x[n-M] $$
両辺をz変換すると、$x[n-M]$ は $z^{-M} X(z)$ なので
$$ Y_C(z) = X(z) – z^{-M} X(z) = X(z)\left(1 – z^{-M}\right) $$
くし形フィルタの伝達関数は
$$ H_C(z) = 1 – z^{-M} $$
これは零点を $z^M = 1$、すなわち単位円上に $M$ 等分して並ぶ $M$ 個の点 $z = e^{j 2\pi k / M}\ (k=0,1,\dots,M-1)$ に持ちます。周波数領域で見ると、これらの零点が等間隔の「くし（comb）の歯」のように振幅応答をゼロに落とすことから「くし形」という名前が付いています。$M=1$ のときは零点が $z=1$（直流）のひとつだけになります。
縦続が移動平均になることの導出
ここからが核心です。積分器とくし形フィルタを縦続接続（直列につなぐ）すると何が起きるでしょうか。LTIシステムの縦続では伝達関数が掛け算になるので、$M=R$（差分遅延を間引き率に合わせた標準的なケース）として、両者を掛け合わせます。
$$ H(z) = H_I(z)\, H_C(z) = \frac{1}{1 – z^{-1}} \cdot \left(1 – z^{-RM}\right) = \frac{1 – z^{-RM}}{1 – z^{-1}} $$
ここで、有限等比級数の和の公式
$$ \frac{1 – r^{L}}{1 – r} = 1 + r + r^2 + \cdots + r^{L-1} $$
を $r = z^{-1}$、$L = RM$ として適用します。すると分数が消えて、きれいな多項式になります。
$$ H(z) = \frac{1 – z^{-RM}}{1 – z^{-1}} = \sum_{k=0}^{RM-1} z^{-k} = 1 + z^{-1} + z^{-2} + \cdots + z^{-(RM-1)} $$
この結果が決定的です。伝達関数の係数（インパルス応答）はすべて $1$ で、長さ $RM$ の 矩形窓（boxcar）FIRフィルタ、すなわち（正規化していない）移動平均フィルタそのものになっています。インパルス応答を書き下すと
$$ h[n] = \begin{cases} 1 & 0 \le n \le RM-1 \\ 0 & \text{otherwise} \end{cases} $$
注目すべきは、積分器の極 $z=1$ と、くし形の零点のうちの $z=1$ がぴったり相殺している点です（分子 $1 – z^{-RM}$ も $z=1$ で $0$ になる）。この極零相殺のおかげで、不安定すれすれだった積分器が、全体としては零点しか持たない安定なFIRフィルタへと生まれ変わります。CICフィルタが「IIR的な積分器を含むのにFIR的に安定」という不思議な性質を持つのは、まさにこのためです。
段数 $N$ のCICフィルタは、この基本ブロックを $N$ 段縦続したものなので、伝達関数は
$$ H(z) = \left(\frac{1 – z^{-RM}}{1 – z^{-1}}\right)^{N} = \left(\sum_{k=0}^{RM-1} z^{-k}\right)^{N} $$
となります。これは長さ $RM$ の矩形窓を $N$ 回畳み込んだものに相当し、$N$ が大きいほど滑らかな（ガウス関数に近づく）インパルス応答を持つ低域通過フィルタになります。
ここまでで、CICフィルタの正体が「乗算器なしで実現された $N$ 重の移動平均」であることがわかりました。次は、このフィルタが周波数領域でどんな顔をしているのか——なぜ $\mathrm{sinc}$ 関数が現れるのか——を見ていきましょう。
周波数応答とsinc特性
振幅応答の導出
周波数応答は、伝達関数 $H(z)$ を単位円上 $z = e^{j\omega}$ で評価したものでした。1段あたりの伝達関数 $H(z) = (1 – z^{-RM})/(1 – z^{-1})$ に $z = e^{j\omega}$ を代入します。
$$ H(e^{j\omega}) = \frac{1 – e^{-j\omega RM}}{1 – e^{-j\omega}} $$
分子と分母をそれぞれ「半角に括り出す」テクニックで整理しましょう。一般に $1 – e^{-j\theta}$ は、$e^{-j\theta/2}$ を括り出すと
$$ 1 – e^{-j\theta} = e^{-j\theta/2}\left(e^{j\theta/2} – e^{-j\theta/2}\right) = e^{-j\theta/2}\cdot 2j\sin\!\left(\frac{\theta}{2}\right) $$
と書けます（オイラーの公式 $2j\sin x = e^{jx} – e^{-jx}$ を使いました）。これを分子（$\theta = \omega RM$）と分母（$\theta = \omega$）の両方に適用すると
$$ H(e^{j\omega}) = \frac{e^{-j\omega RM/2}\cdot 2j\sin\!\left(\dfrac{\omega RM}{2}\right)}{e^{-j\omega/2}\cdot 2j\sin\!\left(\dfrac{\omega}{2}\right)} = e^{-j\omega (RM-1)/2}\cdot \frac{\sin\!\left(\dfrac{\omega RM}{2}\right)}{\sin\!\left(\dfrac{\omega}{2}\right)} $$
指数項 $e^{-j\omega(RM-1)/2}$ は純粋な遅延（線形位相）を表し、振幅には寄与しません。したがって振幅応答は
$$ \left|H(e^{j\omega})\right| = \left|\frac{\sin\!\left(\dfrac{\omega RM}{2}\right)}{\sin\!\left(\dfrac{\omega}{2}\right)}\right| $$
となります。これは ディリクレ核（周期的sinc関数） と呼ばれる形です。分母の $\sin(\omega/2)$ が周期性をもたらしますが、$\omega$ が小さい通過帯域付近では $\sin(\omega/2) \approx \omega/2$ と近似でき、
$$ \left|H(e^{j\omega})\right| \approx \left|\frac{\sin\!\left(\dfrac{\omega RM}{2}\right)}{\omega/2}\right| = RM\left|\frac{\sin\!\left(\dfrac{\omega RM}{2}\right)}{\dfrac{\omega RM}{2}}\right| = RM\,\bigl|\mathrm{sinc}\bigl(\tfrac{\omega RM}{2}\bigr)\bigr| $$
と $\mathrm{sinc}$ 関数に帰着します（ここで $\mathrm{sinc}(x) = \sin(x)/x$ の記法を用いました）。矩形窓のフーリエ変換が $\mathrm{sinc}$ になるという、信号処理の基本事実とぴたりと一致しています。直流（$\omega=0$）での利得は $RM$ で、これは窓長そのものです。
$N$ 段の場合は、振幅応答が単純に $N$ 乗されます。
$$ \left|H(e^{j\omega})\right| = \left|\frac{\sin\!\left(\dfrac{\omega RM}{2}\right)}{\sin\!\left(\dfrac{\omega}{2}\right)}\right|^{N} \;\approx\; \left(RM\right)^{N}\bigl|\mathrm{sinc}\bigl(\tfrac{\omega RM}{2}\bigr)\bigr|^{N} $$
これが「CICフィルタの周波数応答は $|\mathrm{sinc}|^N$ 型になる」という有名な性質です。直流利得は $(RM)^N$ と巨大になるため、実装では後で正規化（ビットシフト）が必要になります。
零点とエイリアス帯の減衰
振幅応答の分子 $\sin(\omega RM/2)$ がゼロになるのは
$$ \frac{\omega RM}{2} = \pi k \quad\Longrightarrow\quad \omega = \frac{2\pi k}{RM}, \quad k = 1, 2, \dots $$
のときです。すなわち $\omega = 2\pi/(RM)$ の整数倍に、深さ無限大（$N$ 段なら $N$ 重）の零点（ヌル）が等間隔に並びます。これがエイリアシング対策の要です。
デシメーション率を $R$ とすると、間引き後のサンプリングレートは $1/R$ になり、新しいナイキスト帯域は $\omega < \pi/R$（元の正規化周波数で）に縮みます。間引きによって、元の信号の $\omega = 2\pi k/R$ 付近の帯域（折り返しを起こす危険帯域、エイリアス帯）が、間引き後の通過帯域 $\omega \in [-\omega_c, \omega_c]$ へ折り返してきます。
ここで $M=1$ の標準ケースを考えると、CICのヌルは $\omega = 2\pi k/R$ にちょうど位置します。つまり 折り返してくる危険帯域の中心に、CICフィルタが深いヌルを配置している のです。これは偶然ではなく、くし形フィルタの差分遅延を間引き率に合わせて設計した必然的な結果です。エイリアス帯の幅（通過帯域端 $\omega_c$ で決まる）にわたる最小減衰量は、近似的に次式で見積もれます。
$$ \text{エイリアス帯の最悪減衰} \approx 20 N \log_{10}\!\left(\frac{\sin\!\left(\dfrac{\pi(1 – f_c/R\cdot R)}{?}\right)}{\cdots}\right)\ \text{[dB]} $$
——と書くと複雑なので、実用的にはもっとシンプルに考えます。エイリアスが最も浅く減衰される（最悪の）周波数は、最初のヌル $\omega = 2\pi/R$ の 手前側のエッジ、すなわち $\omega = 2\pi/R – \omega_c$ です。ここでの減衰量を
$$ A_{\min} \approx -20 N \log_{10}\!\left|\frac{\sin\!\left(\dfrac{\pi(1/R – f_c)\cdot R M}{?}\right)}{\sin(\cdots)}\right| $$
と細かく追うより、設計指針として押さえるべきは次の3点です。
- 段数 $N$ を増やすと、すべてのヌル・阻止帯域の減衰が $N$ 倍（dB単位）深くなる。 エイリアス抑圧を $A$ dB欲しいなら、$N$ を大きくすればよい。
- 通過帯域比 $f_c / f_{s,\text{out}}$ を小さく取る ほど、通過帯域がヌルから遠ざかり、エイリアス帯の減衰が深くなる。一般に通過帯域はナイキストの $1/4$ 程度以下に抑える設計が多い。
- 差分遅延 $M$ を $2$ にすると、ヌルの位置が $\omega = 2\pi k/(RM)$ と $M$ 倍密になる ため、より広いエイリアス帯をカバーできるが、その分通過帯域も狭くなる。
段数 $N$ を増やすことの代償が、次節で扱う「通過帯域ドループ」の悪化です。$\mathrm{sinc}^N$ の頂上付近は $N$ が大きいほど鋭く垂れ下がるため、通過帯域内の信号まで減衰させてしまうのです。
その前に、CICフィルタの真の威力——デシメーションと一体化したときの効率——を見ておきましょう。
デシメーションとの結合（Noble恒等式）
なぜ積分器を高速側に置けるのか
ここまでは「移動平均してから間引く」という順序を暗黙に想定してきました。しかし素朴に考えると、移動平均（積分器＋くし形）を高いレートで全部やってから間引くのでは、せっかくレートを下げても演算量が減りません。CICフィルタが効率的なのは、くし形フィルタを間引き後の低速側に移動できる からです。その理論的根拠が Noble恒等式（多レート恒等式） です。
Noble恒等式は、デシメータ $\downarrow R$ とフィルタの順序を入れ替える規則で、次のように述べられます。
$$ \bigl[\downarrow R\bigr] \to H(z) \quad\Longleftrightarrow\quad H(z^{R}) \to \bigl[\downarrow R\bigr] $$
言葉で言えば、「$R$ で間引いた後に $H(z)$ をかける」ことは、「先に $H(z^R)$ をかけてから $R$ で間引く」ことと等価、ということです。逆向きに使えば、$H(z^R)$ という形のフィルタは間引きの後ろ側（低速側）に移せます。
くし形フィルタの伝達関数は $H_C(z) = 1 – z^{-RM}$ でした。これを $z \to z^R$ の置き換えで見ると、$M=1$ のとき低速側のくし形 $1 – z^{-M}$ を高速側に戻すと $1 – z^{-RM}$ になります。逆に、高速側で $1-z^{-RM}$ として動かしていたくし形は、Noble恒等式によって間引き後の低速側へ「$1 – z^{-M}$」という短い遅延のくし形として移動できるのです。
一方、積分器 $1/(1-z^{-1})$ は $z^{-1}$ という1サンプルの遅延を含むため、間引きをまたいで移動できません（$z^{-R}$ の形ではないため）。よって積分器は必ず高速側に残ります。これがCICの名前の順序「Integrator（高速側）→ デシメーション → Comb（低速側）」を決めています。
結合後の構造と演算コスト
以上を組み合わせると、デシメーション型CICフィルタは次の構造になります。
- 高速レート $f_s$ で、積分器を $N$ 段（各段 $y_I[n] = y_I[n-1] + x[n]$）
- $\downarrow R$ で $R$ サンプルに1個だけ取り出す
- 低速レート $f_s/R$ で、差分遅延 $M$ のくし形を $N$ 段（各段 $y_C[n] = x[n] – x[n-M]$）
積分器は高速側にあるので $f_s$ レートで動きますが、各段の演算は「1回の加算」だけです。くし形は低速側なので $f_s/R$ レートでしか動かず、各段「1回の減算」だけ。合計の演算量は、サンプルあたり高速側で $N$ 加算、低速側で $N$ 減算（実効的に $N/R$ 加算相当）です。窓長 $RM$ がどれだけ大きくても、この演算量は変わりません。乗算は一切登場しません。これがFPGAで数百MHz級の信号をリアルタイム処理できる秘密です。
インターポレーション（アップサンプリング）の場合は順序が逆になり、「低速側でくし形 → $\uparrow R$ → 高速側で積分器」という構造になりますが、原理は同じです。本記事ではデシメーションに焦点を当てます。
効率の代償として、実装上もうひとつ向き合うべき問題があります。積分器は累積和を取るため、内部の値がどんどん大きくなり、オーバーフローの危険があります。次節で、このビット成長を定量的に見積もりましょう。
固定小数点でのビット成長
なぜビット幅が増えるのか
積分器 $y_I[n] = y_I[n-1] + x[n]$ は累積和なので、入力が同符号で続くと出力は際限なく大きくなります。直流利得が $(RM)^N$ にもなることを思い出すと、$N$ 段積分器を通った後の信号は入力の最大 $(RM)^N$ 倍にまで膨れ上がる可能性があります。固定小数点（整数演算）でこれをオーバーフローさせずに扱うには、内部レジスタのビット幅を十分に確保しなければなりません。
ここで重要なのは、CICフィルタの内部はオーバーフローしてよい という驚くべき事実です。積分器は2の補数のモジュロ演算（ラップアラウンド）で値が一周しても、後段のくし形が正確に差分を取れば、極零相殺が完全に成り立つ限り正しい出力が得られます。ただしこれは「全段を通したトータルのダイナミックレンジを収容できるビット幅をすべての段で確保する」ことが条件です。中途半端にビットを削ると壊れます。
Hogenauerのビット幅公式
Hogenauer のオリジナル論文では、各段で必要なビット幅が解析されていますが、最終出力で必要なビット幅 $B_{\text{out}}$ の見積もりとして、入力ビット幅を $B_{\text{in}}$ とすると、最大利得 $(RM)^N$ をビットで表現するのに必要な追加ビット数（ビット成長 $B_{\text{growth}}$）は
$$ B_{\text{growth}} = \left\lceil N \log_2 (RM) \right\rceil $$
で与えられます。したがって全段でオーバーフローを防ぐのに必要なビット幅は
$$ B_{\max} = B_{\text{in}} + \left\lceil N \log_2 (RM) \right\rceil $$
となります。この式の意味は明快です。直流利得 $(RM)^N$ を2進で表すのに $\log_2 (RM)^N = N\log_2(RM)$ ビット必要で、それを入力ビット幅に上乗せするだけです。
たとえば $N=4$、$R=64$、$M=1$、$B_{\text{in}}=16$ ビットなら
$$ B_{\text{growth}} = \lceil 4 \log_2 64\rceil = \lceil 4 \times 6 \rceil = 24\ \text{ビット} $$
となり、内部レジスタは $16 + 24 = 40$ ビット必要です。$R$ が大きいほど、$N$ が大きいほど、内部ビット幅が急速に増えることがわかります。実際の設計では、各積分器段・くし形段で必要なビット幅を個別に計算し、後段ほど削れる場合もある（Hogenauer のプルーニング）ことが知られていますが、安全側の設計としてはすべての段を $B_{\max}$ にするのが簡単です。
最後に正規化として、出力を $(RM)^N$ で割る（あるいは $B_{\text{growth}}$ ビット右シフトする）ことで利得を $1$ 付近に戻します。$RM$ が2のべきなら単純なビットシフトで済むため、ここでも乗算・除算は不要です。
ビット成長の問題はハードウェア固有の話なので、Pythonの浮動小数点シミュレーションでは表面化しません。とはいえ、固定小数点で実装する読者のために重要な設計パラメータなので押さえておきましょう。次は、CICの最大の弱点——通過帯域ドループ——と、その補償方法に進みます。
通過帯域ドループとFIRコンペンセータ
ドループはなぜ起きるか
理想的な低域通過フィルタは、通過帯域内で利得が完全に平坦（$=1$）であってほしいものです。ところがCICの振幅応答は $|\mathrm{sinc}(\omega RM/2)|^N$ という形をしているので、$\omega=0$ から離れるにつれて緩やかに減衰していきます。$\mathrm{sinc}$ 関数は頂上が丸く垂れ下がっているため、通過帯域の端（カットオフ付近）では利得が $1$ より目立って小さくなります。この通過帯域内での利得の垂れ下がりを 通過帯域ドループ（passband droop） と呼びます。
ドループの大きさを定量化しましょう。間引き後のナイキスト周波数の何分の1かを通過帯域端 $f_c$ とすると、正規化通過帯域端は $\omega_c = 2\pi f_c / f_s$ です。出力レート基準で $f = f_c$ における利得は（正規化後、直流利得 $(RM)^N$ で割って）
$$ \frac{|H(e^{j\omega_c})|}{(RM)^N} = \left|\frac{\sin(\omega_c RM/2)}{RM\,\sin(\omega_c/2)}\right|^{N} \approx \bigl|\mathrm{sinc}(\omega_c RM/2)\bigr|^{N} $$
たとえば通過帯域端を出力ナイキストの $1/2$（$\omega_c RM/2 = \pi/2$ 付近）まで取ると、$\mathrm{sinc}(\pi/2) = \sin(\pi/2)/(\pi/2) = 2/\pi \approx 0.637$ なので、$N=4$ 段では $0.637^4 \approx 0.165$、すなわち $-15.6$ dBもの減衰になってしまいます。これでは通過帯域の高い周波数成分がごっそり削られてしまい、使い物になりません。実用上は通過帯域を狭く取る（ヌルから遠ざける）か、後段で補償をかけるかのいずれかが必要です。
CICコンペンセータの設計
ドループ補償の王道は、CICの後段（間引き後の低速側）に、$\mathrm{sinc}^N$ 特性の逆特性を近似する小さなFIRフィルタ——CICコンペンセータ（inverse sinc filter） ——を縦続することです。低速側に置くので、レートが下がっておりFIR乗算のコストが小さくて済みます。
理想的な補償フィルタの振幅特性は、CICのドループを打ち消すものですから
$$ \left|H_{\text{comp}}(e^{j\omega})\right| = \frac{1}{\bigl|\mathrm{sinc}(\omega RM/2)\bigr|^{N}} = \left(\frac{\omega RM/2}{\sin(\omega RM/2)}\right)^{N} $$
を通過帯域内で満たせばよい、ということになります（ただし阻止帯域では暴れないように打ち切る必要があります）。これは通過帯域端に向かって利得が増えていく「高域ブースト」の特性で、$1/\mathrm{sinc}$ 形からインバースsincフィルタとも呼ばれます。
実装上は、間引き後の正規化周波数 $f’ = f/f_{s,\text{out}}$ を用いて、CICの低速側から見た振幅特性を考えます。低速側では1段あたりの応答が $\mathrm{sinc}(\pi M f’ R / R) = \mathrm{sinc}(\pi M f’)$ ……と、ここは丁寧に扱う必要があります。間引き後の周波数 $f’$ は出力レートで正規化されており、CICの $\mathrm{sinc}$ 引数は出力レート基準では $\sin(\pi M f’) / \sin(\pi M f’/R)$ という形になりますが、補償の実用設計では「通過帯域内で $1/\mathrm{sinc}^N$ を目標応答とするFIRフィルタを、周波数サンプリング法や firwin2 / remez で設計する」アプローチが最も手軽で確実です。
代表的な簡易補償として、3タップの対称FIR $[-\alpha,\ 1+2\alpha,\ -\alpha]$（$\alpha$ は $1/16$ や $1/8$ 程度）を使う方法も広く知られています。これは中心を持ち上げ両端をわずかに下げる山型特性で、ちょうど $\mathrm{sinc}$ の逆を1次近似します。$\alpha$ を $N$ と通過帯域幅に応じて調整すると、通過帯域のリプルを数分の1〜十分の1に圧縮できます。
補償をかけるとドループが平坦化される一方、阻止帯域では補償フィルタの利得が暴れないよう、補償FIRの遮断特性も設計しなければなりません。次のPython実装では、CICの応答そのものと、firwin2 による逆sinc補償フィルタを設計し、補償前後で通過帯域がどれだけ平坦になるかを実際に確かめます。
Pythonでの実装
積分器・くし形構造でのCICデシメータ
まずは、CICフィルタを「scipyの一般FIRに丸投げ」せず、積分器とくし形の構造そのままに実装します。これがハードウェアの動作と対応しています。
import numpy as np
def cic_decimator(x, R, M=1, N=4):
"""積分器・くし形構造によるCICデシメータ。
x : 入力信号（高速レート）
R : デシメーション率
M : 差分遅延（通常 1 か 2）
N : 段数
戻り値 : 間引き後の信号（低速レート）
"""
# --- 1) 積分器を N 段（高速レート側）---
# 各段で累積和を取る。np.cumsum がまさに積分器そのもの。
v = x.astype(np.float64)
for _ in range(N):
v = np.cumsum(v)
# --- 2) デシメーション（R サンプルに 1 個取り出す）---
v = v[::R]
# --- 3) くし形を N 段（低速レート側、差分遅延 M）---
for _ in range(N):
delayed = np.concatenate([np.zeros(M), v[:-M]])
v = v - delayed
# --- 4) 利得正規化（直流利得 (R*M)^N で割る）---
v = v / (R * M) ** N
return v
np.cumsum を $N$ 回適用するのが積分器段、v - delayed（$M$ サンプル遅延との差）が くし形段に対応しています。間引き v[::R] を積分器とくし形の間に挟んでいる点が、Noble恒等式で導いた「積分器は高速側、くし形は低速側」という構造をそのまま反映しています。最後に直流利得 $(RM)^N$ で割って正規化しています。
このコードで実際にステップ入力や正弦波を通してみると、移動平均としての平滑化が確認できます。次に、伝達関数から計算した理論応答と、この構造実装が一致するかを周波数領域で検証しましょう。
周波数応答の可視化（sinc特性と零点）
CICフィルタの伝達関数 $H(z) = \left((1-z^{-RM})/(1-z^{-1})\right)^N$ から、scipy.signal.freqz で周波数応答を描きます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
def cic_freq_response(R, M, N, worN=8192):
"""CIC（高速レート側から見た）伝達関数の係数を作り、周波数応答を返す。"""
# 1段分の分子 (1 - z^{-RM}) と分母 (1 - z^{-1}) を N 段畳み込む
b_single = np.zeros(R * M + 1)
b_single[0] = 1.0
b_single[-1] = -1.0 # 1 - z^{-RM}
a_single = np.array([1.0, -1.0]) # 1 - z^{-1}
b, a = np.array([1.0]), np.array([1.0])
for _ in range(N):
b = np.convolve(b, b_single)
a = np.convolve(a, a_single)
w, H = signal.freqz(b, a, worN=worN)
return w, H
plt.figure(figsize=(10, 6))
R, M = 8, 1
for N in [1, 2, 3, 5]:
w, H = cic_freq_response(R, M, N)
mag = np.abs(H) / (R * M) ** N # 直流利得で正規化
mag_db = 20 * np.log10(mag + 1e-12)
plt.plot(w / np.pi, mag_db, label=f"N={N}")
plt.axvline(2 / R, color="gray", ls="--", alpha=0.6, label="first null (2/R)")
plt.xlabel("Normalized frequency (×π rad/sample)")
plt.ylabel("Magnitude [dB]")
plt.title(f"CIC magnitude response (R={R}, M={M})")
plt.ylim(-100, 5)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig("cic_freq_response.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、CICフィルタの本質的な特徴が読み取れます。第一に、$\omega = 2\pi k/(RM)$（ここでは $2/R = 0.25\pi$ の整数倍）に深いヌルが等間隔に並んでおり、まさに $\mathrm{sinc}$ 型の応答であることが確認できます。第二に、段数 $N$ を増やすほど、阻止帯域全体の減衰がdB単位で $N$ 倍深くなっています（$N=1$ のサイドローブが $-13$ dB程度なのに対し、$N=5$ では $-65$ dB程度）。第三に、その代償として、通過帯域（$\omega \approx 0$ 付近）でのカーブの垂れ下がり（ドループ）も $N$ が大きいほど顕著になっていることがわかります。これがエイリアス抑圧とドループのトレードオフです。
次に、最初のヌルが本当に折り返し帯域の中心に来ているか、デシメーション動作で確かめます。
デシメーション動作の確認
複数の正弦波を含む信号にCICデシメータを適用し、間引き後にエイリアスがどう抑えられるかを見ます。
import numpy as np
import matplotlib.pyplot as plt
# cic_decimator は前掲の定義を使う
fs = 8000.0 # 入力サンプリングレート [Hz]
R, M, N = 8, 1, 4 # 出力レートは fs/R = 1000 Hz
n = np.arange(0, 8000)
t = n / fs
# 通過帯域内の 100 Hz（残したい）と、折り返しを起こす 900 Hz（消したい）
# 900 Hz は出力ナイキスト 500 Hz を超えており、間引くと 100 Hz に化ける危険
x = np.sin(2 * np.pi * 100 * t) + 0.8 * np.sin(2 * np.pi * 900 * t)
y = cic_decimator(x, R, M, N)
fs_out = fs / R
f_out = np.fft.rfftfreq(len(y), 1 / fs_out)
Y = np.abs(np.fft.rfft(y)) / len(y) * 2
plt.figure(figsize=(10, 5))
plt.plot(f_out, Y, "b-")
plt.axvline(100, color="g", ls="--", alpha=0.7, label="100 Hz (kept)")
plt.xlabel("Frequency [Hz] (output rate)")
plt.ylabel("Amplitude")
plt.title(f"CIC decimated spectrum (R={R}, N={N}, fs_out={fs_out:.0f} Hz)")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig("cic_decimation.png", dpi=150, bbox_inches="tight")
plt.show()
このスペクトルから、CICフィルタがアンチエイリアシングフィルタとして機能していることが明確に読み取れます。残したかった 100 Hz 成分は出力にしっかり残っている一方、900 Hz の成分はほとんど消えています。もしフィルタなしで単純に間引いていたら、900 Hz は出力レート 1000 Hz で折り返して $1000 – 900 = 100$ Hz の偽信号となり、本物の 100 Hz と区別できなくなっていたはずです。CICのヌルが 900 Hz 付近（$2\pi \times 900/8000 \approx 0.225\pi$、最初のヌル $0.25\pi$ の近傍）に深い減衰を与えたおかげで、この致命的なエイリアスが未然に防がれたわけです。
ここまでで、デシメーション本来の役割は果たせました。最後に、通過帯域ドループを補償フィルタでどこまで平らにできるかを検証します。
CICコンペンセータによるドループ補償
間引き後の低速側に、$1/\mathrm{sinc}^N$ を近似する逆sinc FIRフィルタを firwin2 で設計し、補償前後の通過帯域を比べます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
R, M, N = 8, 1, 4
fs_out = 1.0 # 出力レートを 1 に正規化（fs_out=1）
ntaps = 33 # 補償 FIR のタップ数（奇数→線形位相）
fc = 0.4 # 補償する通過帯域端（出力ナイキスト=0.5 の 0.8 倍）
# --- 目標応答: 通過帯域で 1/sinc^N、阻止帯域で 0 ---
freqs = np.linspace(0, 0.5, 512) # 出力レート基準の正規化周波数 [0, 0.5]
# 出力レートから見た CIC の sinc 引数（間引き後なので RM が効く形）
# 通過帯域内で 1/|sinc(pi*M*f)|^N を目標にする（f は出力正規化周波数）
x_arg = np.pi * M * freqs
sinc_resp = np.ones_like(freqs)
nz = x_arg > 1e-9
sinc_resp[nz] = (np.sin(x_arg[nz]) / x_arg[nz])
target = np.ones_like(freqs)
pb = freqs <= fc
target[pb] = 1.0 / (sinc_resp[pb] ** N) # 通過帯域: 逆 sinc ブースト
target[~pb] = 0.0 # 阻止帯域: 遮断
# firwin2 用に周波数(0..1)と利得を渡す（1 = ナイキスト=0.5*fs_out*2）
f_norm = freqs / 0.5
h_comp = signal.firwin2(ntaps, f_norm, target)
# --- 補償フィルタと、CIC×補償の合成応答を評価 ---
w, Hc = signal.freqz(h_comp, worN=4096)
fw = w / np.pi * 0.5 # 出力レート基準 [0, 0.5]
# CIC 自身の（出力レート基準の）正規化応答
xc = np.pi * M * fw
cic = np.ones_like(fw)
nz2 = xc > 1e-9
cic[nz2] = np.abs(np.sin(xc[nz2]) / xc[nz2]) ** N
plt.figure(figsize=(10, 6))
plt.plot(fw, 20*np.log10(cic + 1e-12), label="CIC only (droop)")
plt.plot(fw, 20*np.log10(np.abs(Hc) + 1e-12), label="Compensator (1/sinc^N)")
plt.plot(fw, 20*np.log10(cic*np.abs(Hc) + 1e-12), 'k-', lw=2,
label="CIC × Compensator")
plt.axvline(fc, color="gray", ls="--", alpha=0.6, label=f"passband edge fc={fc}")
plt.xlabel("Normalized frequency (output rate, 0.5 = Nyquist)")
plt.ylabel("Magnitude [dB]")
plt.title(f"CIC droop compensation (R={R}, M={M}, N={N})")
plt.ylim(-30, 10)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig("cic_compensation.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフが補償の効果を端的に示しています。青のCIC単体の応答は通過帯域端 $f_c = 0.4$ に向かって明らかに垂れ下がっており、これがドループです。オレンジの補償フィルタは逆に通過帯域端へ向かって利得が持ち上がる「高域ブースト」特性を持っています。両者を掛け合わせた黒の合成応答は、通過帯域内でほぼ平坦（$0$ dB付近）になっており、ドループがきれいに打ち消されていることが読み取れます。補償フィルタは間引き後の低速レートで動くので、追加コストは小さく抑えられます。
最後に、補償前後の通過帯域ドループ量を数値で比較してみましょう。
import numpy as np
R, M, N = 8, 1, 4
fc = 0.4 # 通過帯域端（出力ナイキスト 0.5 基準）
# CIC 単体の通過帯域端での減衰
xc = np.pi * M * fc
cic_edge = np.abs(np.sin(xc) / xc) ** N
print(f"CIC only : 通過帯域端 fc={fc} での利得 = "
f"{20*np.log10(cic_edge):.2f} dB")
# 理想補償をかけた場合（1/sinc^N で完全相殺 → 0 dB）
print("Compensated: 通過帯域端での利得 ≈ 0.00 dB（理想補償）")
print(f"改善量 : {abs(20*np.log10(cic_edge)):.2f} dB のドループを補償")
この出力から、$R=8,\ M=1,\ N=4$ のCICでは通過帯域端で約 $-9.7$ dBものドループが生じており、補償フィルタによってこれをほぼ $0$ dBまで平坦化できることが定量的に確認できます。段数 $N=4$ と通過帯域端 $f_c=0.4$ がともに大きいため、$|\mathrm{sinc}|^N$ の垂れ下がりがこれほど深くなる点に注意してください。実際のSDR受信機では、このようにCICで大まかに間引いてから、低速側でコンペンセータと急峻なハーフバンドFIRを組み合わせ、最終的なチャネル選択を行う多段構成が標準的に用いられます。
まとめ
本記事では、乗算器を使わず加算と遅延だけで高効率なデシメーションを実現するCICフィルタについて、構造・伝達関数・周波数特性・実装・補償まで通して解説しました。
- CICの正体は乗算器なしの移動平均: 積分器 $1/(1-z^{-1})$ とくし形 $1-z^{-RM}$ を縦続すると、等比級数和の公式から伝達関数が長さ $RM$ の矩形窓 $\sum_{k=0}^{RM-1}z^{-k}$ になる。積分器の極 $z=1$ とくし形の零点 $z=1$ が相殺し、安定なFIRに化ける
- 周波数応答は $|\mathrm{sinc}|^N$ 型: 半角公式で整理すると $|H(e^{j\omega})| = |\sin(\omega RM/2)/\sin(\omega/2)|^N$ となり、$\omega = 2\pi k/(RM)$ に等間隔の深いヌルを持つ。段数 $N$ を増やすと阻止帯域減衰がdB単位で $N$ 倍深くなる
- エイリアス帯にヌルを配置: 差分遅延 $M$ を間引き率に合わせると、折り返してくる危険帯域の中心にCICのヌルが来て、エイリアスを抑圧する
- Noble恒等式で積分器を高速側・くし形を低速側へ: 窓長 $RM$ によらずサンプルあたり一定の加減算だけで動き、乗算ゼロを達成
- ビット成長: 内部レジスタは $B_{\text{in}} + \lceil N\log_2(RM)\rceil$ ビット必要。CIC内部はオーバーフローしてよいが、全段でダイナミックレンジを確保すること
- 通過帯域ドループとCICコンペンセータ: $\mathrm{sinc}^N$ の垂れ下がりを、低速側に置いた $1/\mathrm{sinc}^N$ 近似の逆sinc FIRで補償し、通過帯域を平坦化する
CICフィルタは、エイリアス抑圧（$N$ を増やす）と通過帯域平坦性（補償フィルタ）のバランスを取りながら、超高速ストリームを実用的なレートまで落とす「第一段」として理想的です。実際の受信機では、CICで粗く間引いた後、低速側でコンペンセータと急峻なFIRを段階的にかける多段デシメーションが定石になっています。
次のステップとして、以下の記事も参考にしてください。