スマートフォンで録音した48 kHzの音声を、通信用に8 kHzへ落としてから送る。あるいは、44.1 kHzで記録されたCD音源を、48 kHzで動くオーディオインターフェースで再生する。こうした「サンプリング周波数を途中で変える」処理は、私たちが意識しないところで日常的に行われています。サンプリング周波数を下げたり上げたり、さらには44.1 kHzから48 kHzのような中途半端な比率で変換したりする一連の技術がマルチレート信号処理（multirate signal processing）です。
なぜわざわざレートを変えるのでしょうか。理由は大きく2つあります。第一に計算量とデータ量の削減です。低い周波数成分しか含まない信号を高いサンプリング周波数で持ち続けるのは、メモリも演算も無駄です。必要十分なレートまで落とせば、後段のフィルタ処理やストレージのコストが劇的に下がります。第二に異なるレートのシステムを橋渡しするためです。サンプリング周波数の異なる機器やプロトコルを接続するには、レート変換が避けて通れません。
マルチレート処理は応用範囲も広く、たとえば次のような場面で本質的な役割を果たします。
- ソフトウェア無線（SDR）と通信: 広帯域でサンプリングした受信信号から目的チャネルだけを抜き出すとき、ダウンサンプリングで処理レートを落とします。送信側ではアップサンプリングと整形フィルタでシンボルをきれいな波形に成形します。
- オーディオとマルチメディア: サンプリングレート変換（SRC）、サブバンド符号化（MP3やAACの基礎であるフィルタバンク）、オーバーサンプリングDAC/ADCなど、ほぼすべての音響処理がマルチレートに支えられています。
本記事では、間引き（デシメーション）とゼロ挿入（補間）がスペクトルにどう作用するかを数式で導き、なぜアンチエイリアスフィルタや補間フィルタが必要なのかを理解します。さらに $L/M$ 倍の有理数レート変換、そして乗算回数を $1/M$ に削減するポリフェーズ分解の原理（Noble恒等式）まで踏み込みます。
本記事の内容
- ダウンサンプリング（間引き）のスペクトルへの作用と折り返しの導出
- デシメーションフィルタ（アンチエイリアスフィルタ）が必要な理由
- アップサンプリング（ゼロ挿入）のイメージスペクトルと補間フィルタ
- $L/M$ 倍の有理数レート変換のカスケード構成
- ポリフェーズ分解とNoble恒等式による効率化の原理
- Pythonでのスペクトル可視化とポリフェーズFIRの計算量比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、離散時間信号のスペクトル（DTFT）が $2\pi$ 周期を持つこと、サンプリング定理とエイリアシング、FIRフィルタの周波数応答について理解していると、本記事の議論がスムーズに追えます。
マルチレートとは — 「コマ送り」と「コマ間挿入」のイメージ
マルチレート処理の2つの基本操作は、映画のフィルムにたとえると直感的に理解できます。
ダウンサンプリング（間引き、decimation）は、フィルムを「コマ送り」で見るようなものです。1秒30コマで撮影した映像を、3コマに1コマだけ抜き出して10コマにすれば、データ量は1/3になります。ただし、速く動く被写体は「飛び飛び」にしか見えなくなり、本来の滑らかな動きが失われます。これが後で説明するエイリアシングに対応します。
アップサンプリング（補間、interpolation）はその逆で、コマとコマの「あいだ」に新しいコマを挿入する操作です。ただし、ただ挿入するだけでは何を入れればよいか分かりません。素朴に「真っ黒なコマ（ゼロ）」を入れると映像はチラついてしまいます。このチラつきを取り除いて滑らかな中間コマを作るのが補間フィルタの役割です。
このように、マルチレート処理は単に「サンプルを抜く・足す」だけでは済まず、必ずフィルタとの組み合わせで初めて意味を持ちます。なぜフィルタが必要なのかを理解するには、間引きやゼロ挿入が周波数領域（スペクトル）でどんな作用をするかを正確に把握する必要があります。まずはダウンサンプリングから、その数式を一歩ずつ導いていきましょう。
ダウンサンプリング（間引き）の数学
間引き器の定義
$M$ を正の整数とします。$M$ 倍のダウンサンプラ（down-sampler）は、入力 $x[n]$ から $M$ サンプルおきに1つだけ取り出して出力 $y[n]$ を作る操作です。
$$ y[n] = x[Mn] $$
たとえば $M=3$ なら、$x[0], x[3], x[6], \dots$ だけを残し、間の $x[1], x[2], x[4], x[5], \dots$ を捨てます。サンプリング周波数で言えば $f_s$ が $f_s/M$ になります。
この操作は一見単純ですが、スペクトルには劇的な変化を引き起こします。サンプルを捨てることが、周波数領域でどう見えるのか。これを正確に知ることが、マルチレート処理を理解する第一歩です。
スペクトルの導出
入力 $x[n]$ のDTFT（離散時間フーリエ変換）を $X(e^{j\omega})$ とします。
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} $$
出力 $y[n] = x[Mn]$ のDTFT $Y(e^{j\omega})$ を $X(e^{j\omega})$ で表すのがゴールです。定義に従って書き下すと
$$ Y(e^{j\omega}) = \sum_{n=-\infty}^{\infty} y[n] e^{-j\omega n} = \sum_{n=-\infty}^{\infty} x[Mn] e^{-j\omega n} $$
ここで $m = Mn$ と置きたいところですが、$m$ は $M$ の倍数だけを動くため、そのままでは $x$ のすべてのサンプルを足し上げる和に書き換えられません。そこで、$M$ の倍数だけを拾い出す「くし（インパルス列）」を使う巧妙な手法を導入します。
任意の整数 $m$ に対して、次の和が「$m$ が $M$ の倍数なら1、そうでなければ0」を返すことを利用します。
$$ \frac{1}{M} \sum_{k=0}^{M-1} e^{j 2\pi k m / M} = \begin{cases} 1 & m \equiv 0 \pmod{M} \\ 0 & \text{otherwise} \end{cases} $$
これは等比級数の和で確かめられます。$m$ が $M$ の倍数なら各項が $1$ なので和は $M$、$1/M$ を掛けて $1$。$m$ が $M$ の倍数でなければ公比 $e^{j2\pi m/M} \neq 1$ の等比級数で、$M$ 項の和は $\frac{1 – e^{j2\pi m}}{1 – e^{j2\pi m/M}} = 0$ となります。この「ふるい」を使うと、$M$ の倍数だけを残す処理を全インデックスの和として書けます。
まず $y[n] = x[Mn]$ を、いったん全インデックス $m$ にわたる和へと「水増し」します。上のふるいを使えば
$$ \sum_{n=-\infty}^{\infty} x[Mn] e^{-j\omega n} = \sum_{m=-\infty}^{\infty} x[m] \left( \frac{1}{M} \sum_{k=0}^{M-1} e^{j 2\pi k m / M} \right) e^{-j\omega m / M} $$
と書けます。ここで $m=Mn$ のとき $e^{-j\omega n} = e^{-j\omega m/M}$ であること、そしてふるいが $m$ が $M$ の倍数のときだけ $1$ を返すことを使いました。$k$ の和と $m$ の和の順序を入れ替えると
$$ Y(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} \sum_{m=-\infty}^{\infty} x[m] \, e^{j 2\pi k m / M} e^{-j\omega m / M} $$
指数をまとめると $e^{-j m (\omega – 2\pi k)/M}$ となります。内側の $m$ の和はまさにDTFTの定義式そのもので、角周波数が $(\omega – 2\pi k)/M$ のときの $X$ にほかなりません。したがって
$$ \boxed{\; Y(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} X\!\left( e^{\,j (\omega – 2\pi k)/M} \right) \;} $$
これがダウンサンプリングの基本公式です。導出を振り返ると、「ふるい」を使って間引きを全インデックスの和に書き換え、和の順序を入れ替えただけで自然にこの形に至りました。
公式が語ること — 周波数軸の伸長と重畳
導いた公式の意味を読み解きましょう。$Y(e^{j\omega})$ は $M$ 個の項の和です。各項を見ると、
- $X(e^{j\omega/M})$ という形は、もとのスペクトルを周波数軸方向に $M$ 倍に引き伸ばしたものです（$\omega$ の代わりに $\omega/M$ を入れるので、横軸が $M$ 倍に拡大される）。
- $2\pi k$ のシフトは、引き伸ばしたスペクトルのコピーを $k=0,1,\dots,M-1$ の分だけ $2\pi$ ごとにずらして重ねることを意味します。
- 全体に $1/M$ の振幅スケーリングがかかります。
つまりダウンサンプリングは、スペクトルを横に $M$ 倍引き伸ばし、$M$ 個のコピーを足し合わせる操作です。引き伸ばしによって、もとは $[-\pi/M, \pi/M]$ に収まっていた帯域が $[-\pi, \pi]$ いっぱいに広がります。問題は、もし $x[n]$ が $|\omega| > \pi/M$ の成分を持っていると、引き伸ばされたコピー同士が重なり合い、互いに足し合わされて分離できなくなる点です。これが折り返し（エイリアシング）です。
この「重なって混ざる」という現象こそ、間引きの前にフィルタが必要になる根本的な理由です。次の節で、なぜ事前にフィルタをかけねばならないのかを具体的に見ていきます。
デシメーションフィルタはなぜ必要か
エイリアシングの発生条件
ダウンサンプリングの公式から、$M$ 個のスペクトルコピーが重ならない条件を考えます。$k=0$ の項は $\omega/M$ で引き伸ばされたメインのスペクトル、$k=1$ の項は $2\pi$ だけずれたコピーです。これらが $[-\pi, \pi]$ の範囲で重ならないためには、もとの信号 $x[n]$ の帯域が
$$ |\omega| < \frac{\pi}{M} $$
に制限されている必要があります。これは新しいナイキスト周波数（レート $f_s/M$ における $\pi$）を、もとの周波数軸に戻したものに対応します。言い換えると、間引き後のナイキスト周波数 $f_s/(2M)$ を超える成分があると、それが低い周波数へ「折り返し」、本来の信号に混ざってしまいます。
たとえば $f_s = 48$ kHz の信号を $M=6$ で間引いて $8$ kHz にする場合、新しいナイキスト周波数は $4$ kHz です。もし元信号に $5$ kHz の成分があれば、それは $8 – 5 = 3$ kHz の偽の信号として現れ、二度と取り除けません。
アンチエイリアスフィルタ + 間引き = デシメーション
この折り返しを防ぐには、間引きの前にローパスフィルタをかけて、$|\omega| \geq \pi/M$ の成分をあらかじめ除去しておけばよいのです。このフィルタをアンチエイリアスフィルタ（anti-aliasing filter）またはデシメーションフィルタと呼びます。理想的にはカットオフ $\omega_c = \pi/M$ のローパスフィルタです。
「アンチエイリアスフィルタ + ダウンサンプラ」をひとまとめにした処理をデシメーション（decimation）と呼びます。構成は次のとおりです。
$$ x[n] \;\longrightarrow\; \boxed{H(z),\ \omega_c = \pi/M} \;\longrightarrow\; v[n] \;\longrightarrow\; \boxed{\downarrow M} \;\longrightarrow\; y[n] $$
フィルタリングで帯域を $\pi/M$ 以内に絞ってから間引けば、引き伸ばされたコピーが重ならず、$Y(e^{j\omega}) = \frac{1}{M} X(e^{j\omega/M})$（$|\omega| < \pi$）という、折り返しのないきれいなスペクトルが得られます。
ここで注意したいのは、デシメーションでは「帯域外の情報を捨てている」という事実です。間引きによって低周波成分を温存する代わりに、高周波成分は永久に失われます。これは情報の損失ですが、目的の信号が低周波数帯に収まっているなら問題なく、むしろデータ量削減という恩恵だけを受けられます。
「捨てる」のがデシメーションなら、「増やす」のはどうなるのか。次はその逆操作であるアップサンプリングを見ていきましょう。スペクトルにはまた別の、しかし鏡のような現象が現れます。
アップサンプリング（ゼロ挿入）の数学
ゼロ挿入器の定義
$L$ を正の整数とします。$L$ 倍のアップサンプラ（up-sampler）は、入力 $x[n]$ のサンプルとサンプルのあいだに $L-1$ 個のゼロを挿入する操作です。
$$ v[n] = \begin{cases} x[n/L] & n \text{ が } L \text{ の倍数} \\ 0 & \text{otherwise} \end{cases} $$
たとえば $L=3$ なら、$x[0], 0, 0, x[1], 0, 0, x[2], \dots$ という列を作ります。サンプル数は $L$ 倍になり、サンプリング周波数は $f_s$ から $L f_s$ になります。
ゼロを挟むだけなので情報は失われませんが、スペクトルには「コピーが増える」という現象が起きます。これを次に導出します。
スペクトルの導出
出力 $v[n]$ のDTFTを計算します。$v[n]$ は $L$ の倍数のインデックスでのみ非ゼロなので、$n = Lm$ とおくと和はこれらの項だけ残ります。
$$ V(e^{j\omega}) = \sum_{n=-\infty}^{\infty} v[n] e^{-j\omega n} = \sum_{m=-\infty}^{\infty} v[Lm] e^{-j\omega L m} $$
ここで $v[Lm] = x[m]$ なので
$$ V(e^{j\omega}) = \sum_{m=-\infty}^{\infty} x[m] e^{-j\omega L m} = X(e^{j\omega L}) $$
非常にシンプルな結果です。
$$ \boxed{\; V(e^{j\omega}) = X(e^{j\omega L}) \;} $$
導出はダウンサンプリングよりずっと短く済みました。ゼロ挿入は単に「インデックスのスケーリング」だけなので、ふるいのような工夫は不要だったのです。
公式が語ること — イメージの出現
$V(e^{j\omega}) = X(e^{j\omega L})$ は何を意味するでしょうか。$\omega$ を $L$ 倍してから $X$ に入れるので、もとのスペクトルは周波数軸方向に $1/L$ に圧縮されます。DTFTは $2\pi$ 周期ですから、もとのスペクトルが $[-\pi, \pi]$ にあったとすると、圧縮後は $[-\pi/L, \pi/L]$ に縮まり、その外側の $[-\pi, \pi]$ には同じスペクトルのコピーが $L$ 個並びます。
このコピーをイメージ（image）と呼びます。具体的には、もとの信号成分が周波数 $\omega_0$ にあったとすると、アップサンプリング後には $\omega_0/L$ に本来の成分が現れ、さらに $\frac{2\pi k \pm \omega_0}{L}$（$k=1,\dots,L-1$）の位置にイメージが現れます。これらのイメージは、ゼロを挿入したことで波形が「ギザギザ」になった結果として生じる、高周波の偽スペクトルです。
イメージをそのまま放置すると、滑らかに補間された信号になりません。コマ送りの例で言えば、真っ黒なコマを挟んだままチラついている状態です。これを除くフィルタが必要になります。
補間フィルタ + ゼロ挿入 = インターポレーション
イメージを取り除くには、ゼロ挿入の後にカットオフ $\omega_c = \pi/L$ のローパスフィルタをかけ、$[-\pi/L, \pi/L]$ のメイン成分だけを残せばよいのです。このフィルタを補間フィルタ（interpolation filter）と呼びます。「ゼロ挿入 + 補間フィルタ」をまとめた処理がインターポレーション（interpolation）です。
$$ x[n] \;\longrightarrow\; \boxed{\uparrow L} \;\longrightarrow\; v[n] \;\longrightarrow\; \boxed{H(z),\ \omega_c = \pi/L} \;\longrightarrow\; y[n] $$
補間フィルタは、挿入したゼロを「埋める」役割も果たします。ローパスフィルタはゼロの位置に隣接サンプルから内挿した値を生成するので、結果として滑らかな高レート信号が得られます。
ここで一つ注意があります。補間フィルタのゲインは通常 $L$ に設定します。なぜなら、ゼロ挿入によってエネルギーが $L$ 個のイメージに分散して振幅が $1/L$ になるため、メイン成分の振幅をもとに戻すには $L$ 倍のゲインが必要だからです。デシメーションフィルタのゲインが $1$ でよかったのと対照的です。
デシメーションとインターポレーションが理解できれば、両者を組み合わせて任意の有理数倍のレート変換ができます。次はその仕組みを見ていきましょう。
有理数倍 $L/M$ のレート変換
なぜカスケードが必要か
整数倍のダウン・アップは扱えるようになりましたが、現実には $44.1$ kHz から $48$ kHz への変換のように、整数比でないレート変換が頻繁に必要です。この比は
$$ \frac{48000}{44100} = \frac{480}{441} = \frac{160}{147} $$
と既約分数 $L/M = 160/147$ で表せます。つまり「$160$ 倍にアップサンプリングしてから $147$ 倍にダウンサンプリングする」ことで、ちょうど $160/147$ 倍のレート変換が実現できます。
一般に、レートを $L/M$ 倍（$L, M$ は互いに素な正の整数）にするには、$L$ 倍アップサンプリング、フィルタリング、$M$ 倍ダウンサンプリングをこの順で行います。
順序が重要 — アップが先、ダウンが後
ここで決定的に重要なのが順序です。必ず「アップサンプリング → ダウンサンプリング」の順にします。逆にすると、先にダウンサンプリングで情報を捨ててしまい、その後にアップサンプリングしても失われた情報は戻りません。アップを先にすればサンプル点が増えるだけで情報は失われず、その後のダウンで安全に間引けます。
1つの共通フィルタに統合する
素朴に考えると、アップサンプリングの後の補間フィルタ（$\omega_c = \pi/L$）と、ダウンサンプリングの前のアンチエイリアスフィルタ（$\omega_c = \pi/M$）の2つが必要に見えます。しかし、この2つのフィルタは直列に接続された2つのローパスフィルタであり、間に他の処理がないので1つのローパスフィルタにまとめられます。
まとめたフィルタのカットオフは、2つのカットオフのうち厳しいほう（小さいほう）を採用します。
$$ \omega_c = \min\!\left( \frac{\pi}{L}, \frac{\pi}{M} \right) = \frac{\pi}{\max(L, M)} $$
ゲインは補間フィルタの分の $L$ を採用します。最終的な $L/M$ レート変換器の構成は次のようになります。
$$ x[n] \;\to\; \boxed{\uparrow L} \;\to\; \boxed{H(z),\ \omega_c = \frac{\pi}{\max(L,M)},\ \text{gain}=L} \;\to\; \boxed{\downarrow M} \;\to\; y[n] $$
この構成は理屈の上では完璧ですが、効率の点で大きな無駄を含んでいます。次の節でその無駄を明らかにし、ポリフェーズ分解という強力な解決策を導きます。
カスケード構成の非効率性
たとえば $L/M = 160/147$ の変換を考えます。$L=160$ 倍にアップサンプリングすると、信号のサンプル数が一気に $160$ 倍になります。その高レート信号に対してフィルタリングを行い、その後 $M=147$ で間引いて結局ほとんどを捨てます。
ここには2つの無駄があります。第一に、アップサンプリングで挿入したゼロにもフィルタの乗算を行っている点です。ゼロに係数を掛けても結果はゼロなので、この乗算は完全に無駄です。第二に、ダウンサンプリングで捨てるサンプルについてもフィルタ出力を計算している点です。捨てるなら計算する意味がありません。
これらの無駄を取り除くのがポリフェーズ分解です。考え方を理解すれば、計算量を劇的に削減できます。
ポリフェーズ分解による効率化
基本アイデア — ゼロとの乗算をなくす
ポリフェーズ分解の核心は、「ゼロに掛ける乗算をそもそも実行しない」という発想です。アップサンプリングで挿入されたゼロは規則的に並んでいるので、フィルタ係数を $M$（または $L$）個のグループに分けておけば、各時刻で実際に非ゼロのサンプルに掛かる係数だけを選んで計算できます。
これを数式で表現するために、まずフィルタの伝達関数 $H(z)$ をポリフェーズ成分に分解します。
ポリフェーズ分解の定義
長さ $N$ のFIRフィルタ $h[n]$ の伝達関数を
$$ H(z) = \sum_{n=0}^{N-1} h[n] z^{-n} $$
とします。これを $M$ 個のポリフェーズ成分に分けます。インデックス $n$ を $n = Mr + \rho$（$\rho = 0, 1, \dots, M-1$ は剰余、$r$ は商）と書き分けると、和を $\rho$ ごとにグループ化できます。
$$ H(z) = \sum_{\rho=0}^{M-1} \sum_{r} h[Mr + \rho] \, z^{-(Mr+\rho)} = \sum_{\rho=0}^{M-1} z^{-\rho} \sum_{r} h[Mr+\rho] \, z^{-Mr} $$
ここで、$\rho$ 番目のポリフェーズ成分（サブフィルタ）を
$$ E_\rho(z) = \sum_{r} h[Mr + \rho] \, z^{-r} $$
と定義すると、$H(z)$ は次のように $M$ 個のサブフィルタの和で書けます。
$$ \boxed{\; H(z) = \sum_{\rho=0}^{M-1} z^{-\rho} \, E_\rho(z^M) \;} $$
各 $E_\rho$ は、もとのフィルタ係数を $M$ 個おきに拾い集めたものなので、長さは約 $N/M$ です。$M$ 個のサブフィルタを合わせれば全係数を過不足なくカバーします。つまり、フィルタを「位相（phase）」ごとに分解したわけで、これが「ポリフェーズ（多相）」の名の由来です。
この分解だけでは計算量は変わりません。威力を発揮するのは、次に説明するNoble恒等式と組み合わせたときです。
Noble恒等式
ポリフェーズの真価は、ダウンサンプラ／アップサンプラとフィルタの順序を交換することにあります。その根拠となるのがNoble恒等式（Noble identities）です。
ダウンサンプリングのNoble恒等式: 「フィルタ $G(z^M)$ をかけてから $M$ 倍ダウンサンプリングする」ことは、「先に $M$ 倍ダウンサンプリングしてからフィルタ $G(z)$ をかける」ことに等しい。
$$ \boxed{\; \big[\,\downarrow M\,\big] \circ G(z^M) = G(z) \circ \big[\,\downarrow M\,\big] \;} $$
アップサンプリングのNoble恒等式: 「$L$ 倍アップサンプリングしてからフィルタ $G(z^L)$ をかける」ことは、「先にフィルタ $G(z)$ をかけてから $L$ 倍アップサンプリングする」ことに等しい。
$$ \boxed{\; G(z^L) \circ \big[\,\uparrow L\,\big] = \big[\,\uparrow L\,\big] \circ G(z) \;} $$
ポイントは、フィルタの引数が $z^M$（または $z^L$）という「$M$ 倍に間延びした」形になっている場合に限り、間引き／ゼロ挿入と交換できるという点です。$G(z^M)$ は係数のあいだに $M-1$ 個のゼロが挟まったフィルタなので、間引きと自然に整合するのです。
なぜこれが効率化につながるのかを、デシメーションの場合で見てみましょう。
デシメーションのポリフェーズ構成
$M$ 倍デシメーション（フィルタ $H(z)$ → $\downarrow M$）を考えます。$H(z)$ をポリフェーズ分解すると
$$ y[n] = \big[\,\downarrow M\,\big]\!\left( \sum_{\rho=0}^{M-1} z^{-\rho} E_\rho(z^M) \, X(z) \right) $$
各 $E_\rho(z^M)$ にダウンサンプリングのNoble恒等式を適用すると、$\downarrow M$ をフィルタの内側に移動できます。すると各枝は「先に $\downarrow M$、その後にサブフィルタ $E_\rho(z)$」という形になります。
$$ Y(z) = \sum_{\rho=0}^{M-1} E_\rho(z) \cdot \big[\,\downarrow M\,\big]\!\big( z^{-\rho} X(z) \big) $$
この構成の決定的な利点は、サブフィルタ $E_\rho(z)$ が低レート側（間引き後）で動作する点です。もとの構成ではフルレートでフィルタリングしてから間引いていたのに対し、ポリフェーズ構成では先に間引いてから各サブフィルタを動かします。つまり捨てるサンプルの計算を一切行いません。
計算量が $1/M$ になる理由
具体的に乗算回数を比較しましょう。タップ数 $N$ のフィルタで $M$ 倍デシメーションする場合を考えます。
素朴な構成（先にフィルタ、後で間引き）: フルレートの全サンプルについて $N$ 回の乗算を行い、その後 $M$ サンプルに $1$ つだけ残します。出力1サンプルあたりの実効乗算回数は $N \cdot M / M = N$… ではなく、出力サンプルを得るのに必要な乗算は $N$ 回ですが、フルレートで計算してしまう実装だと、入力サンプルあたり $N$ 回、すなわち出力サンプルあたり $N \times M$ 回計算してしまう無駄が生じます。捨てるサンプルの出力も計算しているからです。
ポリフェーズ構成: 各サブフィルタ $E_\rho$ のタップ数は約 $N/M$ で、$M$ 個のサブフィルタの出力を合計します。出力1サンプルを得るのに必要な乗算は $M \times (N/M) = N$ 回です。これは「捨てるサンプルを計算しない」分、素朴なフルレート実装に比べて $1/M$ の計算量です。
要するに、ポリフェーズ分解は「アップサンプリングで生じるゼロとの乗算」と「ダウンサンプリングで捨てるサンプルの計算」という2つの無駄を構造的に排除します。これにより、計算は常に低いほうのレートで行われ、乗算回数が劇的に削減されます。
理論はここまでです。次は実際にPythonで信号を作り、間引き・補間がスペクトルにどう作用するかを目で確認し、ポリフェーズFIRを実装して計算量の削減を測定しましょう。
Pythonでの実装
ダウンサンプリングのスペクトルとエイリアシング
まず、フィルタなしで素朴に間引くとエイリアシングが起きること、そしてアンチエイリアスフィルタをかければ防げることを可視化します。元信号として、低周波の正弦波と、新ナイキストを超える高周波の正弦波を混ぜたものを使います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# サンプリング設定
fs = 12000 # 元のサンプリング周波数 [Hz]
M = 4 # 間引き率（新レートは 3000 Hz, 新ナイキスト 1500 Hz）
N = 4096 # サンプル数
t = np.arange(N) / fs
# 元信号: 500 Hz（残したい）と 2500 Hz（新ナイキスト 1500 Hz を超える）
x = np.sin(2*np.pi*500*t) + 0.8*np.sin(2*np.pi*2500*t)
# (A) フィルタなしで素朴に間引く
y_naive = x[::M]
# (B) アンチエイリアスフィルタをかけてから間引く
# 新ナイキスト fs/(2M) に合わせてカットオフ pi/M
h_aa = signal.firwin(101, 1.0/M) # カットオフ = 1/M（ナイキスト基準の正規化）
x_filt = signal.lfilter(h_aa, 1.0, x)
y_dec = x_filt[::M]
print(f"元レート: {fs} Hz, 間引き後レート: {fs//M} Hz, 新ナイキスト: {fs//(2*M)} Hz")
この時点で、2500 Hzの成分は間引き後の新ナイキスト1500 Hzを超えています。理論によれば、フィルタなしで間引くとこの成分が $3000 – 2500 = 500$ Hz に折り返し、もとの500 Hz成分と混ざってしまうはずです。次にスペクトルを描いて確認します。
def spectrum(sig, fs_local):
"""片側振幅スペクトルを返す"""
Nlen = len(sig)
f = np.fft.rfftfreq(Nlen, 1/fs_local)
amp = np.abs(np.fft.rfft(sig)) / Nlen * 2
return f, amp
f_in, A_in = spectrum(x, fs)
f_naive, A_naive = spectrum(y_naive, fs//M)
f_dec, A_dec = spectrum(y_dec, fs//M)
fig, axes = plt.subplots(3, 1, figsize=(10, 9))
axes[0].plot(f_in, A_in, 'b')
axes[0].set_title('Input spectrum (fs=12 kHz): 500 Hz + 2500 Hz')
axes[0].axvline(1500, color='r', ls='--', alpha=0.6, label='new Nyquist 1500 Hz')
axes[0].legend(); axes[0].set_xlim(0, fs/2)
axes[1].plot(f_naive, A_naive, 'm')
axes[1].set_title('Naive downsample by 4 (NO filter) -> aliasing at 500 Hz')
axes[1].set_xlim(0, fs//M/2)
axes[2].plot(f_dec, A_dec, 'g')
axes[2].set_title('Decimation by 4 (with anti-alias filter)')
axes[2].set_xlim(0, fs//M/2)
for ax in axes:
ax.set_xlabel('Frequency [Hz]'); ax.set_ylabel('Amplitude'); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('downsample_aliasing.png', dpi=150, bbox_inches='tight')
plt.show()
3枚のグラフから、理論どおりの結果が読み取れます。一番上の入力スペクトルには500 Hzと2500 Hzの2本のピークがはっきり立っています。真ん中のフィルタなしで間引いた結果では、2500 Hz成分が消えた代わりに500 Hz付近のピークが膨らんでおり、折り返し成分がもとの500 Hzに重なってしまったことが分かります。これがエイリアシングです。一番下のアンチエイリアスフィルタ付きデシメーションでは、2500 Hz成分が間引き前に除去されているため、500 Hzのピークだけがきれいに残り、折り返しが起きていません。フィルタが折り返しを確実に防いでいることが視覚的に確認できます。
アップサンプリングのスペクトルとイメージ
次に、ゼロ挿入によってイメージが現れること、そして補間フィルタでそれを除けることを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 2000 # 元のサンプリング周波数 [Hz]
L = 4 # アップサンプリング率（新レート 8000 Hz）
N = 2048
t = np.arange(N) / fs
x = np.sin(2*np.pi*300*t) # 300 Hz の正弦波
# (A) ゼロ挿入のみ（イメージが残る）
v = np.zeros(N*L)
v[::L] = x # L サンプルごとに元サンプルを置く
# (B) ゼロ挿入 + 補間フィルタ（ゲイン L、カットオフ pi/L）
h_interp = L * signal.firwin(101, 1.0/L) # ゲインを L 倍してエネルギーを補償
y_interp = signal.lfilter(h_interp, 1.0, v)
print(f"元レート: {fs} Hz -> アップ後レート: {fs*L} Hz")
ゼロ挿入したベクトル v のサンプリング周波数は $L f_s = 8000$ Hz です。理論によれば、300 Hzの成分のほかに、$8000 – 300 = 7700$ Hz などのイメージが新しい帯域に現れるはずです。補間フィルタのゲインを $L=4$ 倍にしているのは、ゼロ挿入で振幅が $1/L$ に下がるのを補うためです。スペクトルで確認します。
def spectrum(sig, fs_local):
Nlen = len(sig)
f = np.fft.rfftfreq(Nlen, 1/fs_local)
amp = np.abs(np.fft.rfft(sig)) / Nlen * 2
return f, amp
f_v, A_v = spectrum(v, fs*L)
f_y, A_y = spectrum(y_interp, fs*L)
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
axes[0].plot(f_v, A_v, 'm')
axes[0].set_title('Zero-insertion only (fs=8 kHz): images appear')
for k in range(1, L):
axes[0].axvline(k*fs - 300, color='r', ls='--', alpha=0.4)
axes[0].axvline(k*fs + 300, color='r', ls='--', alpha=0.4)
axes[1].plot(f_y, A_y, 'g')
axes[1].set_title('Interpolation (zero-insert + LPF, gain=L): images removed')
for ax in axes:
ax.set_xlim(0, fs*L/2); ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Amplitude'); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('upsample_imaging.png', dpi=150, bbox_inches='tight')
plt.show()
上のグラフでは、300 Hzの本来の成分に加えて、1700 Hz、2300 Hz、3700 Hzといった位置に複数のイメージ（赤破線の位置）が現れています。これがゼロ挿入によって生じる偽のスペクトルです。下のグラフでは補間フィルタによってこれらのイメージがすべて除去され、300 Hzの成分だけが残り、しかもゲイン補償によって振幅がもとの正弦波と同じ約1.0に復元されています。ゼロ挿入だけでは不十分で、補間フィルタが滑らかな高レート信号を作るために不可欠であることが分かります。
時間波形で見る補間効果
スペクトルだけでなく、時間波形でも補間フィルタが「ゼロを埋めて滑らかにする」様子を見てみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 50
L = 8
N = 16
t = np.arange(N) / fs
x = np.sin(2*np.pi*3*t) # 低周波正弦波（粗いサンプリング）
v = np.zeros(N*L); v[::L] = x # ゼロ挿入
h = L * signal.firwin(8*L+1, 1.0/L) # 補間フィルタ
y = signal.lfilter(h, 1.0, v)
delay = (len(h)-1)//2 # FIRの群遅延を補正
t_up = np.arange(N*L) / (fs*L)
fig, ax = plt.subplots(figsize=(10, 4))
ax.stem(t, x, linefmt='C0-', markerfmt='C0o', basefmt=' ', label='original samples')
ax.plot(t_up - delay/(fs*L), y, 'g-', label='interpolated (smooth)')
ax.set_xlabel('Time [s]'); ax.set_ylabel('Amplitude')
ax.set_title('Interpolation reconstructs a smooth waveform between samples')
ax.legend(); ax.grid(alpha=0.3); ax.set_xlim(0, 0.25)
plt.tight_layout()
plt.savefig('interpolation_time.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフでは、もとの粗いサンプル（青い点）のあいだを、補間後の信号（緑の曲線）が滑らかな正弦波としてつないでいる様子が見えます。挿入したゼロは補間フィルタによって周囲のサンプルから内挿された値に置き換わり、もとの連続信号の形が高い解像度で復元されています。補間が単なるゼロ詰めではなく「サンプル間の値を推定する処理」であることが直感的に理解できます。
有理数倍 $L/M$ のレート変換
次に、44.1 kHzから48 kHzへの変換を scipy.signal.resample_poly で実行します。この関数は内部でポリフェーズ分解を使っており、効率的に有理数倍変換を行います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from math import gcd
fs_in = 44100
fs_out = 48000
g = gcd(fs_in, fs_out)
L, M = fs_out // g, fs_in // g # 既約分数 L/M
print(f"L/M = {L}/{M} -> {fs_in} Hz x {L}/{M} = {fs_in*L/M:.1f} Hz")
# テスト信号: 1 kHz + 5 kHz
dur = 0.05
t_in = np.arange(int(fs_in*dur)) / fs_in
x = np.sin(2*np.pi*1000*t_in) + 0.5*np.sin(2*np.pi*5000*t_in)
# ポリフェーズによる有理数倍リサンプリング
y = signal.resample_poly(x, L, M)
t_out = np.arange(len(y)) / fs_out
print(f"入力サンプル数: {len(x)}, 出力サンプル数: {len(y)}")
print(f"比の確認: {len(y)/len(x):.5f} (理論値 {L/M:.5f})")
出力された $L/M$ は $160/147$ となり、44.1 kHzを $160/147$ 倍すれば正確に48 kHzになることが確認できます。出力サンプル数も入力のおよそ $160/147$ 倍になっています。レート変換が正しい比率で行われたことが、サンプル数の比からも裏づけられます。波形とスペクトルを確認します。
def spectrum(sig, fs_local):
Nlen = len(sig)
f = np.fft.rfftfreq(Nlen, 1/fs_local)
amp = np.abs(np.fft.rfft(sig)) / Nlen * 2
return f, amp
f_in, A_in = spectrum(x, fs_in)
f_out, A_out = spectrum(y, fs_out)
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
axes[0].plot(f_in, A_in, 'b'); axes[0].set_title('Input spectrum @ 44.1 kHz')
axes[0].set_xlim(0, 8000)
axes[1].plot(f_out, A_out, 'g'); axes[1].set_title('Resampled spectrum @ 48 kHz (160/147)')
axes[1].set_xlim(0, 8000)
for ax in axes:
ax.set_xlabel('Frequency [Hz]'); ax.set_ylabel('Amplitude'); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('rational_resample.png', dpi=150, bbox_inches='tight')
plt.show()
2つのスペクトルを比べると、入力にあった1 kHzと5 kHzのピークが、出力でも同じ周波数位置・同じ振幅で保たれていることが分かります。サンプリング周波数は変わったのに、信号そのものの周波数成分は変化していません。これがレート変換の正しい振る舞いで、イメージやエイリアシングが現れていないことから、内部の共通フィルタが適切に働いていることも確認できます。
ポリフェーズFIRの自作実装と計算量比較
最後に、ポリフェーズ分解の効果を測るため、デシメーションを「素朴なフルレート実装」と「ポリフェーズ実装」の両方で書き、出力の一致と乗算回数を比較します。
import numpy as np
from scipy import signal
def decimate_naive(x, h, M):
"""素朴な実装: フルレートで全サンプルをフィルタしてから間引く"""
y_full = np.convolve(x, h) # フルレートでフィルタ（全サンプル計算）
y = y_full[::M] # 間引き（大半を捨てる）
mults = len(x) * len(h) # フルレートでの乗算回数
return y, mults
def decimate_polyphase(x, h, M):
"""ポリフェーズ実装: 間引き後の低レートでサブフィルタを動かす"""
N = len(h)
# フィルタ係数を M 個のポリフェーズ成分に分解
h_pad = np.concatenate([h, np.zeros((-N) % M)]) # 長さを M の倍数に揃える
E = h_pad.reshape(-1, M).T # E[rho] が rho 番目のサブフィルタ
L_out = (len(x) + M - 1) // M
# 各枝の入力は x を rho サンプル遅延させてから M 間引きしたもの（= x[Mk-rho]）
xx = np.concatenate([np.zeros(M), x]) # 遅延分のゼロを前置
y = np.zeros(L_out + E.shape[1])
mults = 0
for rho in range(M):
xr = xx[M::M] if rho == 0 else xx[M-rho::M] # 枝rhoの入力 x[Mk-rho]
sub = np.convolve(xr, E[rho])
y[:len(sub)] += sub
mults += len(xr) * len(E[rho])
return y, mults
ここでは2つの実装を用意しました。decimate_naive はフルレートで畳み込んでから間引くので、捨てるサンプルの分まで計算してしまいます。decimate_polyphase は係数を $M$ 個のサブフィルタに分け、それぞれを間引き後の低レートで動かすので、無駄な計算が省かれます。実際に動かして乗算回数を比べます。
import numpy as np
from scipy import signal
# 上のセルで定義した decimate_naive / decimate_polyphase を使う想定
np.random.seed(0)
M = 8
x = np.random.randn(8000) # 入力信号
h = signal.firwin(128, 1.0/M) # アンチエイリアスフィルタ
y_naive, m_naive = decimate_naive(x, h, M)
y_poly, m_poly = decimate_polyphase(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:.2e}")
print(f"素朴な実装の乗算回数: {m_naive}")
print(f"ポリフェーズ実装の乗算回数: {m_poly}")
print(f"削減率: {m_poly / m_naive:.3f} (理論値 ~ 1/M = {1/M:.3f})")
実行すると、両実装の出力の最大誤差はほぼゼロ（数値誤差レベル）で、ポリフェーズ実装が素朴な実装と同じ結果を与えることが分かります。一方で乗算回数は、ポリフェーズ実装が素朴な実装のおよそ $1/M$ になっています。これは理論で導いた「ポリフェーズ分解は捨てるサンプルの計算を排除し、計算量を $1/M$ にする」という結論と一致します。同じ結果をはるかに少ない計算で得られるため、リアルタイム処理や組み込み実装でポリフェーズ構成が標準的に使われる理由がよく分かります。
scipy.signalによる検証
自作実装の正しさを、scipyの公式デシメーション関数とも比べておきましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 8000
M = 5
N = 4000
t = np.arange(N)/fs
x = np.sin(2*np.pi*200*t) + 0.5*np.sin(2*np.pi*700*t)
# scipy の decimate（内部で FIR/IIR アンチエイリアス + 間引き）
y_scipy = signal.decimate(x, M, ftype='fir')
t_dec = np.arange(len(y_scipy)) / (fs/M)
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(t[:200], x[:200], 'b-', alpha=0.5, label='input @ 8 kHz')
ax.plot(t_dec[:200//M], y_scipy[:200//M], 'go-', ms=4, label=f'decimated @ {fs//M} Hz')
ax.set_xlabel('Time [s]'); ax.set_ylabel('Amplitude')
ax.set_title(f'scipy.signal.decimate by M={M}')
ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('scipy_decimate.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、scipy.signal.decimate で間引いた信号（緑の点）が、もとの信号（青の線）の200 Hz成分をきれいに保持していることが分かります。アンチエイリアスフィルタが内蔵されているため、700 Hz成分は新ナイキスト（800 Hz）に近いものの適切に処理され、低周波の波形が忠実に再現されています。実務では、こうした検証済みのライブラリ関数を使うのが安全ですが、その内部で何が起きているかを本記事の理論で理解しておくことが、パラメータ選択やトラブルシューティングの土台になります。
まとめ
本記事では、マルチレート信号処理の基礎であるデシメーションと補間について、スペクトルへの作用を導出し、効率的な実装まで解説しました。
- ダウンサンプリング: $Y(e^{j\omega}) = \frac{1}{M}\sum_{k=0}^{M-1} X(e^{j(\omega-2\pi k)/M})$ により、スペクトルは横に $M$ 倍引き伸ばされ $M$ 個のコピーが重畳する。帯域が $\pi/M$ を超えると折り返し（エイリアシング）が起きる
- デシメーション: 間引きの前にカットオフ $\pi/M$ のアンチエイリアスフィルタをかけ、折り返しを防ぐ
- アップサンプリング: $V(e^{j\omega}) = X(e^{j\omega L})$ により、スペクトルは $1/L$ に圧縮され、$L-1$ 個のイメージが現れる
- インターポレーション: ゼロ挿入の後にカットオフ $\pi/L$・ゲイン $L$ の補間フィルタをかけ、イメージを除去して滑らかな信号を作る
- 有理数倍 $L/M$ 変換: アップ→共通フィルタ（カットオフ $\pi/\max(L,M)$）→ダウンの順に行う。順序はアップが先
- ポリフェーズ分解: $H(z) = \sum_\rho z^{-\rho} E_\rho(z^M)$ と分解し、Noble恒等式でフィルタと間引きを交換すると、計算が低レート側で行われ乗算回数が $1/M$ になる
マルチレート処理は、サブバンド符号化・フィルタバンク・ソフトウェア無線・オーバーサンプリング変換といった応用の土台です。ポリフェーズの考え方は、QMFフィルタバンクやウェーブレット変換、効率的なFFTベースのレート変換へと自然に発展していきます。
次のステップとして、以下の記事も参考にしてください。