信号再構成と補間：sinc補間・線形補間・スプライン補間のPython実装

信号再構成と補間（Whittaker-Shannon 補間公式・sinc 補間・Lanczos・bicubic・スプライン）を Python 実装で比較。scipy.interpolate.interp1d・scipy.interpolate.CubicSpline・scipy.interpolate.PchipInterpolator・scipy.signal.resample・numpy.sinc を用い、サンプリング定理に基づく理想再構成、アップサンプリング/ダウンサンプリング、打ち切り誤差、画像補間（bilinear/bicubic/Lanczos）まで体系化。

22 min read

はじめに

サンプリング定理 では、ナイキスト周波数より高い周波数で標本化すれば、離散信号から連続信号を完全に復元できることを証明しました（ナイキスト条件の導出やエイリアシングの仕組みはそちらの記事に譲ります）。本記事ではこの結果を前提として、実際に「どう復元するか」——**信号再構成（reconstruction）と、より一般的な補間（interpolation）**の手法——に焦点を絞り、sinc補間・線形補間・キュービック補間・スプライン補間を数式・周波数応答・実行結果の三面から比較します。

実務では理想的なsinc補間ではなく、計算量と精度のトレードオフから線形補間やスプライン補間を使う場面が多くあります。それぞれの特性を周波数領域から理解しておくと、用途に応じた手法選択ができるようになります。

信号再構成の数学的基礎

標本化された信号の表現

サンプリング周期 \(T_s\) 、サンプリング周波数を

\[ f_s = 1/T_s \]

とします。この周期で標本化された信号は、ディラックのデルタ関数を用いて次のように表されます。

\[x_s(t) = \sum_{n=-\infty}^{\infty} x[n]\,\delta(t - nT_s) \tag{1}\]

ここで \(x[n] = x(nT_s)\) は離散サンプル値です。

理想的な再構成：sinc補間

サンプリング定理の証明から、ナイキスト条件を満たす標本化信号は、理想ローパスフィルタ（カットオフ \(f_s/2\) ）を通すことで完全に元の連続信号に戻せます。理想ローパスフィルタの時間領域インパルス応答は sinc関数 です。

\[h(t) = \text{sinc}\!\left(\frac{t}{T_s}\right) = \frac{\sin(\pi t / T_s)}{\pi t / T_s} \tag{2}\]

これと標本化信号の畳み込みにより、再構成された連続信号は以下で表されます。

\[x(t) = \sum_{n=-\infty}^{\infty} x[n]\,\text{sinc}\!\left(\frac{t - nT_s}{T_s}\right) \tag{3}\]

これがWhittaker–Shannonの補間公式であり、ナイキスト条件下では完全な再構成を与えます。

補間公式が「補間」であることの証明

式 \((3)\) が単なる「近似」ではなく厳密な意味での補間（サンプル点上で元の値を正確に再現する）であることは、sinc関数の次の性質から保証されます。

\[ \text{sinc}(k) = \frac{\sin(\pi k)}{\pi k} = \begin{cases} 1 & k = 0 \\ 0 & k \in \mathbb{Z} \setminus \{0\} \end{cases} \]

（\(k=0\) のときはロピタルの定理より極限値 \(1\) 、\(k\) がゼロでない整数のときは \(\sin(\pi k) = 0\) より \(0\) になります。）

この性質を使って式 \((3)\) に \(t = mT_s\) （\(m\) は整数）を代入すると、

\[ x(mT_s) = \sum_{n=-\infty}^{\infty} x[n]\,\text{sinc}(m - n) = x[m] \]

となります。和の中で \(n = m\) の項だけが \(\text{sinc}(0) = 1\) として残り、それ以外のすべての項は \(\text{sinc}(m-n) = 0\) で消えるからです。つまり式 \((3)\) はサンプル点上では入力値をそのまま再現し、サンプル点の間ではサンプリング定理の一意性により定まる帯域制限信号を与えます。これがsinc補間が数学的に「完全な補間」と呼ばれる根拠です。

実用的な補間との関係

式 \((3)\) は無限和であり、各 sinc 関数の裾が長く、計算コストも遅延も大きいため、リアルタイム処理ではそのままは使えません。実用上は次の2点で近似します。

  • sinc関数を有限長で打ち切る（窓関数を掛ける）
  • sinc以外の補間関数（線形、キュービック、スプライン）で近似する

それぞれの近似がどの程度元信号を保存するかは、その補間関数の周波数応答で決まります。

補間手法の比較

線形補間（Linear Interpolation）

