マルチレート信号処理の理論とPython実装：デシメーション・補間・ポリフェーズフィルタ

scipy.signal.decimate / resample_poly / upfirdn で学ぶマルチレート信号処理の理論・設計・Python実装。デシメーション（間引き）と補間の周波数領域解析、有理数比 L/M リサンプリング、ポリフェーズ分解による計算量削減、Noble identities まで体系的に解説します。

7 min read

はじめに

サンプリング定理により、帯域制限された信号はナイキスト周波数以上のサンプリング周波数で標本化すれば完全に復元できます。しかし実際のシステムでは、信号の途中でサンプリング周波数そのものを変えたい場面が頻繁に発生します。

  • オーディオ規格間の変換：CD（44.1 kHz）↔ プロ音響・動画（48 kHz）
  • ΔΣ型ADC：数MHzでオーバーサンプリングした1bit信号を数十kHzへ間引く
  • 通信システム：シンボルレートとDAC/ADCレートの整合、パルス整形
  • ウェーブレット変換・フィルタバンク：各段で1/2に間引く構造そのもの

このように複数のサンプリングレートが混在するシステムを扱う技術がマルチレート信号処理（multirate signal processing）です。本記事では、ダウンサンプリング（デシメーション）とアップサンプリング（補間）の周波数領域での振る舞い、有理数比リサンプリング、そして計算量を劇的に削減するポリフェーズ分解と Noble identities を、数式とPython実装（scipy.signal.decimate / resample / resample_poly / upfirdn）で体系的に解説します。

離散信号処理の全体像は離散信号処理の基礎ロードマップを、連続時間側の補間理論は信号再構成と補間を参照してください。

ダウンサンプリング（デシメーション）

定義と周波数領域の解析

信号 \(x[n]\) を \(M\) サンプルごとに1つだけ残す操作（\(M\) 倍のダウンサンプリング）は次のように定義されます。

\[y[m] = x[Mm] \tag{1}\]

サンプリング周波数は \(f_s \to f_s / M\) に低下します。このときDTFT領域では、スペクトルが \(M\) 倍に引き伸ばされ、さらに \(M\) 個のシフト複製が重なり合います。

\[Y(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} X\!\left(e^{j(\omega - 2\pi k)/M}\right) \tag{2}\]

\(k = 0\) の項が目的のスペクトル（\(M\) 倍に伸長されたもの）で、\(k \neq 0\) の項はエイリアシング成分です。元の信号が \(|\omega| < \pi / M\) （実周波数で \(f < f_s / 2M\) ）に帯域制限されていれば複製同士は重ならず、エイリアシングは発生しません。

アンチエイリアシングフィルタの必要性

実際の信号は \(\pi / M\) 以上の成分を含むのが普通なので、間引きの前にカットオフ \(\pi / M\) のローパスフィルタを必ず入れます。この「LPF + 間引き」のセットを**デシメーション（decimation）**と呼びます。

\[y[m] = \sum_{k} h[k] \, x[Mm - k] \tag{3}\]

ここで \(h[n]\) はアンチエイリアシングフィルタのインパルス応答です。フィルタ設計の基礎はローパスフィルタの設計と比較とFIRフィルタとIIRフィルタの比較を参照してください。

Python実装：フィルタなし間引きとの比較

scipy.signal.decimate はアンチエイリアシングフィルタと間引きをまとめて実行します。フィルタなしの単純間引き x[::M] と比較して、エイリアシングの発生を確認します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import decimate

# --- テスト信号: 50 Hz の目的成分 + 380 Hz の高周波成分 ---
fs = 1000                      # 元のサンプリング周波数 [Hz]
t = np.arange(0, 2.0, 1/fs)
x = np.sin(2*np.pi*50*t) + 0.8*np.sin(2*np.pi*380*t)

M = 4                          # 間引き率: 1000 Hz → 250 Hz（新ナイキスト 125 Hz）
fs_new = fs / M

# (1) フィルタなしの単純間引き → 380 Hz が 120 Hz に折り返す
x_naive = x[::M]

# (2) decimate: アンチエイリアシングフィルタ + 間引き
x_dec = decimate(x, M, ftype='fir', zero_phase=True)

def spectrum_db(sig, fs):
    """Hann窓付き振幅スペクトル [dB]（最大値で正規化）"""
    win = np.hanning(len(sig))
    X = np.abs(np.fft.rfft(sig * win))
    f = np.fft.rfftfreq(len(sig), 1/fs)
    return f, 20*np.log10(X / X.max() + 1e-12)

fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=False)

f0, S0 = spectrum_db(x, fs)
axes[0].plot(f0, S0)
axes[0].set_title('Original (fs=1000 Hz): peaks at 50 Hz and 380 Hz')
axes[0].set_xlim(0, 500)

f1, S1 = spectrum_db(x_naive, fs_new)
axes[1].plot(f1, S1, color='crimson')
axes[1].set_title('Naive x[::4] (fs=250 Hz): 380 Hz aliases to 120 Hz!')
axes[1].set_xlim(0, 125)

f2, S2 = spectrum_db(x_dec, fs_new)
axes[2].plot(f2, S2, color='seagreen')
axes[2].set_title('scipy.signal.decimate (fs=250 Hz): alias suppressed')
axes[2].set_xlim(0, 125)

for ax in axes:
    ax.set_ylabel('Magnitude [dB]')
    ax.set_ylim(-100, 5)
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

単純間引きでは 380 Hz の成分が新しいサンプリング周波数 \(f_s' = 250\) Hz に対して折り返し、\(380 - 250 = 130\) Hz → さらにナイキスト（125 Hz）で折り返して \(250 - 130 = 120\) Hz に偽ピークとして出現し、目的帯域を汚染します。decimate を使えばこの成分は間引き前にフィルタで除去されます。

なお decimate のデフォルトは8次Chebyshev I型IIRフィルタです。zero_phase=True（デフォルト）では filtfilt 相当の零位相フィルタリングが行われ、位相歪みを避けられます。線形位相が必要でフィルタ特性を明示したい場合は ftype='fir' を推奨します。また、\(M > 13\) 程度の大きな間引き率では、decimate を複数回に分けて段階的に適用するのが定石です（例：\(M = 32 \to 8 \times 4\) ）。

アップサンプリング（補間）

ゼロ挿入とイメージ成分

\(L\) 倍のアップサンプリングは、まずサンプル間に \(L - 1\) 個のゼロを挿入します。

\[x_u[n] = \begin{cases} x[n/L] & (n = 0, \pm L, \pm 2L, \ldots) \\ 0 & (\text{otherwise}) \end{cases} \tag{4}\]

このときスペクトルは

\[X_u(e^{j\omega}) = X(e^{j\omega L}) \tag{5}\]

となり、\(1/L\) に圧縮された元スペクトルの複製（イメージ成分）が \(\omega = 2\pi k / L\) の周りに \(L\) 個並びます。

補間フィルタ

イメージ成分を除去するため、ゼロ挿入の後に利得 \(L\) ・カットオフ \(\pi / L\) のローパスフィルタを適用します。

\[H_I(e^{j\omega}) = \begin{cases} L & (|\omega| < \pi / L) \\ 0 & (\text{otherwise}) \end{cases} \tag{6}\]

利得が \(L\) 必要なのは、ゼロ挿入によって信号のエネルギーが \(1/L\) に薄まるためです。この「ゼロ挿入 + LPF」のセットを**補間（interpolation）**と呼びます。理想フィルタは sinc 補間そのものであり、信号再構成と補間で扱った理論と直結しています。

Python実装：ゼロ挿入・upfirdn・resample_poly

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample, resample_poly, upfirdn, firwin

# --- テスト信号: 500 Hz + 1500 Hz ---
fs = 8000
t = np.arange(0, 0.2, 1/fs)
x = np.sin(2*np.pi*500*t) + 0.5*np.sin(2*np.pi*1500*t)

L = 3                          # 8000 Hz → 24000 Hz
fs_up = fs * L

# (1) ゼロ挿入のみ → イメージ成分が残る
x_zero = np.zeros(len(x) * L)
x_zero[::L] = x

# (2) ゼロ挿入 + 補間フィルタ（利得 L, カットオフ π/L）を upfirdn で
h = L * firwin(121, 1/L)       # firwin のカットオフはナイキスト正規化
x_upfirdn = upfirdn(h, x, up=L)

# (3) ポリフェーズ実装の resample_poly（内部でカイザー窓FIRを自動設計）
x_poly = resample_poly(x, L, 1)

# (4) FFTベースの resample
x_fft = resample(x, len(x) * L)

def spectrum_db(sig, fs):
    win = np.hanning(len(sig))
    X = np.abs(np.fft.rfft(sig * win))
    f = np.fft.rfftfreq(len(sig), 1/fs)
    return f, 20*np.log10(X / X.max() + 1e-12)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))

