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

はじめに

サンプリング定理では、ナイキスト周波数より高い周波数で標本化すれば、離散信号から連続信号を完全に復元できることを示しました。本記事ではその「復元」、すなわち**信号再構成（reconstruction）と、より一般的な補間（interpolation）**の手法を、数式とPython実装の両面から解説します。

実務では理想的な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 関数の裾が長く、計算コストも遅延も大きいため、リアルタイム処理ではそのままは使えません。実用上は次の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}\]

これは三角関数（triangle function）と離散信号の畳み込みと等価です。三角関数のフーリエ変換は \(\text{sinc}^2\) なので、周波数領域では穏やかな低域通過特性を持ちますが、サイドローブの抑制は不十分で高周波の漏れ込みが起こります。

キュービック補間（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\) が使われ、線形補間より滑らかでsincより計算が軽いという中庸の選択肢になります。

スプライン補間（Spline Interpolation）

各区間で滑らかな多項式（通常3次）を当てはめ、隣接する区間の境界で関数値・1次・2次微分が連続するように制約します。3次スプラインは時間領域で \(C^2\) 連続な滑らかな曲線を与え、画像処理やオーディオのアップサンプリングで広く使われます。

周波数領域での比較

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

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

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

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

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

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

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

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

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

アンチエイリアスフィルタ：カットオフ \(f_s/(2M)\) のローパスフィルタを掛ける
間引き（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()

ナイキスト周波数（250 Hz）以下の成分しか含まないようにすれば sinc 補間は理論的に完全再構成を与えますが、有限長で打ち切るためエッジ近傍では誤差が増えます。実用的には信号の中央部に対して 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")

resampleはFFTベースで信号長が変わると遅くなりますが、ナイキスト条件を満たす信号には高品質です。resample_polyはマルチレート信号処理で標準的なポリフェーズ実装を行い、長い時系列でも高速に動作します。

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

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")

decimateはアンチエイリアスフィルタと間引きをまとめて行います。ftype="iir"（既定）は8次のChebyshev I型、ftype="fir"は30タップのHamming窓FIRです。リアルタイム性が重要ならIIR、線形位相が必要ならFIRを選びます。

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

各補間手法が高周波成分にどう影響するかを見るには、補間後の信号のスペクトルを確認します。線形補間は \(\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()

まとめ

理想的な信号再構成は sinc 補間（Whittaker–Shannon の補間公式）で与えられるが、無限長和なので実装には近似が必要
線形・キュービック・スプライン補間は計算量と精度のトレードオフを取った実用的な代替手段
周波数応答の観点では、線形補間は \(\text{sinc}^2\) のロールオフを持ち、高周波の漏れ込みが避けられない
アップサンプリングは「ゼロ詰め＋アンチイメージングフィルタ」、ダウンサンプリングは「アンチエイリアスフィルタ＋間引き」が標準手順
scipy.signal.resample は FFT による等価sinc補間、resample_poly はポリフェーズ実装、decimate はアンチエイリアシング込みの間引き

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

サンプリング定理：標本化と再構成の数学的基礎 - 本記事の前提となる、ナイキスト条件と完全再構成の数学的背景を解説しています。
Z変換とデジタルフィルタの伝達関数 - 補間フィルタの伝達関数を理解するためのZ変換と極零点解析を解説しています。
高速フーリエ変換（FFT）の仕組みとPython実装 - 補間結果のスペクトル評価で使うFFTの仕組みを解説しています。
窓関数とパワースペクトル密度（PSD） - sinc補間を有限長で打ち切る際に使う窓関数の特性を解説しています。
FIRフィルタとIIRフィルタの比較 - アンチエイリアス・アンチイメージングに使う各種フィルタの選び方を解説しています。
移動平均フィルタの種類と比較 - 簡易的なローパスフィルタとしての移動平均の特性を解説しています。

参考文献

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.
SciPy Interpolation documentation
SciPy Signal Resampling documentation