Home > On-Demand Archives > Workshops >
The Space Between Samples: Variable Delay with Subsample Precision
Dan Boschen - Watch Now - EOC 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.
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}")