f0, S0 = spectrum_db(x_zero, fs_up)
axes[0].plot(f0, S0, color='crimson')
axes[0].set_title('Zero insertion only: images at 6.5/8.5/9.5 kHz etc.')

f1, S1 = spectrum_db(x_upfirdn, fs_up)
axes[1].plot(f1, S1, color='seagreen')
axes[1].set_title('upfirdn (zero insertion + LPF): images removed')

f2, S2 = spectrum_db(x_poly, fs_up)
axes[2].plot(f2, S2, color='navy')
axes[2].set_title('resample_poly(x, 3, 1): same result, auto-designed filter')

for ax in axes:
    ax.set_ylabel('Magnitude [dB]')
    ax.set_ylim(-100, 5)
    ax.set_xlim(0, fs_up / 2)
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

ゼロ挿入のみでは 24 kHz レートにおいて \(8000 \pm 500\) , \(8000 \pm 1500\) Hz などにイメージ成分が並びますが、補間フィルタを通すと元の 500 / 1500 Hz 成分だけが残ります。スペクトルの読み方はFFTによるスペクトル解析と窓関数とPSDを参照してください。

有理数比リサンプリング（L/M 変換）

サンプリング周波数を有理数比 \(L / M\) 倍に変換するには、**アップサンプリング（\(\uparrow L\) ）→ ローパスフィルタ → ダウンサンプリング（\(\downarrow M\) ）**の順に接続します。先に間引くと情報が失われるため、必ずアップサンプリングが先です。

中央のフィルタは、補間フィルタ（カットオフ \(\pi / L\) ）とアンチエイリアシングフィルタ（カットオフ \(\pi / M\) ）を1つで兼ねるため、

\[H(e^{j\omega}) = \begin{cases} L & \left(|\omega| < \min\left(\dfrac{\pi}{L}, \dfrac{\pi}{M}\right)\right) \\ 0 & (\text{otherwise}) \end{cases} \tag{7}\]

とします。代表例が CD → 48 kHz 変換で、\(48000 / 44100 = 160 / 147\) です。

import numpy as np
import matplotlib.pyplot as plt
from math import gcd
from scipy.signal import resample_poly

# --- 44.1 kHz → 48 kHz 変換 ---
fs_in, fs_out = 44100, 48000
g = gcd(fs_out, fs_in)
L, M = fs_out // g, fs_in // g
print(f"L/M = {L}/{M}")        # → 160/147

# 1 kHz 正弦波を変換して精度を確認
duration = 0.05
t_in = np.arange(0, duration, 1/fs_in)
x_in = np.sin(2*np.pi*1000*t_in)

x_out = resample_poly(x_in, L, M)
t_out = np.arange(len(x_out)) / fs_out

# 理論値（連続正弦波を 48 kHz で直接サンプリングしたもの）と比較
x_ref = np.sin(2*np.pi*1000*t_out)
edge = 200                     # 端点の過渡部分を除いて評価
err = np.max(np.abs(x_out[edge:-edge] - x_ref[edge:-edge]))
print(f"最大誤差: {err:.2e}")   # → 1e-5 前後（フィルタ精度で決まる）

plt.figure(figsize=(10, 4))
plt.plot(t_in[:100]*1000, x_in[:100], 'o-', label='44.1 kHz input',
         markersize=4, alpha=0.7)
plt.plot(t_out[:109]*1000, x_out[:109], '.-', label='48 kHz output',
         markersize=4, alpha=0.7)