連続する2サンプル \((t_n, x[n])\) と \((t_{n+1}, x[n+1])\) を直線で結びます。

\[x(t) = x[n] + \frac{x[n+1] - x[n]}{T_s}(t - nT_s), \quad nT_s \leq t < (n+1)T_s \tag{4}\]

三角カーネルとsinc²ロールオフの証明

式 \((4)\) はカーネルによる畳み込みとして書き直せます。三角関数（triangle function）を

\[ \Lambda(u) = \begin{cases} 1 - |u| & |u| \leq 1 \\ 0 & |u| > 1 \end{cases} \]

と定義すると、線形補間は次の畳み込み和に一致します。

\[ x_{\text{lin}}(t) = \sum_{n=-\infty}^{\infty} x[n]\, \Lambda\!\left(\frac{t - nT_s}{T_s}\right) \]

（\(t\) が \(nT_s\) と \((n+1)T_s\) の間にあるとき、この和は \(n\) と \(n+1\) の2項だけが非ゼロとなり、式 \((4)\) の直線の式に一致することが確認できます。）

\(\Lambda(u)\) が \(\text{sinc}^2\) のロールオフを持つことは、\(\Lambda\) が幅1の矩形関数 \(\text{rect}(u)\) （\(|u|\leq 1/2\) で \(1\) 、それ以外で \(0\) ）どうしの自己畳み込みであることから証明できます。

\[ \Lambda(u) = (\text{rect} * \text{rect})(u) = \int_{-\infty}^{\infty} \text{rect}(\tau)\,\text{rect}(u - \tau)\, d\tau \]

（矩形どうしが重なる区間の長さが \(|u|\) に対して線形に減少するため、畳み込み結果は三角関数になります。）フーリエ変換の畳み込み定理 \(\mathcal{F}\{f * g\} = \mathcal{F}\{f\}\cdot\mathcal{F}\{g\}\) と \(\mathcal{F}\{\text{rect}\}(f) = \text{sinc}(f)\) を使うと、

\[ \mathcal{F}\{\Lambda\}(f) = \text{sinc}(f)^2 \]

が導かれます。これが「線形補間のロールオフは \(\text{sinc}^2\) 」という主張の厳密な証明です。\(\text{sinc}^2(f)\) は \(\text{sinc}(f)\) より高周波での減衰は速い（絶対値を2乗するため）ものの、多項式的な減衰（\(O(1/f^2)\) ）にとどまり、理想LPFの矩形状のスペクトルには遠く及ばないため、高周波の漏れ込みが残ります。

キュービック補間（Cubic Interpolation）

4つのサンプルから3次多項式を構成し、補間値を計算します。代表的な手法として Catmull–Rom spline や cubic convolution があります。Keys (1981) のcubic convolution kernelは次式です。

\[h(t) = \begin{cases} (a+2)|t|^3 - (a+3)|t|^2 + 1 & |t| \leq 1 \\ a|t|^3 - 5a|t|^2 + 8a|t| - 4a & 1 < |t| < 2 \\ 0 & |t| \geq 2 \end{cases} \tag{5}\]

通常 \(a = -0.5\) が使われます。この値はどこから来るのでしょうか。Keys (1981) は、補間カーネルが2次関数（放物線）を厳密に再現するという条件から \(a=-0.5\) を導出しました（3次までの多項式ではなく2次関数を正確に再現できる時点で、テイラー展開の3次の項まで一致する高い近似次数が得られます）。この性質を数値的に検証してみます。

import numpy as np


def keys_kernel(t, a):
    """Keys (1981) の cubic convolution kernel。"""
    t = np.abs(t)
    h = np.zeros_like(t)
    m1 = t <= 1
    m2 = (t > 1) & (t < 2)
    h[m1] = (a + 2) * t[m1] ** 3 - (a + 3) * t[m1] ** 2 + 1
    h[m2] = a * t[m2] ** 3 - 5 * a * t[m2] ** 2 + 8 * a * t[m2] - 4 * a
    return h


def cubic_conv_interp(x_of_n, t_query, a, n_min, n_max):
    """整数格子上の関数 x_of_n を cubic convolution で補間する。"""
    out = np.zeros_like(t_query, dtype=float)
    for i, t in enumerate(t_query):
        n0 = int(np.floor(t))
        acc = 0.0
        for k in range(n0 - 1, n0 + 3):
            if n_min <= k <= n_max:
                acc += x_of_n(k) * keys_kernel(np.array([t - k]), a)[0]
        out[i] = acc
    return out


