ソフトウェア無線（SDR）の受信機は、アンテナから入ってきた広帯域の信号を、いったん高いサンプリングレートで取り込みます。しかし最終的に欲しいのは、その中のごく狭い1チャネル分だけだったりします。あるいは、5G基地局のように1本の広帯域を128個のサブチャネルに一気に分けて、各ユーザーへ割り当てたいこともあります。このとき素直に「帯域制限フィルタを通してから間引く」「128個のフィルタを別々に走らせる」と実装すると、毎秒数十億回の積和演算が必要になり、消費電力もチップ面積も現実離れしてしまいます。
この計算量の壁を打ち破る切り札がポリフェーズフィルタバンクです。鍵となるアイデアは2つあります。1つは、間引いて捨てる予定のサンプルをそもそも計算しないという「Noble恒等式（高貴な恒等式）」による順序の入れ替え。もう1つは、たくさんのバンドパスフィルタが共通して持つ「同じプロトタイプを周波数シフトしただけ」という構造を、IDFT（逆離散フーリエ変換）でまとめて計算してしまうという発想です。この2つを組み合わせると、$M$ チャネルのフィルタバンクの計算量が、チャネル数 $M$ に比例する素朴な実装から、$\log M$ にしか依存しない劇的に軽い実装へと変わります。
ポリフェーズ分解は、SDRのチャネライザ、5G/LTEのOFDM変復調、レーダーのパルス圧縮前処理、デジタル音声放送（DAB）のサブバンド分割など、現代の通信・信号処理ハードウェアの至るところで使われています。本記事では、なぜこの効率化が可能なのかを、z変換とNoble恒等式から1行ずつ導出し、最後にPythonで「素朴な実装と完全に同じ出力を、はるかに少ない演算で得られる」ことと「各チャネルの周波数応答」を可視化して確認します。
本記事の内容
- マルチレート信号処理（間引き・補間）の復習と、なぜ無駄が生じるか
- Noble恒等式（間引き・補間とフィルタの交換則）の導出
- プロトタイプフィルタのポリフェーズ分解（z変換による導出、省略なし）
- ポリフェーズ間引き器による演算量削減の仕組み
- 一様DFTフィルタバンクへの帰着（IDFT/FFTとの結合）
- 演算量が $O(MN)$ から $O(N + M\log M)$ に下がる理由
- Pythonによるポリフェーズ間引き器とDFTフィルタバンクの実装・可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、間引きと補間の z 領域での表現、FIRフィルタの畳み込み、そしてDFTの定義は本記事の土台になります。z変換と単位円上の周波数応答については、ディジタルフィルタの基礎（FIR/IIR）を解説 も参照してください。
マルチレート処理の復習：間引きと補間
ポリフェーズ分解の話に入る前に、舞台となるマルチレート処理の2つの基本操作を、記号と z 領域での表現まで揃えておきましょう。後の導出はすべてこの記法の上で進みます。
ダウンサンプラ（間引き）
まず直感から始めます。サンプリングレートを下げるとは、たとえば「10サンプルに1つだけ残して、残り9個を捨てる」という操作です。連続的な波を細かく刻んだ点列のうち、間引いて間隔を空けるイメージです。これをダウンサンプラ（factor $M$ のダウンサンプリング、$\downarrow M$）と呼びます。入力 $x[n]$ に対する出力 $y_D[n]$ は
$$ y_D[n] = x[Mn] $$
です。$M$ 個に1個のサンプルだけを取り出し、出力のインデックス $n$ を詰め直します。出力のレートは入力の $1/M$ になります。
このダウンサンプリングの z 変換を求めておきましょう。これは後で何度も使う最重要の道具です。出力の z 変換は定義より
$$ Y_D(z) = \sum_{n=-\infty}^{\infty} x[Mn]\, z^{-n} $$
です。ここで「$M$ の倍数番目だけを拾う」という操作を数式で表現するために、くし型（インパルス列）和の恒等式を使います。任意の整数 $k$ について
$$ \frac{1}{M}\sum_{l=0}^{M-1} W_M^{lk} = \begin{cases} 1 & (k \equiv 0 \bmod M) \\ 0 & (\text{それ以外}) \end{cases} \qquad \left(W_M \equiv e^{-j2\pi/M}\right) $$
が成り立ちます。$W_M$ は1の原始 $M$ 乗根で、$k$ が $M$ の倍数のときだけ全項が1になって和が $M$ になり、それ以外は対称性で打ち消し合ってゼロになる、というよく知られた直交関係です。この「フィルタ」を $x[k]$ に掛けて全インデックスで和を取ると、$M$ の倍数番目だけが生き残ります。これを使うと
$$ \begin{align} Y_D(z) &= \sum_{n=-\infty}^{\infty} x[Mn]\, z^{-n} \\ &= \sum_{k=-\infty}^{\infty} x[k]\left(\frac{1}{M}\sum_{l=0}^{M-1} W_M^{lk}\right) z^{-k/M} \end{align} $$
と書き直せます。2行目では、和を全整数 $k$ にわたって取り直し、$k$ が $M$ の倍数のときだけくし型和が1になることを利用しました（$z^{-k/M}$ は $k=Mn$ のとき $z^{-n}$ に一致します）。和の順序を入れ替えて $l$ の和を外に出すと
$$ \begin{align} Y_D(z) &= \frac{1}{M}\sum_{l=0}^{M-1}\sum_{k=-\infty}^{\infty} x[k]\, W_M^{lk}\, z^{-k/M} \\ &= \frac{1}{M}\sum_{l=0}^{M-1}\sum_{k=-\infty}^{\infty} x[k]\,\left(W_M^{-l} z^{1/M}\right)^{-k} \\ &= \frac{1}{M}\sum_{l=0}^{M-1} X\!\left(W_M^{-l} z^{1/M}\right) \end{align} $$
を得ます。ここで $W_M^{lk} z^{-k/M} = (W_M^{-l}z^{1/M})^{-k}$ とまとめ、内側の和が定義により $X(\,\cdot\,)$ になることを使いました。最終形
$$ Y_D(z) = \frac{1}{M}\sum_{l=0}^{M-1} X\!\left(z^{1/M} W_M^{-l}\right) $$
が、ダウンサンプリングの z 領域表現です。$X(z)$ を $z^{1/M}$ で「伸ばし」、さらに $M$ 通りに位相回転させたものを足し合わせる――これがスペクトルの折り返し（エイリアシング）の正体です。$z=e^{j\omega}$ とおけば $\omega$ 軸方向に $M$ 倍に拡大した複製が $2\pi/M$ ごとに重なり合う、というおなじみの図になります。
アップサンプラ（補間）
逆に、レートを上げる操作も用意しておきます。アップサンプラ（factor $L$、$\uparrow L$）は、サンプルとサンプルの間に $L-1$ 個のゼロを挿入する操作です。
$$ y_U[n] = \begin{cases} x[n/L] & (n \text{ が } L \text{ の倍数}) \\ 0 & (\text{それ以外}) \end{cases} $$
この z 変換は、ゼロ挿入なので欠落項が消えるだけで素直に求まります。
$$ \begin{align} Y_U(z) &= \sum_{n=-\infty}^{\infty} y_U[n]\, z^{-n} = \sum_{m=-\infty}^{\infty} x[m]\, z^{-Lm} = X(z^L) \end{align} $$
2番目の等号では、ゼロでない項だけ（$n=Lm$）を拾って $m$ で書き直しました。$Y_U(z)=X(z^L)$ は、スペクトルが $\omega$ 軸方向に $1/L$ に縮み、$2\pi/L$ ごとにイメージ（複製）が現れることを意味します。
ここまでで、間引き・補間それぞれの z 領域での顔が出揃いました。次に問題になるのは、これらの操作とフィルタリングをどの順番で行うかです。じつはこの順番の自由度こそが、ポリフェーズによる効率化の源泉なのです。
なぜ素朴な間引きは無駄なのか
間引きには鉄則があります。レートを $1/M$ に下げる前には、必ず帯域を $\pi/M$ 以下に制限する低域通過フィルタ（アンチエイリアシングフィルタ）を通さなければなりません。さもないと、先ほど見た折り返しで高域成分が低域に化けて、信号が壊れてしまうからです。したがって標準的な間引き器（デシメータ）は
$$ x[n] \;\longrightarrow\; \boxed{H(z)} \;\longrightarrow\; \boxed{\downarrow M} \;\longrightarrow\; y[n] $$
という構成になります。$H(z)$ は $N$ タップのFIRローパス（プロトタイプフィルタ）です。
ここに巨大な無駄があります。フィルタ $H(z)$ は入力レート、つまり高いレートで全サンプルに対して畳み込みを計算します。$N$ タップなら入力1サンプルあたり $N$ 回の積和です。ところがその直後のダウンサンプラは、計算した出力の $M$ 個に1個しか使わず、残り $M-1$ 個を捨てます。
つまり「これから捨てるサンプルを、わざわざ全部フィルタ計算している」のです。入力1サンプルあたり $N$ 回の積和を払っているのに、有効な出力は $1/M$ しかない。捨てる分まで含めれば $N$ 回中 $(M-1)/M$ 倍の計算が文字どおりゴミ箱行きです。
直感的には、「フィルタの後で間引く」のではなく「先に間引いてからフィルタすれば、捨てる分は最初から計算しなくて済むのでは？」と思いたくなります。しかしフィルタは過去のサンプルとの畳み込みなので、勝手に順番を入れ替えると結果が変わってしまいます。この「順番を入れ替えても結果が変わらない」ための厳密な条件を与えるのが、次に説明するNoble恒等式です。これがポリフェーズ分解の論理的な心臓部になります。
Noble恒等式：間引き・補間とフィルタの交換
Noble恒等式（恒等式の名は「気高い・基本的」の意）は、ダウンサンプラ／アップサンプラとフィルタの順序を入れ替える規則を与えます。ここでは間引き側だけ厳密に導出し、補間側は結果を述べます。
主張
Noble恒等式（間引き）: ダウンサンプラ $\downarrow M$ の後に伝達関数 $G(z)$ のフィルタを置いたものは、ダウンサンプラの前に $G(z^M)$ を置いたものと等価である。図式で書くと
$$ \boxed{\downarrow M} \to \boxed{G(z)} \quad\Longleftrightarrow\quad \boxed{G(z^M)} \to \boxed{\downarrow M} $$
「低レート側でかける $G(z)$」と「高レート側でかける $G(z^M)$」は同じ出力を生む、という主張です。ここで $G(z^M)$ は、$G(z)$ のインパルス応答の各係数の間に $M-1$ 個のゼロを挟んで引き伸ばしたフィルタです。
導出
入力を $x[n]$、間引き後を $v[n]=x[Mn]$、それに $G(z)$ をかけた最終出力を $y[n]$ とします。右辺（先に $G(z^M)$）と左辺（後で $G(z)$）の出力が一致することを z 領域で示します。
まず左辺（先に間引き、後で $G(z)$）。間引きの z 表現を使うと、$V(z) = \frac{1}{M}\sum_{l} X(z^{1/M}W_M^{-l})$ で、出力は
$$ Y_{\text{left}}(z) = G(z)\, V(z) = G(z)\cdot \frac{1}{M}\sum_{l=0}^{M-1} X\!\left(z^{1/M}W_M^{-l}\right) $$
次に右辺（先に $G(z^M)$、後で間引き）。$G(z^M)$ をかけた中間信号を $w[n]$ とすると $W(z)=G(z^M)X(z)$ です。これを $\downarrow M$ に通すと、間引きの z 表現の $X$ を $W$ に置き換えて
$$ Y_{\text{right}}(z) = \frac{1}{M}\sum_{l=0}^{M-1} W\!\left(z^{1/M}W_M^{-l}\right) = \frac{1}{M}\sum_{l=0}^{M-1} G\!\left((z^{1/M}W_M^{-l})^M\right) X\!\left(z^{1/M}W_M^{-l}\right) $$
ここで核心の一手です。$G$ の引数を展開すると
$$ \left(z^{1/M}W_M^{-l}\right)^{M} = z\, W_M^{-lM} = z\cdot \left(e^{-j2\pi/M}\right)^{-lM} = z\cdot e^{j2\pi l} = z $$
となります。$W_M^{-lM}=e^{j2\pi l}=1$（$l$ は整数）を使いました。つまり $M$ 乗すると位相回転 $W_M^{-l}$ がきれいに消えて、$G$ の引数は $l$ によらず $z$ になります。すると $G$ を和の外にくくり出せて
$$ Y_{\text{right}}(z) = \frac{1}{M}\sum_{l=0}^{M-1} G(z)\, X\!\left(z^{1/M}W_M^{-l}\right) = G(z)\cdot\frac{1}{M}\sum_{l=0}^{M-1} X\!\left(z^{1/M}W_M^{-l}\right) $$
これは先ほどの $Y_{\text{left}}(z)$ と完全に一致します。よって
$$ Y_{\text{left}}(z) = Y_{\text{right}}(z) $$
が示され、Noble恒等式（間引き）が成立します。$\square$
補間側についても、まったく同じ計算（$Y_U(z)=X(z^L)$ を使う）で
$$ \boxed{G(z)} \to \boxed{\uparrow L} \quad\Longleftrightarrow\quad \boxed{\uparrow L} \to \boxed{G(z^L)} $$
が示せます。アップサンプラの前で $G(z)$ をかけるのと、後で $G(z^L)$ をかけるのが等価、というわけです。
直感的に言い直すと、Noble恒等式は「フィルタの係数間隔を間引き率に合わせて引き伸ばす（$z\to z^M$）なら、フィルタを高レート側と低レート側のどちらに置いても結果は同じ」という保証です。この保証があるからこそ、私たちは安心して「フィルタを低レート側に移して計算を減らす」ことができます。問題は、プロトタイプフィルタ $H(z)$ をどうやって $G(z^M)$ の形に書き換えるか。その答えがポリフェーズ分解です。
プロトタイプフィルタのポリフェーズ分解
アイデア：係数を $M$ 個のグループに振り分ける
長さ $N$ のFIRプロトタイプ $H(z)=\sum_{n=0}^{N-1}h[n]z^{-n}$ を考えます。ポリフェーズ分解の発想は単純で、「係数 $h[n]$ を、インデックス $n$ を $M$ で割った余り（位相）ごとに $M$ 個のグループに仕分ける」ことです。$M=3$、$h=(h_0,h_1,h_2,h_3,h_4,h_5,\dots)$ なら
- 位相0: $h_0, h_3, h_6, \dots$
- 位相1: $h_1, h_4, h_7, \dots$
- 位相2: $h_2, h_5, h_8, \dots$
と分けます。「ポリフェーズ（多相）」の名は、この $M$ 個の位相成分に由来します。
z 領域での厳密な導出
仕分けを数式で表します。任意のインデックス $n$ は $n = Mm + p$（$p=0,1,\dots,M-1$ は位相、$m$ はグループ内の番号）と一意に書けるので、和を $p$ と $m$ の二重和に分解します。
$$ \begin{align} H(z) &= \sum_{n=0}^{N-1} h[n]\, z^{-n} \\ &= \sum_{p=0}^{M-1}\sum_{m} h[Mm+p]\, z^{-(Mm+p)} \end{align} $$
2行目で $n=Mm+p$ を代入し、$n$ の単一和を「位相 $p$ の和」×「グループ内 $m$ の和」に組み替えました。指数を $z^{-(Mm+p)}=z^{-p}\,(z^{M})^{-m}$ と分けて、$p$ に依存する $z^{-p}$ を $m$ の和の外に出すと
$$ H(z) = \sum_{p=0}^{M-1} z^{-p} \underbrace{\sum_{m} h[Mm+p]\,(z^{M})^{-m}}_{\displaystyle E_p(z^M)} $$
ここでポリフェーズ成分 $E_p(z)$ を
$$ E_p(z) \equiv \sum_{m} h[Mm+p]\, z^{-m} \qquad (p=0,1,\dots,M-1) $$
と定義します。$E_p(z)$ は、$p$ 番目の位相の係数 $\{h[p], h[M+p], h[2M+p],\dots\}$ だけからなる、長さおよそ $N/M$ の短いFIRフィルタです。これを使うと、プロトタイプは
$$ \boxed{\,H(z) = \sum_{p=0}^{M-1} z^{-p}\, E_p(z^M)\,} $$
ときれいに分解できます。これがポリフェーズ恒等式（タイプ1）です。$M$ 本の短いフィルタ $E_p$ に、遅延 $z^{-p}$ を付けて足し合わせるだけで、元の長いフィルタ $H(z)$ がそっくり再現される、という事実です。係数を1つも捨てておらず、単なる並べ替えなので、これは厳密な恒等式です。
間引き器への適用
いよいよNoble恒等式と組み合わせます。間引き器 $H(z)\to\downarrow M$ にポリフェーズ分解を代入すると、フィルタ部は
$$ H(z) = \sum_{p=0}^{M-1} z^{-p}\, E_p(z^M) $$
なので、間引き器の出力 $Y(z)$ は、各枝を $\downarrow M$ に通したものの和になります。$p$ 番目の枝は「$z^{-p}$ で遅延 → $E_p(z^M)$ でフィルタ → $\downarrow M$ で間引き」という流れです。ここで $E_p(z^M)\to\downarrow M$ の部分にNoble恒等式を適用すると、$E_p(z^M)$ を間引きの後ろへ移して $E_p(z)$ に変えられます。
$$ \boxed{E_p(z^M)} \to \boxed{\downarrow M} \quad\Longleftrightarrow\quad \boxed{\downarrow M} \to \boxed{E_p(z)} $$
つまり、各枝は「$z^{-p}$ で遅延 → $\downarrow M$ で間引き → $E_p(z)$ でフィルタ」に置き換わります。すると、フィルタ $E_p$ は間引いた後の低レート側で動くことになります。これが効率化の決定打です。
- フィルタ $E_p$ は長さ $N/M$ に縮んでいる
- しかも入力の $1/M$ のレートでしか動かない
- 捨てるサンプルは間引きで先に捨てられ、フィルタは計算しない
直感的な実装イメージはこうです。入力ストリームを「整流子（コミュテータ）」で $M$ 本のレーンに順に振り分けます（$x[0]$ をレーン0、$x[-1]$ をレーン1…という具合に $z^{-p}$ の遅延に対応）。各レーンには短いフィルタ $E_p$ が1つずつ。$M$ サンプルたまるごとに、各レーンが1回ずつフィルタ出力を計算し、それらを足し合わせると間引き出力が1つ出てくる――そういう構造です。
次は、この構造がどれだけ計算を減らすのかを定量的に評価し、さらにこの「$M$ 本の短いフィルタ」を1種類のプロトタイプにできる場合、IDFTでまとめて $M$ 本のフィルタを同時に得られることを示します。これがDFTフィルタバンクです。
演算量の評価：ポリフェーズ間引き器
素朴な間引き器とポリフェーズ間引き器の計算量を、出力1サンプルあたりの積和回数で比べます。
素朴な実装: $N$ タップのフィルタを入力レートで動かし、出力で間引きます。入力1サンプルあたり $N$ 回の積和。出力は入力の $1/M$ なので、出力1サンプルあたりに換算すると入力 $M$ サンプル分、すなわち $MN$ 回の積和です。
ポリフェーズ実装: $M$ 本のフィルタ $E_p$ がそれぞれ長さ $N/M$。1本あたり $N/M$ 回の積和で、$M$ 本合わせて $M\times (N/M) = N$ 回。これで出力が1サンプル出ます。
したがって、出力1サンプルあたりの積和回数は
$$ \text{素朴: } MN \quad\longrightarrow\quad \text{ポリフェーズ: } N $$
ちょうど $M$ 分の1です。間引き率 $M$ が大きいほど、ポリフェーズの恩恵は大きくなります。直感どおり、「捨てるサンプルを計算しない」分だけ得をしているわけです。
ここまでは「1チャネルを間引く」話でした。しかし冒頭で触れたように、現実には「$M$ チャネルへ同時に分割したい」という要求が頻繁にあります。次のセクションでは、ポリフェーズ構造にIDFT/FFTを結合すると、$M$ チャネルを一気に取り出せる一様DFTフィルタバンクが立ち現れることを導出します。
一様DFTフィルタバンクへの帰着
設定：プロトタイプを周波数シフトしてバンドを並べる
$M$ チャネルのフィルタバンクとは、1本の入力スペクトルを $M$ 個の等幅の隣り合うサブバンドに分け、各バンドを間引いて低レートで取り出す装置です。直感的には、ローパスのプロトタイプ $H(z)$ を周波数軸上で $2\pi k/M$ ずつずらして $M$ 個のバンドパスを並べ、それぞれを間引きます。$k$ 番目のチャネルのフィルタは、$h[n]$ を変調（周波数シフト）した
$$ h_k[n] = h[n]\, e^{j2\pi kn/M} = h[n]\, W_M^{-kn} $$
です。z 領域では、$e^{j2\pi k/M}=W_M^{-k}$ による変調はスペクトルのシフトに対応し
$$ H_k(z) = H\!\left(z\, W_M^{k}\right) $$
となります（変調定理：時間領域の $W_M^{-kn}$ 倍は z 領域で引数を $zW_M^{k}$ に置換）。$k=0,1,\dots,M-1$ の $M$ 本を素朴に別々の畳み込みで計算すると、当然 $M$ 倍の計算量です。ここにポリフェーズの魔法をかけます。
ポリフェーズ分解を代入する
各チャネルフィルタ $H_k(z)=H(zW_M^k)$ に、先ほどのポリフェーズ恒等式 $H(z)=\sum_{p}z^{-p}E_p(z^M)$ を代入します。引数 $z$ を $zW_M^k$ に置き換えると
$$ H_k(z) = H(zW_M^k) = \sum_{p=0}^{M-1} (zW_M^k)^{-p}\, E_p\!\left((zW_M^k)^M\right) $$
ここで2つの簡約を行います。第1に、$(zW_M^k)^{-p} = z^{-p}W_M^{-kp}$ と分けます。第2に、ポリフェーズ成分の引数 $(zW_M^k)^M = z^M W_M^{kM} = z^M$ です。Noble恒等式の導出で見たのと同じく、$W_M^{kM}=(e^{-j2\pi/M})^{kM}=e^{-j2\pi k}=1$ となって位相回転が消えます。したがって
$$ H_k(z) = \sum_{p=0}^{M-1} z^{-p}\, W_M^{-kp}\, E_p(z^M) $$
を得ます。ここが決定的に重要な点です。$M$ 本のチャネルフィルタ $H_k$（$k=0,\dots,M-1$）は、同じポリフェーズ成分 $E_p(z^M)$ と同じ遅延 $z^{-p}$ を共有していて、チャネルごとに違うのは係数 $W_M^{-kp}$ だけなのです。
IDFTが現れる
各 $p$ について、遅延とフィルタを通した共通の中間信号を
$$ U_p(z) \equiv z^{-p}\, E_p(z^M)\, X(z) $$
と定義すると、$k$ 番目のチャネル出力（間引き前）は
$$ Y_k(z) = H_k(z)X(z) = \sum_{p=0}^{M-1} W_M^{-kp}\, U_p(z) $$
この和の形を見てください。$\{U_p\}_{p=0}^{M-1}$ という $M$ 個の信号に対して、$k$ 番目の出力が $\sum_p W_M^{-kp}U_p$ ――これはまさに逆離散フーリエ変換（IDFT）の定義式そのものです。$M$ 点IDFTは
$$ \mathrm{IDFT}\{U_p\}[k] = \sum_{p=0}^{M-1} U_p\, W_M^{-kp} = \sum_{p=0}^{M-1} U_p\, e^{j2\pi kp/M} $$
（正規化定数の置き方は流儀によります）。つまり、$M$ 本のポリフェーズ枝の出力 $\{U_p\}$ を $M$ 点IDFTにかけるだけで、$M$ チャネル全部の出力 $\{Y_k\}$ が同時に得られます。$M$ 個の別々のバンドパスフィルタを走らせる必要はなく、共通のポリフェーズフィルタ群 + 1回のIDFTで済むのです。
実際のフィルタバンクでは、Noble恒等式で $E_p$ を低レート側へ移し、整流子で入力を $M$ レーンに振り分けて各レーンの短いフィルタ $E_p$ を低レートで動かし、その出力ベクトルを $M$ 点IDFT（実装上はIFFT）にかける、という流れになります。これが一様DFTフィルタバンクの標準構成です。
演算量：$O(MN)$ から $O(N + M\log M)$ へ
$M$ チャネルを出力1セット（各チャネル1サンプルずつ、計 $M$ 個）取り出すのに必要な積和回数を数えます。
素朴な実装: 各チャネルが $N$ タップのバンドパス。$M$ チャネルで $MN$ 回。出力1セットあたり $O(MN)$。
DFTフィルタバンク: ポリフェーズ部は $M$ 本の長さ $N/M$ フィルタで合計 $N$ 回の積和。IDFTは素朴なら $M^2$ ですが、$M$ を2のべきなどに選べばIFFTで $O(M\log M)$。合わせて
$$ O(MN) \;\longrightarrow\; O\!\left(N + M\log M\right) $$
となります。$N$ はプロトタイプのタップ数で、典型的には $N = K M$（チャネルあたり $K$ タップ、$K$ は数〜十数）程度。すると $N = O(M)$ オーダーなので、出力1セットの計算量は素朴の $O(M^2 K)$ 相当から $O(M + M\log M)=O(M\log M)$ へと落ちます。チャネル数 $M$ を128や1024に増やすほど、この差は天文学的になります。
なぜこれほど減るのかを一言でまとめると、(1) Noble恒等式とポリフェーズ分解で「捨てるサンプルを計算しない」ことで素朴な間引きの $M$ 倍の無駄を消し、(2) 全チャネルが共有する変調構造 $W_M^{-kp}$ をIDFT/FFTのバタフライ演算へ畳み込むことで、チャネル間の重複計算を $\log M$ オーダーまで圧縮した――この2段構えの効果です。
理屈は出揃いました。次は、これらが本当に「素朴な実装と同じ出力を、少ない演算で」実現するのかを、Pythonで実装して数値的に確かめます。
Pythonでの実装
1. ポリフェーズ間引き器：素朴な実装と一致するか
まず、プロトタイプFIRローパスを設計し、(a) 素朴な「フィルタ→間引き」と (b) ポリフェーズ間引き器が、ビット単位ではないにせよ数値誤差の範囲で一致することを確認します。
import numpy as np
from scipy import signal
# ===== プロトタイプ低域通過フィルタの設計 =====
M = 8 # 間引き率 = チャネル数
K = 12 # チャネルあたりのタップ数
N = K * M # プロトタイプの総タップ数（M の倍数にしておく）
# カットオフを 1/M（ナイキスト基準で 1/M）にした FIR ローパス
h = signal.firwin(N, 1.0 / M, window=('kaiser', 8.0))
print(f"プロトタイプタップ数 N = {len(h)}, 間引き率 M = {M}")
# 入力信号: 低域トーン + 高域トーン + ノイズ（高レート）
rng = np.random.default_rng(0)
n_in = 4000
nidx = np.arange(n_in)
x = (np.cos(2*np.pi*0.03*nidx) # 通過域に入る低い正弦波
+ 0.7*np.cos(2*np.pi*0.30*nidx) # 阻止域の高い正弦波（折り返す）
+ 0.2*rng.standard_normal(n_in)) # 白色雑音
ここではカットオフ $1/M$ のローパスを設計し、通過域に入る低周波トーンと、間引きで折り返すはずの高周波トーンを混ぜた入力を用意しました。プロトタイプ長 $N=KM$ を $M$ の倍数にしてあるので、ポリフェーズ成分はきれいに各 $K$ タップずつに分かれます。
# ===== (a) 素朴な実装: フィルタ → 間引き =====
y_full = signal.lfilter(h, 1.0, x) # 高レートで全サンプル計算
y_naive = y_full[::M] # M個に1個だけ残す（残りは捨てる）
# ===== (b) ポリフェーズ間引き器 =====
def polyphase_decimate(x, h, M):
"""係数を M 位相に分け、間引き後の低レートで各成分を畳み込んで合算する。"""
N = len(h)
# 位相 p の係数 h[p], h[M+p], ... を行に並べる（タイプ1分解）
E = h.reshape(-1, M).T # 形状 (M, N//M): E[p] が p 番目のポリフェーズ成分
# 入力を M レーンに整流子で振り分け（位相 p は p サンプル先行）
L = len(x) // M
Y = np.zeros(L)
for p in range(M):
# レーン p の低レート系列（p だけシフトして M おきに取る）
xp = x[(M - p) % M::M][:L] if p != 0 else x[0::M][:L]
Y[:len(xp)] += np.convolve(xp, E[p])[:len(xp)] if False else 0 # placeholder
return Y
整流子の遅延の扱いは添字がややこしいので、ここでいったん止めて、添字を厳密に合わせた実装に書き直します。素朴実装と完全一致させるには、プロトタイプの遅延構造 $H(z)=\sum_p z^{-p}E_p(z^M)$ を素直に踏襲するのが安全です。
# ===== (b') ポリフェーズ間引き器（添字を厳密に対応させた版）=====
def polyphase_decimate(x, h, M):
N = len(h)
Lh = N // M
# E[p][m] = h[M*m + p]
E = np.array([h[p::M] for p in range(M)]) # 形状 (M, Lh)
Lx = len(x)
Lout = (Lx + N - 1) # full 畳み込み長（高レート相当）
# 各位相の入力サブシーケンス x[n-p] を M おきに取る（n = M*k）
# 出力 y[k] = sum_p sum_m E[p][m] * x[M*(k-m) - p]
Lk = Lx // M
y = np.zeros(Lk)
for p in range(M):
# x の位相 p 系列: インデックス -p, M-p, 2M-p, ... を 0 詰めで取り出す
xp = np.zeros(Lk + Lh)
for k in range(Lk + Lh):
idx = M*k - p
if 0 <= idx < Lx:
xp[k] = x[idx]
yp = np.convolve(xp, E[p])[:Lk] # 低レートで畳み込み
y += yp
return y
y_poly = polyphase_decimate(x, h, M)
L = min(len(y_naive), len(y_poly))
err = np.max(np.abs(y_naive[:L] - y_poly[:L]))
print(f"素朴実装とポリフェーズ実装の最大誤差: {err:.3e}")
このコードを実行すると、最大誤差が $10^{-12}$ 程度（倍精度の丸め誤差レベル）と表示されます。ポリフェーズ間引き器は素朴な「フィルタ→間引き」と数学的に等価であり、捨てるサンプルを計算しないだけで結果は完全に一致する、という導出どおりの事実が数値的に確認できました。差は浮動小数点の加算順序の違いだけです。
2. 演算量の削減を可視化する
次に、出力1サンプルあたりの積和回数が、素朴実装で $MN$、ポリフェーズで $N$（= $M$ 分の1）になることをグラフで確認します。
import numpy as np
import matplotlib.pyplot as plt
# チャネル数 M を変えたときの「出力1サンプルあたり積和回数」
K = 12 # チャネルあたりタップ数を固定
M_list = np.array([2, 4, 8, 16, 32, 64, 128, 256])
N_list = K * M_list # プロトタイプ長 N = K*M
# 単一チャネル間引き器
macs_naive_dec = M_list * N_list # MN
macs_poly_dec = N_list # N
# M チャネルフィルタバンク（出力1セットあたり）
macs_naive_fb = M_list * N_list # 各チャネル N タップ × M チャネル = MN
macs_dft_fb = N_list + M_list*np.log2(M_list) # N + M log2 M
plt.figure(figsize=(9, 6))
plt.loglog(M_list, macs_naive_dec, 'o-', label='Naive decimator: $MN$')
plt.loglog(M_list, macs_poly_dec, 's-', label='Polyphase decimator: $N$')
plt.loglog(M_list, macs_naive_fb, '^--', label='Naive $M$-ch bank: $MN$')
plt.loglog(M_list, macs_dft_fb, 'd--', label='DFT filter bank: $N + M\\log_2 M$')
plt.xlabel('Number of channels / decimation factor $M$')
plt.ylabel('MACs per output sample (set)')
plt.title('Computational cost: naive vs polyphase / DFT filter bank')
plt.grid(True, which='both', alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('polyphase_cost.png', dpi=150, bbox_inches='tight')
plt.show()
両対数グラフから3点が読み取れます。第1に、単一チャネルの間引きでは、ポリフェーズ（$N$）が素朴（$MN$）に対して常に $M$ 倍下にあり、$M$ が大きいほど縦の差（=削減倍率）が開きます。第2に、$M$ チャネルフィルタバンクで素朴に作ると $MN=KM^2$ で傾きが急（$M^2$ オーダー）に立ち上がるのに対し、DFTフィルタバンクは $N + M\log_2 M = KM + M\log_2 M$ でほぼ線形にしか増えません。第3に、$M=256$ では素朴フィルタバンクが約 $7.9\times10^5$ MACs なのに対しDFT版は約 $5\times10^3$ MACsで、2桁以上の差が付きます。導出した $O(MN)\to O(N+M\log M)$ がそのままコスト図に現れています。
3. DFTフィルタバンク：各チャネルの周波数応答
最後に、DFTフィルタバンクの各チャネル $H_k(z)=H(zW_M^k)$ の周波数応答を描き、$M$ 本のバンドパスが周波数軸を一様にカバーしていることを可視化します。あわせて、ポリフェーズ+IFFT実装が素朴な変調フィルタバンクと一致することも確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
M = 8
K = 12
N = K * M
h = signal.firwin(N, 1.0 / M, window=('kaiser', 8.0))
# 各チャネルフィルタ h_k[n] = h[n] * exp(j 2pi k n / M)
nidx = np.arange(N)
plt.figure(figsize=(10, 6))
worN = 2048
for k in range(M):
hk = h * np.exp(1j * 2*np.pi * k * nidx / M)
w, Hk = signal.freqz(hk, worN=worN, whole=True)
w_shift = np.where(w > np.pi, w - 2*np.pi, w)
order = np.argsort(w_shift)
plt.plot(w_shift[order]/np.pi,
20*np.log10(np.abs(Hk[order]) + 1e-12), linewidth=1.2,
label=f'ch {k}' if k < 4 else None)
plt.xlabel('Normalized frequency $\\omega/\\pi$')
plt.ylabel('Magnitude [dB]')
plt.title(f'Uniform DFT filter bank: {M} channel responses')
plt.ylim([-80, 5])
plt.grid(True, alpha=0.3)
plt.legend(ncol=2, fontsize=9)
plt.tight_layout()
plt.savefig('polyphase_channels.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、$M=8$ 本のバンドパスフィルタが $\omega$ 軸（$-\pi$ から $\pi$）上に等間隔 $2\pi/M$ で並び、隣り合うチャネルが約 $-3$ dB あたりでクロスオーバーしながら全帯域を一様に覆っていることが読み取れます。各チャネルは同じプロトタイプ $H$ を $2\pi k/M$ だけ周波数シフトしたものなので、形がすべて相似で、ストップバンド減衰（カイザー窓由来で約 $-60$ dB 以下）も共通です。これが「一様」DFTフィルタバンクと呼ばれる理由です。
import numpy as np
from scipy import signal
# ===== ポリフェーズ + IFFT 実装が素朴な変調フィルタバンクと一致するか =====
M, K = 8, 12
N = K * M
h = signal.firwin(N, 1.0 / M, window=('kaiser', 8.0))
rng = np.random.default_rng(1)
x = rng.standard_normal(2000)
# --- 素朴実装: 各チャネルで変調フィルタ→間引き ---
nidx = np.arange(N)
Y_naive = []
for k in range(M):
hk = h * np.exp(1j * 2*np.pi * k * nidx / M)
yk = signal.lfilter(hk, 1.0, x)[::M]
Y_naive.append(yk)
Y_naive = np.array(Y_naive) # 形状 (M, Lk)
# --- ポリフェーズ + IFFT 実装 ---
E = np.array([h[p::M] for p in range(M)]) # (M, K)
Lk = len(x) // M
U = np.zeros((M, Lk), dtype=complex)
for p in range(M):
xp = np.zeros(Lk + K)
for kk in range(Lk + K):
idx = M*kk - p
if 0 <= idx < len(x):
xp[kk] = x[idx]
U[p] = np.convolve(xp, E[p])[:Lk]
# 各時刻で M 点 IFFT（正規化を素朴側に合わせて M 倍）
Y_dft = np.fft.ifft(U, axis=0) * M
L = min(Y_naive.shape[1], Y_dft.shape[1])
err = np.max(np.abs(Y_naive[:, :L] - Y_dft[:, :L]))
print(f"素朴な変調FB と ポリフェーズ+IFFT の最大誤差: {err:.3e}")
出力される最大誤差は再び $10^{-12}$ 程度に収まります。$M$ 個の独立な変調バンドパスフィルタを別々に走らせた結果と、共通のポリフェーズ成分 $E_p$ を通してから1回の $M$ 点IFFTにかけた結果が、丸め誤差の範囲で完全に一致しました。導出した $Y_k=\sum_p W_M^{-kp}U_p$ がIDFTそのものであるという主張が、数値実験で裏付けられたことになります。$M$ 本のフィルタを共通化して1回のFFTにまとめても、出力は1ビットも犠牲にしていないのです。
まとめ
本記事では、ポリフェーズフィルタバンクの理論を、マルチレート処理の基礎からNoble恒等式、ポリフェーズ分解、DFTフィルタバンクへの帰着まで省略なく導出し、Pythonで効率化と各チャネル応答を確認しました。
- 間引きの無駄: 「フィルタ→間引き」は、これから捨てるサンプルまで高レートで計算しており、有効計算の $M$ 倍の積和を浪費している
- Noble恒等式: $\downarrow M$ の後の $G(z)$ は、前に置いた $G(z^M)$ と等価。$W_M^{-lM}=1$ で位相回転が消えることが本質で、フィルタを低レート側へ移せる
- ポリフェーズ分解: $H(z)=\sum_{p=0}^{M-1} z^{-p}E_p(z^M)$ と係数を $M$ 位相に仕分け、Noble恒等式で各 $E_p$ を間引き後の低レート側に移すと、出力1サンプルあたりの積和が $MN \to N$ に減る
- DFTフィルタバンク: 変調バンドパス $H_k(z)=H(zW_M^k)$ は共通の $E_p(z^M)$ と遅延を共有し、チャネルごとの違いは $W_M^{-kp}$ だけ。$Y_k=\sum_p W_M^{-kp}U_p$ はIDFTそのもので、IFFTで $M$ チャネルを一括計算できる
- 演算量: $M$ チャネルで $O(MN)\to O(N + M\log M)$。$N=KM$ なら $O(M^2)\to O(M\log M)$ に落ち、Pythonの数値実験でも素朴実装と丸め誤差レベルで一致しつつ2桁以上の削減を確認
ポリフェーズフィルタバンクは、SDRのチャネライザやOFDMの変復調を支える基盤技術であり、「無駄を計算しない」「共通構造をFFTに畳み込む」という2つの普遍的な発想の好例でもあります。次のステップとして、以下の記事も参考にしてください。