Multirate DSP
Est. read time: 5 minutes | Last updated: December 17, 2024 by John Gentile
Contents
Multirate signal processing is the study and implementation of concepts, algorithms, and architectures that embed sample rate changes at one or more places in the signal data flow. Multirate signal processing is a subset of Digital Signal Processing (DSP) as we primarily deal with sampled, discrete systems (rather than analog, continuous-time systems). The two main reasons to implement Multirate signal processing in a design are:
- Reducing the cost of design implementation
- Enhancing the performance of the implementation
Some of the common tenets of Multirate DSP are:
- A processing task should always be performed at the lowest rate commensurate with the signal bandwidth (the Nyquist rate of the signal of interest (SOI))
- To achieve the above, Multirate DSP often deals with interchanging traditional DSP processing stages, and accepting intentional aliasing, by use of the Noble Identity
For example, a common DSP task is to reduce the bandwidth of a sampled signal to isolate a SOI (e.x. using a low-pass filter (LPF)), and then reduce the sample rate to match the SOI’s bandwidth. However with the Noble Identity, we can reduce the sample rate first, and then perform the filtering at the lower sample rate (and often with a lower order filter compared to it processing at the higher sample rate). This intentional aliasing can be unwrapped by the subsequent Multirate processing. This powerful tool can even be used to intentionally alias a signal from one center frequency to another (e.g. from an Intermediate Frequency (IF) to baseband) while reducing sample rate.
from IPython.display import YouTubeVideo
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from scipy import signal
from rfproto import filter, measurements, multirate, plot, sig_gen, utils
Sample Rate Conversion
Downsampling
Also known as decimation (though technically the “deci” was a historical term to “remove every tenth one”), downsampling is the process of reducing the sample rate of a system. To downsample by an integer factor , we can simply keep every th sample from an input stream, discarding the rest.
f_start = 20e3
f_end = 25e3
fs = 100e3
num_samples = 16384
lfm_chirp_sig = sig_gen.cmplx_dt_lfm_chirp(1, f_start, f_end, fs, num_samples)
lfm_chirp_sig = lfm_chirp_sig.real + 1j*np.zeros(num_samples)
plot.spec_an(lfm_chirp_sig, fs, fft_shift=True, show_SFDR=False)
plt.show()
M = 4
y_ds = multirate.decimate(lfm_chirp_sig, M)
plot.spec_an(y_ds, fs/M, fft_shift=True, show_SFDR=False)
plt.show()
Upsampling
Upsampling is the opposite of downsampling and is the process of increasing the sample rate of a system. To upsample by an integer factor , we can insert zeros between every sample from an input stream.
f_start = -5e3
f_end = 5e3
fs = 100e3
num_samples = 16384
lfm_chirp_sig = sig_gen.cmplx_dt_lfm_chirp(1, f_start, f_end, fs, num_samples)
plot.spec_an(lfm_chirp_sig, fs, fft_shift=True, show_SFDR=False)
plt.show()
L = 4
y_us = multirate.interpolate(lfm_chirp_sig, L)
plot.spec_an(y_us, fs*L, fft_shift=True, show_SFDR=False)
plt.show()
Noble Identities
The noble identity in multirate DSP allows us to commute resampling operations. Since down and upsamplers are Linear Time-Varying (LTV) processes, the order of operations matter- moving a downsampler through a filter is equivalent to a filter followed by the downsampling operation. The primary use of this identity is to allow processing of multirate functions (e.x. filtering) at the ideal/efficient sample rate:
Memoryless operations (e.x. adders, multipliers) may be commuted across down and upsamplers as well:
Polyphase Filtering
Polyphase Filters
M = 5
N = 10 * M # Number of filter taps
fs = 48e3 # sampling frequency (Hz)
w_delta = np.pi/20 # transition bandwidth
wp = ((np.pi/M) - w_delta)/np.pi # passband
ws = ((np.pi/M) + w_delta)/np.pi # stopband
# NOTE: vs MATLAB `firpm` which specifies band egdes from 0->fs,
# scipy.signal.remez() uses bands 0->fs/2
taps = signal.remez(N, [0, wp*fs/2, ws*fs/2, fs/2], [1, 0], fs=fs)
plot.filter_response(taps)
plot.plt.show()
pb_ripple, sb_atten = filter.measure_filter_response(taps, 0.0, wp, ws, 1.0)
print(f"Passband ripple: {pb_ripple:0.2f} dB | Stopband attenuation: {sb_atten:0.2f} dB")
Passband ripple: 0.08 dB | Stopband attenuation: -46.45 dB
h_poly = taps.reshape(len(taps)//M, M).T
# generate x indices
xi = np.zeros_like(h_poly)
color = iter(cm.rainbow(np.linspace(0, 1, M)))
for i in range(M):
xi[i] = list(range(i, N, M))
plt.stem(xi[i], h_poly[i], cm.viridis(i/M), label=f"h[{i}]")
plt.legend()
plt.show()
It can be seen that the response of each polyphase leg is allpass with linear phase from .
fig, axs = plt.subplots(2, 1, layout='constrained')
for i in range(M):
wi, hi = signal.freqz(h_poly[i])
axs[0].plot(wi / np.pi, utils.mag_to_dB(M * hi), linewidth=0.5, label=f"h[{i}]")
axs[1].plot(wi / np.pi, np.unwrap(np.angle(hi)), linewidth=0.5, label=f"h[{i}]")
axs[0].set_ylabel("Amplitude (dB)")
axs[0].set_xlabel(r"Normalized Frequency ($\times \pi$ rad/sample)", fontsize=12)
axs[0].margins(x=0)
axs[0].legend()
axs[1].set_ylabel("Angle (radians)")
axs[1].set_xlabel(r"Normalized Frequency ($\times \pi$ rad/sample)", fontsize=12)
axs[1].margins(x=0)
axs[1].legend()
plt.show()
f_start = 0
f_end = 0.1 * fs / 2.0
print(f"Passband from {f_start} -> {f_end} Hz")
num_samples = 100000
lfm_chirp_sig = np.real(sig_gen.cmplx_dt_lfm_chirp(1.0, f_start, f_end, fs, num_samples))
plot.spec_an(lfm_chirp_sig, fs, show_SFDR=False)
plot.plt.show()
Passband from 0 -> 2400.0 Hz
poly_filt_decim = multirate.polyphase_filter(M, taps * M, decimating=True)
poly_out = []
for x in lfm_chirp_sig:
tmp_out, input_consumed = poly_filt_decim.step(x)
if tmp_out is not None:
poly_out.append(tmp_out)
plot.spec_an(poly_out, fs/M, show_SFDR=False)
plot.plt.show()
plot.spec_an(multirate.interpolate(lfm_chirp_sig, M), fs*M, show_SFDR=False)
plot.plt.show()
poly_filt_interp = multirate.polyphase_filter(M, taps, decimating=False)
poly_out = []
input_idx = 0
while input_idx < len(lfm_chirp_sig):
#for _ in range(24):
tmp_out, input_consumed = poly_filt_interp.step(lfm_chirp_sig[input_idx])
if tmp_out is not None:
poly_out.append(tmp_out)
if input_consumed:
input_idx += 1
plot.spec_an(poly_out, fs*M, show_SFDR=False)
plot.plt.show()
Polyphase Filter Banks (PFBs) and Channelizers
PFB References
- Xilinx XAPP1161 Polyphase Filter Bank Channelizer, v1.0, Application Note
- Spectrometers and Polyphase Filterbanks in Radio Astronomy - Danny C. Price - arXiv
YouTubeVideo('WRO9VD6MYy4')
Polyphase Filter References
- Polyphase Filters - Section 12.4 Porat Notes
- A Review of Polyphase Filter Banks and Their Application - DTIC AFRL
- VHDL whiz: Part 5 - Polyphase FIR Filters
- Polyphase Filters and Filterbanks - Kyle
YouTubeVideo('Zmjk9NE-3k0')
YouTubeVideo('afU9f5MuXr8')
Cascaded Integrator–Comb (CIC) Filter
Benefits:
- Multiplier-free implementation for economical (light resource utilization) design in digital HW systems which can handle arbitrary, and large, rate changes.
- High interpolation CIC filters can push very narrowband signals (e.x. TT&C ~2MS/s max) go to sufficient ADC/DAC sample rates (e.x. 125MS/s)
- Data reduction to not have to pass as much data to/from the front end (e.g. having to push full sample bandwidth data to/from a remote radio head)
CIC filters can be used as decimation (decrease sample rate) and interpolation (increase sample rate) multirate filters:
The basic building blocks are integrator and comb sections (hence the name), along with an interpolator (e.g. zero-stuffing expander) or decimator, increasing or decreasing the output sample rate by times (respectively).
Integrator Section
An integrator is simply a single-pole IIR filter with unity feedback coefficient:
This is also commonly known as an accumulator, which has the -transform transfer function of:
The power response is basically a low-pass filter with a −20 dB per decade (−6 dB per octave) rolloff, but with infinite gain at DC. This is due to the single pole at z = 1; the output can grow without bound for a bounded input. In other words, a single integrator by itself is unstable. - CIC Filter Introduction - Matthew P. Donadio (DSP Guru)
Note the pipelined version of the integrator which allows for more efficient digital hardware implementation (only one added between register stages):
Comb Section
The comb filter runs at the highest sample rate, , in either decimation or interpolation filter form. It looks opposite of an integrator section as it subtracts the current sample value from a value sample periods prior; is the differential delay design parameter, and is often limited to or :
The corresponding transfer function at is:
When and , the power response is a high-pass function with 20dB per decade gain (inverse of the integrator response). When , the power response looks like a familiar raised cosince form with cycles from .
Bit Growth
Due to the cascaded adders/subtracters in the CIC filter, each fixed-point, two’s complement arithmetic operation requires an additional bit of output than input, to prevent loss of precision. Given an input sample bitwidth of , the output bitwidth required can be found to be:
At the expense of added quantization noise, bit growth can be controlled by rounding/scaling at some points within the CIC stages.
#TODO:
- Look at the fred harris paper on multiplier-less CIC with sharpening that doesn’t need compensation FIR
- Also mentioned in https://www.dsprelated.com/showarticle/1337.php
- For bit growth, look at CIC filter register pruning
- Also shown in http://www.jks.com/cic/cic.html and https://github.com/jks-prv/cic_prune
Discrete-Time Test of CIC
f_start = 0
f_end = 5e3
fs = 100e3
num_samples = 100000
bit_width = 16
lfm_chirp_sig = np.real(sig_gen.cmplx_dt_lfm_chirp(2**(bit_width-5), f_start, f_end, fs, num_samples))
freq, y_PSD = measurements.PSD(lfm_chirp_sig, fs, norm=True)
plot.freq_sig(freq, y_PSD, "LFM Chirp Input Spectrum", scale_noise=False)
plt.show()
N = 3 # number of stages
R = 8 # interp/decim factor
M = 1 # differential delay in comb stages
cic_bit_width = 21 # np.ceil(N * np.log2(R*M) + bit_width)
integ_stages = [multirate.integrator(cic_bit_width) for i in range(N)]
comb_stages = [multirate.comb(M) for i in range(N)]
integ_out = np.zeros(len(lfm_chirp_sig))
for i in range(len(integ_out)):
temp_val = int(lfm_chirp_sig[i])
for j in range(N):
temp_val = integ_stages[j].step(temp_val)
integ_out[i] = temp_val
decimate_out = multirate.decimate(integ_out, R)
cic_out = np.zeros(len(decimate_out))
for i in range(len(cic_out)):
temp_val = decimate_out[i]
for j in range(N):
temp_val = comb_stages[j].step(temp_val)
cic_out[i] = temp_val
freq, y_PSD = measurements.PSD(cic_out, fs / R, norm=True)
plot.freq_sig(freq, y_PSD, "CIC Output Spectrum", scale_noise=False)
plt.show()
CIC References
- CIC Filter Introduction - Matthew P. Donadio (DSP Guru)
- A Beginner’s Guide to Cascaded Integrator-Comb (CIC) Filters
- Small Tutorial on CIC Filters
- CIC Filter - Wikipedia
- CIC Compiler v4.0 - Xilinx
- CIC Filter - OpenCores