# x(n) = n^2 という2次関数を整数格子上でサンプリングし、非整数点で補間する
x_of_n = lambda n: float(n) ** 2
t_query = np.array([0.5, 1.3, 2.7, -1.4, 3.5])
x_true = t_query**2

for a in [-1.0, -0.75, -0.5, -0.25, 0.0]:
    x_hat = cubic_conv_interp(x_of_n, t_query, a, -5, 5)
    err = np.max(np.abs(x_hat - x_true))
    print(f"a={a:+.2f}: max abs error on quadratic = {err:.3e}")

実行結果:

a=-1.00: max abs error on quadratic = 6.300e-01
a=-0.75: max abs error on quadratic = 3.150e-01
a=-0.50: max abs error on quadratic = 3.109e-15
a=-0.25: max abs error on quadratic = 3.150e-01
a=+0.00: max abs error on quadratic = 6.300e-01

\(a=-0.5\) のときだけ誤差が浮動小数点の丸め誤差レベル（\(3.1 \times 10^{-15}\) ）まで下がり、2次関数を厳密に再現できていることが数値的に確認できます。それ以外の \(a\) では、\(a=-0.5\) からのずれに比例して系統的な誤差が残ります（\(a=-1.0\) と \(a=0.0\) で同じ \(0.63\) 、\(a=-0.75\) と \(a=-0.25\) で同じ \(0.315\) という対称性も見て取れます）。これが cubic convolution で \(a=-0.5\) が事実上の標準値として使われる理由です。

スプライン補間（Spline Interpolation）

各区間で滑らかな多項式（通常3次）を当てはめ、隣接する区間の境界で関数値・1次・2次微分が連続するように制約します。3次スプラインは時間領域で \(C^2\) 連続な滑らかな曲線を与え、画像処理やオーディオのアップサンプリングで広く使われます。後述の実行結果で見るように、scipy.interpolate.interp1d(kind="cubic") は内部的に3次スプライン補間を実装しているため、scipy.interpolate.CubicSpline とほぼ同一の結果を返します（既定の端点条件がわずかに異なるため厳密には別物ですが、内部領域では数値的に一致します）。

周波数領域での比較

補間手法カーネル計算量通過帯域平坦性阻止帯域減衰
最近傍矩形最小最悪最悪
線形三角小中中（sinc²）
キュービックKeys (\(a=-0.5\) ) など中良良
スプラインB-spline（3次）中非常に良非常に良
sinc（理想）sinc最大完全完全

選択指針：オーディオは品質優先で sinc または高次スプライン、画像はキュービック、リアルタイム制御は線形が定番です。

補間手法の落とし穴とエッジケース

補間手法にはそれぞれ、実務で見過ごされがちな限界があります。

Runge現象（高次多項式補間の暴れ）：単一の多項式で全区間を近似しようとする大域的高次多項式補間（全サンプル点を通る次数 \(N-1\) のラグランジュ補間多項式など）は、サンプル数を増やすと区間の両端で振動が発散するRunge現象を起こします。本記事で扱う3次スプラインやキュービック補間は、区間ごとに低次多項式を当てはめる区分的（piecewise）多項式であるため、この問題を回避しています。精度を上げたい場合は多項式の次数を上げるのではなく、区分（ノット）の数を増やすのが定石です。

Gibbs現象（sinc補間打ち切りのリンギング）：既に述べた通り、sinc補間の無限和を有限長で打ち切ると、打ち切り境界付近で誤差が拡大します。これは、不連続点近傍でのフーリエ級数の打ち切りが生むGibbs現象と本質的に同じ機構（矩形窓とのスペクトル畳み込み）です。後述の数値実験で、打ち切り誤差が信号の両端に集中し、中央部の誤差よりも桁違いに大きくなることを実測します。

外挿（extrapolation）の危険性：scipy.interpolate.interp1d の fill_value="extrapolate" オプションは、サンプル範囲外の点に対してカーネルの多項式をそのまま延長します。線形補間の外挿は直線の延長なので比較的安全ですが、キュービック補間やスプラインの外挿は3次多項式の急激な発散を招きやすくなります。実際、後掲の再構成グラフでも \(t > 38\) ms の範囲でLinear/Cubicが真の信号から急速に乖離しています。外挿は原理的に「情報のない領域を予測する」行為であり、補間（内挿）とは異なるリスクを伴うことを常に意識する必要があります。

