Home > On-Demand Archives > Workshops >
The Space Between Samples: Variable Delay with Subsample Precision
Dan Boschen - Watch Now - DSP Online Conference 2025 - Duration: 02:02:45
This workshop offers a practical, application-driven introduction to implementing programmable fractional delay lines with exceptional precision. We'll begin by reviewing fundamental interpolation approaches, then move to a structured design methodology for developing efficient, high-performance fractional delays, with an emphasis on **polyphase** and **Farrow filter architectures**.
Real-world use cases, including **beamforming**, **timing recovery**, and **high-precision EVM measurement**, will illustrate how these techniques can be applied in practice. The session will also feature a live Python coding demonstration to show you how to implement practical fractional delay filters.
Participants with their own laptops will have the option to follow along, while others can simply observe the walkthrough. All attendees will receive installation instructions and example scripts to support live experimentation and future use.
This guide was created with the help of AI, based on the presentation's transcript. Its goal is to give you useful context and background so you can get the most out of the session.
What this presentation is about and why it matters
The Space Between Samples: Variable Delay with Subsample Precision is a practical workshop on creating programmable fractional delays — the ability to delay a digital waveform by an arbitrary fraction of a sample. This capability is fundamental in many engineering systems: array beamforming (steering by time delays across elements), receiver timing recovery, precise waveform alignment for EVM and distortion measurements, and digital resampling. Fractional delays let you shift the underlying analog-equivalent waveform at sub-sample resolution without changing the system sample rate, which is essential when you need high-precision alignment or low-distortion interpolation.
Who will benefit the most from this presentation
- DSP engineers working on wireless communications, radar, or audio spatial processing who need precise timing/phase alignment.
- Test and measurement engineers implementing high-dynamic-range waveform comparison and EVM analysis.
- Students and practitioners learning multirate signal processing, interpolation, and filter design who want efficient implementation patterns (polyphase and Farrow).
- Engineers choosing trade-offs between memory and computation in real-time systems (FPGA/ASIC/embedded DSP implementations).
What you need to know
This talk assumes basic DSP familiarity but focuses on intuition and practical design steps. To get the most out of the presentation, be comfortable with the following concepts:
- Sampling and the Nyquist concept: why signals above half the sample rate alias and why anti-alias filtering matters for resampling.
- FIR filters and impulse response: an FIR filter is a weighted sum of delayed samples; the coefficients are the impulse response.
- Sinc interpolation: the ideal low-pass (reconstruction) filter has a sinc impulse response; truncated+windowed sinc FIRs are the practical interpolators used for fractional delay.
- Windowing: truncating a sinc introduces ripples; a taper (e.g., Kaiser window) reduces sidelobes and passband ripple.
- Polyphase concept: upsample-by-N then filter is equivalent to computing N sub-filters (polyphase components) running at the input rate — a very efficient way to implement interpolation and fractional delay without processing at the full upsample rate.
- Farrow structure: a method to compress a large polyphase coefficient table by modeling each coefficient column with a polynomial (e.g., Chebyshev, Lagrange). This trades extra arithmetic for much less memory while providing continuous fractional delay.
- Group delay: the effective time delay of a linear-phase filter is the derivative of phase versus frequency. In MathJax form: $\tau_g = -\dfrac{d\phi}{d\omega}$. For fractional-delay filters, you want flat group delay across the passband where your signal energy lies.
- Design trade-offs: attenuation (stopband rejection), transition width, number of taps, and number of polyphase banks (precision) determine memory, compute, and distortion performance.
Glossary
- Fractional delay — delaying a signal by a non-integer number of samples (subsample shift of the underlying waveform).
- Sinc interpolation — using samples of the sinc function (ideal low-pass kernel) to reconstruct intermediate sample values.
- Polyphase — decomposition of an upsample+filter into N parallel filters (phases) operating at the original sample rate.
- Farrow filter — polynomial-based structure that generates polyphase coefficients on-the-fly, reducing memory at the cost of computation.
- Windowed sinc — a sinc function multiplied by a window (e.g., Kaiser) to reduce truncation effects and control sidelobes.
- Group delay — the frequency-dependent delay of a filter; approximate linear phase produces constant group delay across the passband.
- Taps — the number of coefficients in an FIR filter; more taps yield steeper transitions and higher stopband rejection (at cost of delay and compute).
- Upsampling / zero-insertion — inserting zeros between samples to increase the digital sampling rate before filtering to interpolate.
- Anti-aliasing / pre-filter — filter applied before downsampling to prevent high-frequency content from folding into the passband.
- Chebyshev polynomial — a polynomial basis useful for minimizing peak error when approximating coefficient curves in Farrow-like compression.
Final notes — why watch this talk
Dan Boschen delivers a very practical, implementation-focused tour of fractional-delay techniques. The talk is strong on intuition (visualizing shifted sinc kernels and polyphase mapping), provides clear trade-offs (memory versus computation, odd vs even tap choices), and demonstrates how to go from prototype filter to either a lookup-table polyphase implementation or a compressed Farrow form. If you like hands-on learning, the included Jupyter notebook and step-by-step recipe (Kaiser design parameters, taps-per-bank, Chebyshev/Farrow fitting) are particularly valuable: you can reproduce and experiment with the exact designs shown. Expect a concise mix of theory, engineering heuristics, and actionable Python code — a good blend for both learners and practitioners.
Enjoy the workshop — it will sharpen your intuition about "the space between samples" and give you practical tools to implement sub-sample delays reliably in real systems.
Oh very nice! I am familiar with MATLAB's firpm function but have never seen this pm_remez.remez function you show. That's a nice feature.
Yes I love Python since I started using it 8 years ago. Consider jumping in on my Python for DSP course when it starts again next August. It's very easy now to have ChatGPT show you all the code you need, but even with that it's very good to understand what it is showing you as well as the best practice approaches and common gotchas that frustrate most newcomers when they don't understand. That part I can get you through for DSP very efficiently. Further, you can sign up now and get early access to all the videos (and have up to a year of access for that) and still get in on the live workshop sessions next year). See details here: https://dsprelated.com/courses
Hi Dan. Thank you for an excellent talk. You talked about minimizing the aliasing of the stopbands into the passband when decimating by not using a flat stopband. Prof. Harris does this in the Remez design by inserting his own custom cost function, myfrf.m. This gives a stopband rolloff of 1/f. Is there a way to do the same thing as myfrf.m does but in Python Remez design? If so, can you give an example syntax of the normal Remez function call vs calling with the custom 1/f stopband rolloff function?
I asked chatGPT to port the Moerder and harris myfrf.m available from here: https://github.com/OpenResearchInstitute/polyphase-filter-bank/blob/master/polyphase_filter_bank_resources_from_fred_harris/MATLAB_scripts/myfrf.m
or here:
https://www.dsprelated.com/thread/9593/fft-channelizer-and-pfb
I haven't tested this. If you do, please let us know how well it works!
from typing import Tuple, Union
def myfrf(N: Union[int, str],
F: np.ndarray,
GF: np.ndarray,
W: np.ndarray,
A: np.ndarray,
diff_flag: int = 0) -> Tuple[np.ndarray, np.ndarray]:
"""
Python port of fred harris / Moerder 'myfrf.m' for REMEZ.
Parameters
----------
N : int or str
Filter order (ignored here except for the 'defaults' query).
If N == 'defaults', returns ('even', None) to match MATLAB query.
F : (2*nbands,) array_like
Band edges (normalized 0..1). Pairs (F[0],F[1]), (F[2],F[3]), ...
GF : (K,) array_like
Interpolated grid frequencies (normalized 0..1).
W : (nbands,) array_like
Weights, one per band.
A : (2*nbands,) array_like
Desired amplitudes at band edges F.
diff_flag : int
Must be 0. Differentiator option is explicitly disallowed.
Returns
-------
DH : (K,) ndarray
Desired frequency response on GF (magnitude; phase is linear-phase external to this routine).
DW : (K,) ndarray
Weights on GF (positive).
Notes
-----
- Matches the MATLAB logic:
* Linear interpolation for desired magnitude within each band.
* If A(l+1) > 1e-4 (treated as passband): DW = W(band)
else (stopband): DW = W(band) * (GF / GF_band_start)
* Disallows adjacent bands (F[k] == F[k+1]).
* Requires len(F) == len(A).
- Symmetry query: if called as myfrf('defaults', F), returns 'even'.
"""
# Support the symmetry query used by REMEZ in MATLAB
if isinstance(N, str) and N == 'defaults':
# Return symmetry default (even or odd). MATLAB returns just 'even'.
return 'even', None # type: ignore[return-value]
# Basic checks
F = np.asarray(F, dtype=float).ravel()
GF = np.asarray(GF, dtype=float).ravel()
W = np.asarray(W, dtype=float).ravel()
A = np.asarray(A, dtype=float).ravel()
if diff_flag != 0:
raise ValueError("Differentiator option is not allowed with myfrf.")
if F.size != A.size:
raise ValueError("Frequency (F) and amplitude (A) vectors must be the same length.")
if F.size % 2 != 0:
raise ValueError("F and A must contain pairs of band edges; length must be even.")
# Prevent discontinuities in desired function (adjacent bands with zero width)
for k in range(1, F.size - 2, 2):
if F[k] == F[k + 1]:
raise ValueError("Adjacent bands not allowed (F[k] equals F[k+1]).")
nbands = F.size // 2
if W.size != nbands:
raise ValueError("W must have one weight per band (len(W) == len(F)/2).")
DH = np.zeros_like(GF, dtype=float)
DW = np.zeros_like(GF, dtype=float)
l = 0
while (l + 1) // 2 < nbands:
f0, f1 = F[l], F[l + 1]
a0, a1 = A[l], A[l + 1]
band_idx = (l + 1) // 2 # 0-based band index
sel = np.where((GF >= f0) & (GF <= f1))[0]
if sel.size > 0:
if f1 != f0:
slope = (a1 - a0) / (f1 - f0)
# DH(sel) = slope*GF(sel) + (a0 - slope*f0)
DH[sel] = slope * GF[sel] + (a0 - slope * f0)
else:
# Zero-bandwidth band: average the amplitudes
DH[sel] = 0.5 * (a0 + a1)
# Weights:
# Passband if a1 > 1e-4, else stopband with frequency-scaled weight
if a1 > 1e-4:
DW[sel] = W[band_idx]
else:
# MATLAB: DW(sel) = W(band) * (GF(sel)/GF(sel(1)))
# Guard against division by zero if GF[sel][0] == 0
gf0 = GF[sel][0]
if gf0 == 0.0:
# If the first grid point in this band is 0 Hz,
# approximate scaling by setting the first point's scale to 1
# and others relative to max(very small, gf0) to avoid /0.
# This mirrors the original intent without throwing.
scale = np.ones_like(GF[sel])
nz = GF[sel] > 0
if np.any(nz):
scale[nz] = GF[sel][nz] / GF[sel][nz][0]
else:
scale = GF[sel] / gf0
DW[sel] = W[band_idx] * scale
l += 2
return DH, DWThanks Michael. I don't believe the Python remez function supports that, mirroring what we experience in MATLAB's firpm and Octaves remez functions. See this post where I provide a link to myfrf.m if you didn't already have it: https://dsp.stackexchange.com/q/95103/21048 as well as my own observations of using Parks McClellan for these applications. I did ask ChatGPT to port it to Python and posted in a second comment (but haven't tested it).
I mention in the StackExchange post that myfrf.m was inferior in my experience to using least squares. Also noted, I added an updated answer recently showing one example that Professor harris recently gave me where Parks McClellan (remez) would actually be a superior design approach for the design of differentiators specifically (I have been searching for years and in all other cases, the least squares and kaiser windowed design approaches were superior). That said, I really like MattL's answer showing his compromise between least squares and Parks McClellan as a constrained least squares.
However, bottom line for me if you need a linear-phase lowpass filter with stopband roll-off and superior noise performance, use firls (least squares), and if that doesn't converge for some corner cases, then fall-back to a Kaiser windowed Sinc function. If you (or anyone else) have an actual practical example where Parks McClellan would actually be superior, kindly add that to the StackExchange post I linked.
Thanks for the great talk! Could you maybe elaborate a bit more on the advantage of using Chebyshev polynomials rather than Farrow's original formulation?
Great question. The consideration would be for applications such as what I demonstrated in the notebook for very high precision where higher order polynomials are needed, and comes down to approximation techniques for any function: we can do a power series as a Taylor series expansion where our basis “functions” used for interpolation are 1, x, x^2, x^3 etc up to the order of the polynomial used. This is what we typically associate with polynomial interpolation and is what you find in the classic Farrow structure: we estimate a function as a sum of weighted powers. But we can also use other basis functions to approximate functions (the Fourier Series is a great example where we use the sum of weighted sinusoids instead of the sum of weighted powers). So with the Chebyshev approach we use Chebyshev polynomials limited to the range of +/-1 to approximate a function over that range. Why? The error is uniform over that entire range (while with a Taylor series the error is small local to the center of the range), Runge’s phenomenon (ringing) can be severe with Taylor and much supressed with Chebyshev, and Chebyshev polynomials are orthogonal while a power series is not - which leads to high numerical stability even with large orders. So Chebyshev is a great choice for any applications for approximating waveforms over finite intervals with high accuracy. Further as I show in the appendix, any of the values can be determined easily through iteration so it is computationally efficient. The other polynomial choices such as Lagrange, Hermite, Legendre are also good choices as summarized in the table I gave.
Thank you Dan. I do have the myfrf.m function, generously provided to me by Prof harris, years ago. I need to try Python DSP as you show in your lectures. I haven't up to now because I have access to Matlab. However, I didn't realize how powerful Python was for DSP and it is something I need to get running. After I do, I will try that ChatGPT code and see if i can get the myfrf cost function to work in the Python remez as well.
In Googling this, it seems that the pm-remez function will accept a callable function like myfrf as the weights. Here is the response from Google AI "To use the pm-remez function with a custom 1/f stopband cost function, you need to pass a Python callable (function) as the weight argument for the stopband region when calling the pm_remez function. The pm-remez library is designed to accept arbitrary functions for both the desired response and the weights."
The code it suggested:
import numpy as np
import pm_remez
def stopband_weight_1f(f):
"""
Calculates 1/f weight for the stopband.
Avoids division by zero or near-zero frequencies if they are in the stopband start.
"""
Add a small offset (epsilon) to avoid issues at f=0 if necessary
Filter parameters
num_taps = 51 # Example filter order
bands = [0.0, 0.1, 0.2, 0.5] # Normalized frequency bands (0 to 0.5)
Desired amplitude in each band: [Passband_amplitude, Stopband_amplitude]
desired = [1.0, 0.0]
Weights: Passband weight (e.g., 1.0) and Stopband weight function
The function is passed as a callable object
weights = [1.0, stopband_weight_1f]
Design the filter
h = pm_remez.remez(num_taps, bands, desired, weights)
'h' now contains the filter coefficients
print(f"Filter coefficients: {h}")