Signal Reconstruction and Interpolation: Sinc, Linear, and Spline in Python

Introduction

The Sampling Theorem showed that a continuous signal can be perfectly recovered from its samples if it is sampled above the Nyquist rate. This article covers how to perform that recovery — signal reconstruction — and the more general topic of interpolation for resampling and upscaling.

In practice we rarely use ideal sinc interpolation. The trade-offs between linear, cubic, and spline interpolation matter, and understanding each method in the frequency domain lets you pick the right tool for the job.

Mathematical Foundation

The Sampled Signal

A signal sampled at period \(T_s\) (rate \(f_s = 1/T_s\) ) is represented using Dirac deltas:

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

where \(x[n] = x(nT_s)\) .

Ideal Reconstruction: Sinc Interpolation

The sampling theorem proof shows that an ideal lowpass filter with cutoff \(f_s/2\) recovers the original signal. Its impulse response is the sinc function:

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

Convolving the sampled signal with this sinc gives the Whittaker–Shannon interpolation formula:

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

This is perfect reconstruction under the Nyquist condition.

Why Practical Interpolators Differ

Equation (3) is an infinite sum with long sinc tails — too costly for real-time use. We approximate by:

truncating the sinc to finite length (windowing it), or
replacing sinc with a different kernel (linear, cubic, spline).

Each replacement has its own frequency response that determines how well the original signal is preserved.

Comparison of Interpolation Methods

Linear Interpolation

Connect adjacent samples \((t_n, x[n])\) and \((t_{n+1}, x[n+1])\) with a straight line:

\[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}\]

Equivalent to convolving with a triangular kernel. Its Fourier transform is \(\text{sinc}^2\) , so attenuation in the stopband is mediocre — high-frequency leakage is unavoidable.

Cubic Interpolation

Fit a 3rd-order polynomial through four samples. A common choice is Keys’ cubic convolution kernel (1981):

\[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}\]

Typically \(a = -0.5\) . A good middle ground between linear and sinc.

Spline Interpolation

Fit smooth piecewise polynomials (usually cubic) so that function value and first/second derivatives match at each knot. Cubic B-splines yield \(C^2\) -continuous curves and excel for image and audio upsampling.

Frequency-Domain Comparison

Method	Kernel	Cost	Passband Flatness	Stopband Attenuation
Nearest	Rectangle	min	worst	worst
Linear	Triangle	small	medium	medium (\(\text{sinc}^2\) )
Cubic	Keys (\(a=-0.5\) )	mid	good	good
Spline	Cubic B-spline	mid	very good	very good
Sinc (ideal)	sinc	max	perfect	perfect

Practical guidance: audio prefers sinc or high-order spline, images prefer cubic, real-time control loops prefer linear.

Upsampling and Downsampling

Upsampling by \(L\)

Zero-stuffing: insert \(L-1\) zeros between samples → new rate \(L f_s\)
Anti-imaging filter: lowpass with cutoff \(f_s/2\) to remove spectral images

Together these are equivalent to sinc interpolation.

Downsampling by \(M\)

Anti-aliasing filter: lowpass with cutoff \(f_s/(2M)\)
Decimation: keep one out of every \(M\) samples

The order matters — filtering before decimation prevents aliasing that cannot be undone afterward.

Python: Comparing Interpolation Methods

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 interpolation (finite-length truncation)."""
    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


# --- High-resolution "true" continuous signal ---
fs_high = 2000
T = 0.04
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))

# --- Sample at 500 Hz (Nyquist = 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))

# --- Interpolate ---
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)

# --- RMSE evaluation ---
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}")

# --- Visualize ---
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()

When the signal contains only frequencies below the Nyquist rate, sinc interpolation gives near-perfect reconstruction — but truncation produces edge errors. A common trick is to use sinc in the bulk of the signal and a smoother method (linear or spline) near the edges.

Upsampling: scipy.signal.resample

scipy.signal.resample performs FFT-based zero padding in the frequency domain — essentially sinc interpolation.

from scipy.signal import resample, resample_poly

# Upsample by L = 4: 500 Hz → 2000 Hz
L = 4
x_upsampled_fft = resample(x_samples, len(x_samples) * L)

# resample_poly: polyphase FIR — much faster for long signals
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 is high quality but slow on long signals; resample_poly uses a polyphase implementation and stays fast.

Downsampling: Preventing Aliasing

from scipy.signal import butter, sosfiltfilt, decimate

fs_orig = 1000
M = 4
fs_new = fs_orig // M
t = np.arange(0, 1, 1 / fs_orig)
x = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 300 * t)
# 300 Hz exceeds the new Nyquist (125 Hz) — aliasing must be prevented.

# --- Manual: anti-alias filter + decimate ---
sos = butter(8, fs_new / 2 - 5, fs=fs_orig, output="sos")
x_filtered = sosfiltfilt(sos, x)
x_decimated_manual = x_filtered[::M]

# --- scipy.signal.decimate (recommended) ---
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 does both steps. ftype="iir" (default) uses an 8th-order Chebyshev I; ftype="fir" uses a 30-tap Hamming-window FIR. Pick IIR for speed, FIR for linear phase.

Spectrum Comparison

To see how each method affects high frequencies, plot the spectrum of the interpolated signal. Linear interpolation leaves spectral images near integer multiples of \(f_s\) due to its \(\text{sinc}^2\) rolloff; spline and sinc suppress them effectively.

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

Summary

The ideal reconstruction is sinc interpolation (Whittaker–Shannon), but it is an infinite sum and must be approximated.
Linear / cubic / spline are practical approximations with different cost-accuracy trade-offs.
In the frequency domain, linear interpolation has a \(\text{sinc}^2\) rolloff and cannot suppress all high-frequency leakage.
Upsampling = zero-stuff + anti-imaging filter. Downsampling = anti-alias filter + decimate. Order matters.
scipy.signal.resample is FFT-based; resample_poly is polyphase; decimate bundles anti-aliasing with decimation.

The key insight is that interpolation is filtering — every interpolator can be viewed as convolution with a kernel having a specific frequency response.

Sampling Theorem — the Nyquist condition and perfect reconstruction proof that this article builds on.
Z-Transform and Digital Filter Design — pole-zero analysis of interpolation filters.
Fast Fourier Transform (FFT) with Python — the FFT used to evaluate interpolation spectra.
Window Functions and Power Spectral Density — windowing for finite-length sinc truncation.
FIR vs IIR Filters — choosing anti-alias / anti-imaging filters.
Moving Average Filters Compared — moving average as a simple lowpass.

References

Keys, R. G. (1981). “Cubic Convolution Interpolation for Digital Image Processing.” IEEE Trans. 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
SciPy Signal Resampling