非一様サンプリングでは古典手法が破綻する：本記事で扱った sinc・線形・キュービック・スプラインの各手法は、すべてサンプル間隔が一定であることを前提とした定式化です（式 \((3)(4)(5)\) のいずれも \(nT_s\) という等間隔グリッドを仮定しています）。センサーの欠測やイベント駆動サンプリングなど、サンプル間隔が不均一な場合はこれらの式をそのまま適用できません。Pakiyarajah, Pavez, & Ortega (2024) “Irregularity-Aware Bandlimited Approximation for Graph Signal Interpolation”（ICASSP 2024）は、グラフ信号処理の文脈でノード配置の不規則性を考慮した帯域制限近似を提案しており、均一グリッドを前提とする古典的補間理論を非正則な配置へ拡張する方向性の一例です。

既にエイリアシングした信号は補間で復元不可能：補間はあくまで「サンプル値から連続信号を再構成する」操作です。サンプリング時点で既にナイキスト条件が破られエイリアシングが発生していた場合、失われた情報をどの補間手法でも復元することはできません（詳細は サンプリング定理 を参照）。補間手法の選択は、ナイキスト条件を満たした上での「復元の質」を左右するものであり、サンプリング自体の誤りを補うものではない点に注意してください。

アップサンプリングとダウンサンプリング

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

\(L\) 倍にアップサンプリングする標準的な手順は次の2ステップです。

  1. ゼロ詰め（zero stuffing）：各サンプル間に \(L-1\) 個のゼロを挿入し、サンプリングレートを \(L f_s\) に上げる
  2. アンチイメージングフィルタ：カットオフ \(f_s/2\) のローパスフィルタで、ゼロ詰めで生じたイメージ成分を除去

ゼロ詰め後にローパスフィルタを通すことで、結果として時間領域でのsinc補間と等価な操作になります。

ダウンサンプリング（間引き）

\(M\) 倍にダウンサンプリングする際は、エイリアシングを防ぐため事前にアンチエイリアスフィルタを通す必要があります。

  1. アンチエイリアスフィルタ：カットオフ \(f_s/(2M)\) のローパスフィルタを掛ける
  2. 間引き（decimation）：\(M\) サンプルごとに1点を残し、それ以外を削除

順序を逆にすると、信号の高域成分がエイリアシングして低域に折り返し、後処理では取り除けません。

Python実装：複数補間手法の比較

scipy.interpolate と scipy.signal を使い、線形・キュービック・スプライン・sincの各補間を比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, CubicSpline


def sinc_interp(x_samples, t_samples, t_query):
    """Whittaker–Shannonの補間公式（有限長近似）。"""
    Ts = t_samples[1] - t_samples[0]
    # t_query を列、t_samples を行とする行列を作る
    T_grid, N_grid = np.meshgrid(t_query, t_samples, indexing="ij")
    sinc_matrix = np.sinc((T_grid - N_grid) / Ts)
    return sinc_matrix @ x_samples


# --- 1. 元信号（高解像度） ---
fs_high = 2000          # 「真の」連続信号の近似
T = 0.04                # 40 ms
t_dense = np.arange(0, T, 1 / fs_high)
x_true = (np.sin(2 * np.pi * 50 * t_dense)
          + 0.5 * np.sin(2 * np.pi * 180 * t_dense))

# --- 2. 標本化（fs = 500 Hz、ナイキスト周波数 250 Hz） ---
fs = 500
t_samples = np.arange(0, T, 1 / fs)
x_samples = (np.sin(2 * np.pi * 50 * t_samples)
             + 0.5 * np.sin(2 * np.pi * 180 * t_samples))

# --- 3. 補間 ---
linear = interp1d(t_samples, x_samples, kind="linear", fill_value="extrapolate")
cubic = interp1d(t_samples, x_samples, kind="cubic", fill_value="extrapolate")
spline = CubicSpline(t_samples, x_samples)

x_linear = linear(t_dense)
x_cubic = cubic(t_dense)
x_spline = spline(t_dense)
x_sinc = sinc_interp(x_samples, t_samples, t_dense)

# --- 4. 誤差評価（RMSE） ---
methods = {
    "Linear": x_linear,
    "Cubic": x_cubic,
    "Spline": x_spline,
    "Sinc": x_sinc,
}
for name, x_hat in methods.items():
    rmse = np.sqrt(np.mean((x_true - x_hat) ** 2))
    print(f"{name:8s}: RMSE = {rmse:.4f}")