plt.xlabel('Time [ms]')
plt.ylabel('Amplitude')
plt.title('Rational resampling 44.1 kHz -> 48 kHz (L/M = 160/147)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

resample_poly(x, 160, 147) の1行で高品質なレート変換ができます。ただし内部では「160倍に上げてから147で間引く」処理をポリフェーズ構成で効率化しており、素朴に実装した場合の \(1/160\) 程度の計算量で済んでいます。その仕組みを次節で見ていきます。

ポリフェーズ分解

素朴な実装の無駄

式 (3) のデシメーションを素朴に実装すると、「入力レート \(f_s\) で全サンプルにFIRフィルタを適用 → \(M\) 個に1個だけ残す」となり、計算した出力の \(M - 1 / M\) を捨てることになります。タップ数 \(N\) のFIRフィルタなら、毎秒の積和演算数は \(N f_s\) です。

Type-1 ポリフェーズ分解

この無駄を排除するのがポリフェーズ分解です。フィルタ \(H(z) = \sum_n h[n] z^{-n}\) を、係数を \(M\) 個おきに取り出した \(M\) 本の**サブフィルタ（ポリフェーズ成分）**に分解します。

\[H(z) = \sum_{p=0}^{M-1} z^{-p} E_p(z^M) \tag{8}\] \[E_p(z) = \sum_{n} h[nM + p] \, z^{-n} \qquad (p = 0, 1, \ldots, M-1) \tag{9}\]

例えば \(M = 3\) , \(h = [h_0, h_1, h_2, h_3, h_4, h_5]\) なら、\(E_0 = [h_0, h_3]\) , \(E_1 = [h_1, h_4]\) , \(E_2 = [h_2, h_5]\) です。

デシメータのポリフェーズ構成

式 (8) を式 (3) に代入して整理すると、デシメーション出力は

\[y[m] = \sum_{p=0}^{M-1} \sum_{n} e_p[n] \, x_p[m - n], \qquad x_p[m] = x[mM - p] \tag{10}\]

と書けます。つまり「入力を \(M\) 本の低レート系列 \(x_p\) に振り分け（コミュテータ）、各系列を短いサブフィルタ \(E_p\) に通して加算する」構造です。すべてのフィルタ演算が出力レート \(f_s / M\) で動くため、毎秒の積和演算数は

\[C_{\text{direct}} = N f_s \quad \longrightarrow \quad C_{\text{poly}} = M \cdot \frac{N}{M} \cdot \frac{f_s}{M} = \frac{N f_s}{M} \tag{11}\]

と \(1/M\) に削減されます。補間器も同様にポリフェーズ化でき、計算量は \(1/L\) になります（ゼロとの乗算をスキップできるため）。有理数比 \(L/M\) 変換ではこの両方の効果が効きます。

Python実装：ポリフェーズデシメータの検証

式 (10) を直接実装し、素朴な「フィルタ → 間引き」および scipy.signal.upfirdn と一致することを確認します。

import numpy as np
from scipy.signal import firwin, upfirdn
import time

def polyphase_decimate(x, h, M):
    """ポリフェーズ構成によるデシメーション（式(10)の直接実装）"""
    N_out = len(x) // M
    y = np.zeros(N_out)
    for p in range(M):
        e_p = h[p::M]                        # サブフィルタ E_p
        # 低レート入力系列 x_p[m] = x[mM - p]
        idx = np.arange(N_out) * M - p
        x_p = np.zeros(N_out)
        valid = idx >= 0
        x_p[valid] = x[idx[valid]]
        y += np.convolve(x_p, e_p)[:N_out]   # 低レートでの畳み込み
    return y

# --- 検証 ---
rng = np.random.default_rng(42)
fs, M = 48000, 4
x = rng.standard_normal(fs)                  # 1秒分のホワイトノイズ
h = firwin(128, 1/M)                         # カットオフ π/M のFIR

y_naive = np.convolve(x, h)[::M][:len(x)//M] # 素朴: 全レートでフィルタ→間引き
y_poly = polyphase_decimate(x, h, M)
y_scipy = upfirdn(h, x, up=1, down=M)[:len(x)//M]

print("polyphase == naive :", np.allclose(y_poly, y_naive))   # True
print("polyphase == scipy :", np.allclose(y_poly, y_scipy))   # True

# --- 計算量の比較（理論値） ---
N = len(h)
print(f"直接形の積和演算 : {N * fs / 1e6:.1f} M MAC/s")
print(f"ポリフェーズ形   : {N * fs / M / 1e6:.1f} M MAC/s (1/{M})")

3つの実装の出力が一致することが確認できます。scipy.signal.upfirdn・resample_poly はこのポリフェーズ構成をC実装しており、ゼロ挿入した中間信号を実際には生成しません。

Noble identities

ポリフェーズ構成の正当性を支えるのが **Noble identities（ノーブル恒等式）**です。フィルタとレート変換器の順序交換に関する次の2つの等価関係を指します。

ダウンサンプリング側：\(M\) 倍間引きの後に \(H(z)\) を適用するのは、\(H(z^M)\) を適用してから間引くのと等価です。

\[\left(\downarrow M\right) \;\to\; H(z) \quad \equiv \quad H(z^M) \;\to\; \left(\downarrow M\right) \tag{12}\]

アップサンプリング側：\(H(z)\) を適用してから \(L\) 倍ゼロ挿入するのは、ゼロ挿入後に \(H(z^L)\) を適用するのと等価です。

\[H(z) \;\to\; \left(\uparrow L\right) \quad \equiv \quad \left(\uparrow L\right) \;\to\; H(z^L) \tag{13}\]

ここで \(H(z^M)\) は、インパルス応答の各係数の間に \(M - 1\) 個のゼロを挿入したフィルタです。ポリフェーズ分解（式 8）の各項 \(z^{-p} E_p(z^M)\) に式 (12) を適用すると、\(E_p(z^M)\) を間引き器の後ろ（低レート側）へ移動でき、式 (10) の効率的な構造が導かれます。**「フィルタは常に低レート側で動かす」**というマルチレート設計の大原則は、この恒等式によって保証されています。ウェーブレットのフィルタバンクやΔΣ変調器のデシメータ設計でも、この書き換えが多段構成の効率化に使われます。

decimate / resample / resample_poly の使い分け

SciPyの3つのリサンプリング関数の特徴を整理します。

関数方式対応レート比位相特性端点の挙動適した用途
scipy.signal.decimateLPF（IIR/FIR）+ 間引き整数分の1のみzero_phase=True で零位相フィルタ過渡あり整数比のダウンサンプリング専用
scipy.signal.resampleFFTベース（周波数領域で切断/拡張）任意（実数可）零位相周期信号を仮定、非周期端でリンギング周期的・帯域制限された信号、オフライン解析
scipy.signal.resample_polyポリフェーズFIR（↑L → LPF → ↓M）有理数 \(L/M\)線形位相（群遅延補正済み）padtype で制御可能実用の第一選択。ストリーミング・実データ全般
scipy.signal.upfirdnポリフェーズFIR（フィルタは自前）有理数 \(L/M\)使用フィルタに依存ゼロパディングカスタムフィルタでのレート変換、低レベル制御

実務での選択指針は次の通りです。

  • 迷ったら resample_poly：非周期の実データでもアーティファクトが少なく、有理数比に対応
  • resample は信号が周期的か十分に帯域制限されている場合のみ。端点のリンギングに注意
  • decimate は整数比ダウンサンプリング専用の便利関数。大きな \(M\) では多段適用
  • フィルタの阻止域減衰量や遷移帯域を自分で設計したい場合は firwin などで設計して upfirdn へ（設計指針はデジタルフィルタ設計指針ハブを参照）

実用上の応用と設計のポイント

応用分野典型的なレート変換ポイント
オーディオSRC44.1 kHz ↔ 48 kHz（160/147）高精度FIR（阻止域 −100 dB級）をポリフェーズで実装
ΔΣ型ADCのデシメータ数MHz → 数十kHz（\(M\) 大）CICフィルタ + 多段FIRの縦続、Noble identitiesで効率化
通信のパルス整形シンボルレートの4〜8倍へ補間ルートレイズドコサインフィルタをポリフェーズ補間器として実装
ウェーブレット変換各段で \(\downarrow 2\)2チャネルフィルタバンク＝ポリフェーズ行列で表現
スペクトル解析の前処理対象帯域までダウンサンプルFFT点数と周波数分解能の節約

また、バンドパスフィルタと組み合わせて特定帯域だけを低レートに落とす帯域デシメーションは、狭帯域信号の解析で計算量を大きく削減する定番テクニックです。

まとめ

  • ダウンサンプリングはスペクトルを \(M\) 倍に伸長して複製を重ねるため、事前にカットオフ \(\pi/M\) のアンチエイリアシングフィルタが必須（= デシメーション）
  • アップサンプリングはゼロ挿入によりイメージ成分を生むため、利得 \(L\) ・カットオフ \(\pi/L\) の補間フィルタで除去する
  • 有理数比 \(L/M\) 変換は「\(\uparrow L\) → LPF → \(\downarrow M\) 」の順で構成し、フィルタはカットオフ \(\min(\pi/L, \pi/M)\) の1本で兼ねる
  • ポリフェーズ分解はフィルタを \(M\) 本のサブフィルタに分けて低レート側で動かす構造で、計算量を \(1/M\) （補間では \(1/L\) ）に削減する
  • Noble identities はフィルタとレート変換器の順序交換を保証し、「フィルタは低レート側で」という設計原則の理論的基盤となる
  • SciPyでは resample_poly が実用の第一選択。decimate は整数比ダウンサンプリング、resample は周期信号向け、upfirdn はカスタムフィルタ用と使い分ける

関連記事

参考文献

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox Frequency Converter Hz/kHz/MHz/rad/s 単位を相互変換 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション

関連記事