05_04_03 - Interpolation Kernels

Overview

Interpolation kernels reconstruct continuous signal values between discrete samples. This is fundamental for any DSP operation requiring sub-sample precision:

• Fractional delays (chorus, flanger, vibrato)
• Sample rate conversion (resampling)
• Wavetable synthesis (oscillators)
• Pitch shifting (time-domain)
• Modulated delays (tape effects)

Implemented Methods

1. Linear Interpolation (DRAFT Quality)

Function: lerp_kernel(a, b, frac)

// Simplest: straight line between two points
float a = 1.0f, b = 2.0f;
float result = lerp_kernel(a, b, 0.5f);  // → 1.5

Mathematical Definition:

y = a + frac · (b - a)
  = a · (1 - frac) + b · frac

Properties: - ✅ Fastest (2 samples, 3 operations) - ✅ C0 continuous (no value jumps) - ❌ NOT C1 continuous (derivative discontinuities) - ❌ Introduces high-frequency aliasing

When to use: - Real-time LFOs with gentle modulation - Prototyping/quick iteration - CPU-constrained embedded systems

---

2. Cubic Interpolation (LOW Quality)

Function: cubic_interp_kernel(xm1, x0, x1, x2, frac)

// 4-point Lagrange polynomial
float xm1 = 0.0f, x0 = 1.0f, x1 = 2.0f, x2 = 3.0f;
float result = cubic_interp_kernel(xm1, x0, x1, x2, 0.5f);

Mathematical Definition: Uses Lagrange 4-point interpolation to fit a cubic polynomial through samples at positions [n-1, n, n+1, n+2].

Properties: - ✅ C0 continuous - ✅ Much smoother than linear - ✅ Reduced aliasing - ⚠️ Requires 4 samples (boundary issues) - ⚠️ Moderate computational cost

When to use: - Chorus/flanger effects - Wavetable oscillators (basic quality) - General-purpose fractional delay

---

3. Hermite Interpolation (MEDIUM Quality - Recommended)

Function: hermite_interp_kernel(xm1, x0, x1, x2, frac)

// Cubic with derivative constraints
float xm1 = 0.0f, x0 = 1.0f, x1 = 2.0f, x2 = 3.0f;
float result = hermite_interp_kernel(xm1, x0, x1, x2, 0.5f);

Mathematical Definition: Hermite polynomial ensuring C1 continuity by matching both values and derivatives (using centered differences):

m0 = (x1 - xm1) / 2   // slope at x0
m1 = (x2 - x0) / 2    // slope at x1

Properties: - ✅ C1 continuous (smooth derivatives) - ✅ Very musical sound - ✅ Same cost as cubic - ✅ Excellent for wavetable synthesis

When to use: (DEFAULT CHOICE) - Wavetable oscillators (production quality) - Vibrato/chorus with smooth modulation - General-purpose high-quality interpolation

---

4. Windowed-Sinc Interpolation (HIGH/ULTRA Quality)

Function: sinc_interp_kernel(buffer, position, buffer_size, window_size)

// Near-perfect band-limited interpolation
float buffer[64] = { /* samples */ };
float result = sinc_interp_kernel(buffer, 16.5f, 64, 16);  // 16-tap sinc

Mathematical Definition: Convolves signal with windowed sinc function (ideal interpolation kernel):

sinc(x) = sin(π·x) / (π·x)
windowed_sinc(x) = sinc(x) · blackman_window(x)

Properties: - ✅ C∞ continuous - ✅ Near-perfect band-limiting (minimal aliasing) - ✅ Theoretically ideal (Shannon sampling theorem) - ❌ Computationally expensive (window_size × cost) - ⚠️ Requires large window (8-32 samples)

Window Size Recommendations: - window_size = 8: Good quality, moderate cost - window_size = 16: Excellent quality (recommended for mastering) - window_size = 32: Near-perfect (offline processing only)

When to use: - Offline sample rate conversion - Mastering-quality pitch shifting - Academic research / reference implementations

---

Quality vs Performance Trade-off

Method	Samples	Operations	Quality	Aliasing	Use Case
Linear	2	~3	⭐	High	Real-time LFOs, embedded
Cubic	4	~15	⭐⭐⭐	Moderate	Chorus, basic wavetables
Hermite	4	~20	⭐⭐⭐⭐	Low	Default choice
Sinc-8	8	~40	⭐⭐⭐⭐⭐	Very low	High-quality effects
Sinc-16	16	~80	⭐⭐⭐⭐⭐	Minimal	Mastering, offline

---

Usage Examples

Basic Interpolation

#include "interpolation_kernels.h"

using namespace audiolab::kernels::l0;

float samples[] = {1.0f, 2.0f, 3.0f, 4.0f};

// Linear interpolation at position 1.5
float linear = lerp_buffer_kernel(samples, 1.5f, 4);
// → 2.5

// Hermite interpolation (requires 4+ samples)
float hermite = hermite_buffer_kernel(samples, 1.5f, 4);
// → ~2.5 (slightly smoother)