# --- 5. 可視化 ---
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(t_dense * 1000, x_true, "k-", alpha=0.4, label="True signal")
ax.plot(t_samples * 1000, x_samples, "ko", label="Samples")
ax.plot(t_dense * 1000, x_linear, "--", label="Linear")
ax.plot(t_dense * 1000, x_cubic, "--", label="Cubic")
ax.plot(t_dense * 1000, x_sinc, "-", label="Sinc")
ax.set_xlabel("Time [ms]")
ax.set_ylabel("Amplitude")
ax.set_title("Signal Reconstruction with Different Interpolation Methods")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

実行結果:

Linear  : RMSE = 0.2283
Cubic   : RMSE = 0.0634
Spline  : RMSE = 0.0634
Sinc    : RMSE = 0.0607

サンプル数20点（\(f_s = 500\) Hz、区間40 ms）に対して、Linear のRMSEが最も大きく（\(0.2283\) ）、Sinc が最小（\(0.0607\) ）という理論通りの序列が確認できます。注目すべきは Cubic と Spline のRMSEが小数第4位まで完全に一致している点です。これは前節で述べた通り、interp1d(kind="cubic") と CubicSpline がどちらも内部的に3次スプライン補間を実装しているため、境界条件が影響しない範囲ではほぼ同一の結果を返すという事実の実測的な裏付けです。

sinc・線形・キュービック補間による信号再構成の比較。サンプル点(黒丸)から真の信号(灰色)をどれだけ正確に復元できているかを示す

図から、Linear（青破線）がサンプル間で直線的にカクカクと近似しているのに対し、Cubic（緑破線）とSinc（紫実線）は真の信号（灰色）にほぼ重なっていることが視覚的に確認できます。また \(t > 38\) ms の右端では、範囲外への外挿によってLinearとCubicが真の信号から急速に乖離しており、外挿の危険性が図からも読み取れます。

ナイキスト周波数（250 Hz）以下の成分しか含まないようにすれば sinc 補間は理論的に完全再構成を与えますが、有限長で打ち切るためエッジ近傍では誤差が増えます。実用的には信号の中央部に対して sinc を、両端に対して線形/スプラインを使うハイブリッドアプローチが取られます。

補間誤差のオーバーサンプリング比依存性

前節のRMSEは単一のサンプリング周波数（500 Hz）での比較でした。ここでは観測窓の長さを40 msに固定したまま、サンプリング周波数 \(f_s\) を変化させたときに各補間手法のRMSEがどう変化するかを調べます。オーバーサンプリング比を

\[ r = f_s / 2f_{\max} \]

（\(f_{\max}=180\) Hz）と定義します。

import numpy as np
from scipy.interpolate import interp1d, CubicSpline


def sinc_interp(x_samples, t_samples, t_query):
    Ts = t_samples[1] - t_samples[0]
    T_grid, N_grid = np.meshgrid(t_query, t_samples, indexing="ij")
    sinc_matrix = np.sinc((T_grid - N_grid) / Ts)
    return sinc_matrix @ x_samples


T = 0.04
fs_high = 8000  # 真の信号を高精度に近似するための基準レート
t_dense = np.arange(0, T, 1 / fs_high)
x_true = (np.sin(2 * np.pi * 50 * t_dense)
          + 0.5 * np.sin(2 * np.pi * 180 * t_dense))

fmax = 180.0
fs_list = np.array([380, 450, 500, 600, 700, 800, 1000, 1200, 1500, 2000, 3000])
results = {name: [] for name in ["Linear", "Cubic", "Spline", "Sinc"]}

for fs in fs_list:
    t_samples = np.arange(0, T, 1 / fs)
    x_samples = (np.sin(2 * np.pi * 50 * t_samples)
                 + 0.5 * np.sin(2 * np.pi * 180 * t_samples))

    linear = interp1d(t_samples, x_samples, kind="linear", fill_value="extrapolate")
    cubic = interp1d(t_samples, x_samples, kind="cubic", fill_value="extrapolate")
    spline = CubicSpline(t_samples, x_samples)
    x_sinc = sinc_interp(x_samples, t_samples, t_dense)

    for name, x_hat in zip(["Linear", "Cubic", "Spline", "Sinc"],
                           [linear(t_dense), cubic(t_dense), spline(t_dense), x_sinc]):
        rmse = np.sqrt(np.mean((x_true - x_hat) ** 2))
        results[name].append(rmse)

r = fs_list / (2 * fmax)
for name, vals in results.items():
    print(name, [f"{v:.5f}" for v in vals])
print("oversampling ratio r:", [f"{x:.2f}" for x in r])

