Multirate DSP
Est. read time: 3 minutes | Last updated: September 23, 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, 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 = 20 * M # Number of filter taps
fs = 48e3 # sampling frequency (Hz)
w_delta = np.pi/100 # 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} dB | Stopband attenuation: {sb_atten} dB")
Passband ripple: 1.0545818295681215 dB | Stopband attenuation: -24.335793669259637 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()
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
- A Review of Polyphase Filter Banks and Their Application - DTIC AFRL
- Polyphase Filters and Filterbanks - Kyle
YouTubeVideo('Zmjk9NE-3k0')
YouTubeVideo('afU9f5MuXr8')