Quality Selection

// Automatic quality selection
float buffer[100];
float position = 45.7f;

// Draft quality (linear)
float draft = interpolate_kernel(buffer, position, 100, InterpQuality::DRAFT);

// Medium quality (hermite - recommended)
float medium = interpolate_kernel(buffer, position, 100, InterpQuality::MEDIUM);

// Ultra quality (sinc-16)
float ultra = interpolate_kernel(buffer, position, 100, InterpQuality::ULTRA);

Fractional Delay (Chorus Effect)

// Simulate chorus with time-varying delay
float delay_buffer[48000];  // 1 second at 48kHz
float lfo_phase = 0.0f;
float output_buffer[1024];

for (size_t i = 0; i < 1024; ++i) {
    // LFO modulates delay time
    float lfo = std::sin(2.0f * M_PI * lfo_phase);
    float delay_samples = 100.0f + 20.0f * lfo;  // 100±20 samples

    // Hermite interpolation for smooth delay
    output_buffer[i] = hermite_buffer_kernel(delay_buffer, delay_samples, 48000);

    lfo_phase += 1.0f / 1024.0f;  // LFO frequency
}

Wavetable Oscillator

// High-quality wavetable playback
float wavetable[2048];  // One cycle of waveform
float phase = 0.0f;
float frequency = 440.0f;
float sample_rate = 48000.0f;

for (size_t i = 0; i < 1024; ++i) {
    // Phase in wavetable
    float table_position = phase * 2048.0f;

    // Hermite interpolation for smooth playback
    float sample = hermite_buffer_kernel(wavetable, table_position, 2048);

    // Advance phase
    phase += frequency / sample_rate;
    if (phase >= 1.0f) phase -= 1.0f;
}

Batch Processing (SIMD-friendly)

// Interpolate multiple positions at once
float buffer[1000];
float positions[128];  // Positions to interpolate
float output[128];

// Fill positions...
for (size_t i = 0; i < 128; ++i) {
    positions[i] = 100.0f + i * 5.5f;
}

// Batch linear interpolation (SIMD-optimizable)
lerp_batch_kernel(buffer, positions, output, 128, 1000);

---

Frequency Response Comparison

Linear Interpolation

• Passband: DC to ~0.4 × Nyquist (flat)
• Stopband: Heavy aliasing above 0.5 × Nyquist
• Transition: Sharp rolloff with severe aliasing

Cubic Interpolation

• Passband: DC to ~0.6 × Nyquist (flat)
• Stopband: Moderate aliasing
• Transition: Smoother than linear

Hermite Interpolation

• Passband: DC to ~0.7 × Nyquist (flat)
• Stopband: Low aliasing
• Transition: Very smooth

Sinc Interpolation (16-tap)

• Passband: DC to ~0.95 × Nyquist (flat)
• Stopband: Minimal aliasing (<-80dB)
• Transition: Near-ideal brick-wall

---

Implementation Notes

Boundary Handling

All buffer-based interpolation functions (*_buffer_kernel) handle boundaries by: 1. Clamping position to valid range 2. Fallback to simpler method for short buffers (< required samples)

Example:

// Hermite requires 4 samples
float buffer[2] = {1.0f, 2.0f};

// Falls back to linear automatically
float val = hermite_buffer_kernel(buffer, 0.5f, 2);  // Uses lerp

SIMD Optimization

Batch functions are designed for auto-vectorization:

// Compiler can vectorize this loop (4-8 interpolations per iteration)
lerp_batch_kernel(buffer, positions, output, 1024, buffer_size);

Enable with: - GCC/Clang: -O3 -march=native - MSVC: /O2 /arch:AVX2

---

Testing

Run Tests

cd build
cmake --build .
ctest -V -R interpolation

Test Coverage

• ✅ Endpoint correctness (frac=0, frac=1)
• ✅ Fractional positions (0.25, 0.5, 0.75)
• ✅ Smoothness validation (derivative continuity)
• ✅ Frequency response (sine wave reconstruction)
• ✅ Edge cases (short buffers, boundary clamping)
• ✅ Quality comparison (all methods on same data)
• ✅ Batch processing
• ✅ Performance (large buffer stress test)

---

References

• Linear: Basic interpolation theory
• Cubic: Lagrange polynomial interpolation
• Hermite: Hermite spline interpolation (graphics/animation)
• Sinc: Shannon sampling theorem, Whittaker–Shannon interpolation formula
• Blackman window: F.J. Harris, "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform"

---

Future Enhancements

• Optimal fractional delay coefficients (Thiran allpass)
• Polynomial approximations for sinc (faster than actual sinc)
• SIMD intrinsics for batch interpolation
• Lookup tables for sinc window coefficients
• Additional window types (Kaiser, Hamming, Hann)

---

Recommendation: Use Hermite (InterpQuality::MEDIUM) as default. It offers the best quality/performance balance for most audio applications.