実行結果:

Linear ['0.21605', '0.31903', '0.25184', '0.15235', '0.09654', '0.06678', '0.04028', '0.02856', '0.01899', '0.01093', '0.00483']
Cubic ['0.16042', '0.41156', '0.07786', '0.23704', '0.16533', '0.09354', '0.02657', '0.00736', '0.00086', '0.00030', '0.00009']
Spline ['0.16042', '0.41156', '0.07786', '0.23704', '0.16533', '0.09354', '0.02657', '0.00736', '0.00086', '0.00030', '0.00009']
Sinc ['0.06711', '0.07921', '0.06810', '0.05617', '0.04945', '0.04489', '0.03875', '0.03455', '0.03000', '0.02474', '0.01802']
oversampling ratio r: ['1.06', '1.25', '1.39', '1.67', '1.94', '2.22', '2.78', '3.33', '4.17', '5.56', '8.33']

補間誤差(RMSE、対数軸)とオーバーサンプリング比 r=fs/2fmaxの関係。CubicとSplineは完全に重なっており、rが大きくなるほどSincより急速に誤差が減少する

この結果には、素朴な直感（「sincは理想的だから常に最良」）を裏切る興味深い事実が現れています。オーバーサンプリング比 \(r\) が大きくなるほど、Cubic/SplineのRMSEはSincより急速に減少し、\(r=8.33\) では Cubic/Spline のRMSE（\(0.00009\) ）がSinc（\(0.01802\) ）よりも2桁小さくなります。この逆転が起きる理由を、誤差の空間分布を調べて検証します。

import numpy as np
from scipy.interpolate import interp1d


def sinc_interp(x_samples, t_samples, t_query):
    Ts = t_samples[1] - t_samples[0]
    T_grid, N_grid = np.meshgrid(t_query, t_samples, indexing="ij")
    sinc_matrix = np.sinc((T_grid - N_grid) / Ts)
    return sinc_matrix @ x_samples


T = 0.04
fs_high = 8000
t_dense = np.arange(0, T, 1 / fs_high)
x_true = (np.sin(2 * np.pi * 50 * t_dense)
          + 0.5 * np.sin(2 * np.pi * 180 * t_dense))

fs = 3000  # 高オーバーサンプリング比 (r=8.33)
t_samples = np.arange(0, T, 1 / fs)
x_samples = (np.sin(2 * np.pi * 50 * t_samples)
             + 0.5 * np.sin(2 * np.pi * 180 * t_samples))

x_sinc = sinc_interp(x_samples, t_samples, t_dense)
x_cubic = interp1d(t_samples, x_samples, kind="cubic", fill_value="extrapolate")(t_dense)

err_sinc = np.abs(x_true - x_sinc)
err_cubic = np.abs(x_true - x_cubic)

n = len(t_dense)
lo, hi = int(n * 0.2), int(n * 0.8)  # 中央60%の区間

print(f"Sinc  RMSE (全体):   {np.sqrt(np.mean(err_sinc ** 2)):.5f}")
print(f"Sinc  RMSE (中央60%): {np.sqrt(np.mean(err_sinc[lo:hi] ** 2)):.5f}")
print(f"Cubic RMSE (全体):   {np.sqrt(np.mean(err_cubic ** 2)):.5f}")
print(f"Cubic RMSE (中央60%): {np.sqrt(np.mean(err_cubic[lo:hi] ** 2)):.5f}")
print(f"Sinc  最大誤差(端10%): {max(err_sinc[:n // 10].max(), err_sinc[-n // 10:].max()):.5f}")

実行結果:

Sinc  RMSE (全体):   0.01802
Sinc  RMSE (中央60%): 0.00072
Cubic RMSE (全体):   0.00009
Cubic RMSE (中央60%): 0.00001
Sinc  最大誤差(端10%): 0.29210

Sincの誤差は中央60%区間では \(0.00072\) まで下がり、Cubicの中央区間誤差（\(0.00001\) ）と比べても大きな差ではありません。しかし全体RMSEは \(0.01802\) と、中央区間の25倍に膨らんでいます。原因は端部10%区間の最大誤差が \(0.29\) にも達していることです。これは本記事の実装（sinc_interp）が観測窓（40 ms）の外側のサンプルを一切使わずに打ち切っているため、窓の端では sinc の裾が窓外に切り落とされ、Gibbs現象と同じ機構でリンギング誤差が生じるからです。観測窓の長さは \(f_s\) を上げても変わらないため、この打ち切り誤差は\(f_s\) を上げても縮小しません。一方 Cubic/Spline の誤差は「サンプル間の局所的な曲率をどれだけ細かく近似できるか」で決まり、サンプル間隔が縮む（\(f_s\) が上がる）ほど誤差が減り続けます。つまり「打ち切り誤差が支配的な有限窓では、理論上完全なsinc補間が、局所的な区分的多項式補間に負けることがある」というのが、この実験が示す実務上の重要な教訓です。

アップサンプリングの実装：scipy.signal.resample

scipy.signal.resample は内部でFFTを用い、周波数領域でのゼロパディング＋逆FFTにより等価的にsinc補間を行います。

from scipy.signal import resample, resample_poly

# 元信号: fs = 500 Hz
# L=4 倍にアップサンプリング → fs = 2000 Hz
L = 4
x_upsampled_fft = resample(x_samples, len(x_samples) * L)

# resample_poly はマルチレートFIRフィルタを使い、長い信号でも高速
x_upsampled_poly = resample_poly(x_samples, up=L, down=1)

print(f"Original: {len(x_samples)} samples")
print(f"Upsampled: {len(x_upsampled_fft)} samples")

実行結果:

Original: 20 samples
Upsampled: 80 samples

resampleはFFTベースで信号長が変わると遅くなりますが、ナイキスト条件を満たす信号には高品質です。resample_polyはマルチレート信号処理で標準的なポリフェーズ実装を行い、長い時系列でも高速に動作します。ポリフェーズ分解の内部実装や、resample/resample_poly/decimateの詳しい使い分けは マルチレート信号処理の理論とPython実装 で実測ベンチマークとともに扱っているので、そちらを参照してください。本記事では「補間はゼロ詰め＋アンチイメージングフィルタと等価」という原理面に絞ります。

ダウンサンプリングの実装：エイリアシング防止

from scipy.signal import butter, sosfiltfilt, decimate

# fs = 1000 Hz の信号を fs = 250 Hz に間引く（M=4）
fs_orig = 1000
M = 4
fs_new = fs_orig // M
t = np.arange(0, 1, 1 / fs_orig)
# 50 Hz と 300 Hz を含む信号（300 Hz は新しいナイキスト周波数 125 Hz を超える）
x = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 300 * t)

# --- 方法 1: アンチエイリアスフィルタ + 間引き（手動） ---
sos = butter(8, fs_new / 2 - 5, fs=fs_orig, output="sos")  # カットオフ 120 Hz
x_filtered = sosfiltfilt(sos, x)
x_decimated_manual = x_filtered[::M]

# --- 方法 2: scipy.signal.decimate（推奨） ---
x_decimated = decimate(x, M, ftype="iir")

print(f"Original: {len(x)} samples at {fs_orig} Hz")
print(f"Decimated: {len(x_decimated)} samples at {fs_new} Hz")

実行結果:

Original: 1000 samples at 1000 Hz
Decimated: 250 samples at 250 Hz

decimateはアンチエイリアスフィルタと間引きをまとめて行います。ftype="iir"（既定）は8次のChebyshev I型、ftype="fir"は30タップのHamming窓FIRです。リアルタイム性が重要ならIIR、線形位相が必要ならFIRを選びます。decimateとresample_polyのどちらを使うべきかの判断基準や、大きな間引き比での多段適用のコツは、前述の マルチレート信号処理 で詳しく解説しています。

補間結果のスペクトル比較

各補間手法が高周波成分にどう影響するかを見るには、補間後の信号のスペクトルを確認します。線形補間は \(\text{sinc}^2\) のロールオフによりサンプリング周波数の整数倍付近にイメージ成分が残り、高品質の用途では問題になります。スプラインや sinc 補間ではこれが大幅に抑制されます。

from scipy.fft import rfft, rfftfreq

methods = {"Linear": x_linear, "Cubic": x_cubic, "Sinc": x_sinc}
fig, ax = plt.subplots(figsize=(10, 4))
for name, x_hat in methods.items():
    X = np.abs(rfft(x_hat))
    f = rfftfreq(len(x_hat), 1 / fs_high)
    ax.semilogy(f, X / X.max(), label=name)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("|X(f)| (normalized)")
ax.set_title("Spectrum after Interpolation")
ax.set_xlim(0, fs_high / 2)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

補間手法ごとのスペクトル比較（対数軸）。Linearは200Hz以上で緩やかにしか減衰せずピークが残るが、Sincは急峻に減衰する

図から、Sinc（紫）が200 Hz以上でほぼ単調に、しかも最も急峻に減衰しているのに対し、Linear（青）は600〜900 Hz付近にも目に見えるピークが残っており、\(\text{sinc}^2\) ロールオフの緩やかさ（前述の証明の通り \(O(1/f^2)\) にとどまる）が高周波の漏れ込みとして視覚化されています。Cubic（緑）はLinearとSincの中間の減衰特性を示しています。

最新研究動向

古典的な補間理論（sinc・線形・多項式・スプライン）は、あくまで規則的な時間グリッド上の帯域制限信号を前提としています。近年はこの前提を拡張する研究が進んでいます。

Shabanov et al. (2024) “BANF: Band-limited Neural Fields for Levels of Detail Reconstruction”（CVPR 2024、 プロジェクトページ ）は、3次元形状やシーンを表現する**ニューラルフィールド（暗黙的ニューラル表現）**に対して、古典的な意味でのフーリエ解析やローパスフィルタリングを直接適用できないという問題に取り組んでいます。ニューラルフィールドのアーキテクチャに単純な変更を加えることで、フィールド自体を周波数帯域ごとに分解可能にし、Level-of-Detail（詳細度）表現や規則的グリッド上でのサンプリング（marching cubesによるメッシュ化など）を可能にしました。sinc補間や線形補間の背後にある「カーネルによる畳み込みとしての帯域制限」という考え方を、明示的な信号表現を持たないニューラル表現の世界に拡張する試みとして注目されます。

また、Pakiyarajah, Pavez, & Ortega (2024) “Irregularity-Aware Bandlimited Approximation for Graph Signal Interpolation”（ICASSP 2024）は、前節でも触れた通り、グラフ信号処理の文脈でノード配置の不規則性を考慮した帯域制限近似を提案しています。本記事で扱った等間隔グリッド上の補間理論を、任意のグラフ構造・非一様配置に拡張する方向性の一例です。

まとめ

  • 理想的な信号再構成は sinc 補間（Whittaker–Shannon の補間公式）で与えられ、sinc関数の \(\text{sinc}(k)=\delta[k]\) という性質からサンプル点での値を厳密に再現する（証明済み）
  • 線形補間は三角カーネル（矩形関数どうしの自己畳み込み）であり、そのスペクトルが \(\text{sinc}^2\) になることは畳み込み定理から厳密に導出できる
  • キュービック補間の標準パラメータ \(a=-0.5\) は、2次関数を厳密に再現するという条件から導かれ、数値実験でも誤差が浮動小数点精度まで下がることを確認した
  • 実測RMSE（\(f_s=500\) Hz）は Linear \(0.2283\) > Cubic \(\approx\) Spline \(0.0634\) > Sinc \(0.0607\) の順で、理論的な序列と一致する
  • ただし観測窓が固定長の場合、sinc補間の誤差は打ち切り（Gibbs現象類似）に支配されて頭打ちになり、オーバーサンプリング比を上げると局所補間（Cubic/Spline）の方が最終的に高精度になる（\(r=8.33\) で2桁の差）——「sincは常に最良」という直感には注意が必要
  • Runge現象・Gibbs現象・外挿・非一様サンプリング・既存のエイリアシングは、補間手法を選ぶ際に見落とされがちな落とし穴である
  • アップサンプリングは「ゼロ詰め＋アンチイメージングフィルタ」、ダウンサンプリングは「アンチエイリアスフィルタ＋間引き」が標準手順

実装上のポイントは、補間とフィルタリングは表裏一体であるという点です。どの補間手法も「ある周波数応答を持つフィルタによる畳み込み」と見なせます。

関連記事

参考文献

  • Keys, R. G. (1981). “Cubic Convolution Interpolation for Digital Image Processing.” IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(6), 1153-1160.
  • Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
  • Smith, J. O. (2002). Digital Audio Resampling Home Page. CCRMA, Stanford University.
  • Shabanov, A., Govindarajan, S., Reading, C., Goli, L., Rebain, D., Yi, K. M., & Tagliasacchi, A. (2024). “BANF: Band-limited Neural Fields for Levels of Detail Reconstruction.” CVPR 2024, 20571-20580. Project page
  • Pakiyarajah, D., Pavez, E., & Ortega, A. (2024). “Irregularity-Aware Bandlimited Approximation for Graph Signal Interpolation.” ICASSP 2024. arXiv:2312.09405
  • SciPy Interpolation documentation
  • SciPy Signal Resampling documentation

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

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

関連記事