Skip to content
Snippets Groups Projects
Commit eb19d4cc authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Add up(), down(), downsample_bpf().

parent 8f6e81aa
No related branches found
No related tags found
1 merge request!419Resolve RTSD-265
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/rtdsp/multirate.py
+ 228
− 93
View file @ eb19d4cc
......@@ -23,14 +23,42 @@
# Purpose: Multirate functions for DSP
# Description:
#
# Usage:
# . Use [2] to verify this code.
#
# References:
# [1] dsp_study_erko.txt
# [2] http://localhost:8888/notebooks/pfb_os/up_downsampling.ipynb
#
# Books:
# . HARRIS 6 title Figure
# . LYONS
# . CROCHIERE
import numpy as np
from scipy import signal
from .utilities import c_rtol, c_atol, ceil_div
def down(x, D, phase=0):
"""Downsample x[n] by factor D, xD[m] = x[m D], m = n // D
Keep every D-th sample in x and skip every D-1 samples after that sample.
"""
xD = x[phase::D]
return xD
def up(x, U):
"""Upsample x by factor U, xU[m] = x[n] when m = n U, else 0.
After every sample in x insert U-1 zeros.
"""
xU = np.zeros(len(x) * U)
xU[0::U] = x
return xU
###############################################################################
# Polyphase filter
###############################################################################
......@@ -55,25 +83,19 @@ class PolyPhaseFirFilterStructure:
of data block with newest sample at last index. Flipping polyCoefs once
here at init is probably more efficient, then flipping inData and outData
for every block.
. No need to flip the order of the Ntaps columns (axis=1), because the
filter sums the Ntaps.
Filter delay memory polyDelays:
. Shift in data blocks per column at column 0, last sample in data block
and in column is newest.
"""
def __init__(self, Nphases, coefs, commutator):
self.Nphases = Nphases
self.coefs = coefs
self.Ncoefs = len(coefs)
self.Ntaps = ceil_div(self.Ncoefs, Nphases)
self.Nphases = Nphases # branches, rows
self.Ntaps = ceil_div(self.Ncoefs, Nphases) # taps, columns
# coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
# polyCoefs = [[ 3, 7, 11],
# [ 2, 6, 10],
# [ 1, 5, 9],
# [ 0, 4, 8]]), Nphases = 4 rows (axis=0)
# Ntaps = 3 columns (axis=1)
# Apply polyCoefs branches in range(Nphases) order, because coefs
# mapping in poly_coeffs(commutator) already accounts for down or up.
self.polyCoefs = self.poly_coeffs(commutator)
# oldest sample[0]
......@@ -92,18 +114,36 @@ class PolyPhaseFirFilterStructure:
"""Polyphase structure for FIR filter coefficients coefs in Nphases
Input:
. coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
. coefs = prototype FIR filter coefficients (impulse response)
. commutator: 'up' or 'down'
Return:
. commutator: 'up', see title Figure in [HARRIS 6], Figure 10.22 and
10.23 in [LYONS] seem to have commutator arrows in wrong direction.
. polyCoefs:
. E.g. coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
. If Ncoefs < Ntaps * Nphases, then pad polyCoefs with zeros.
. Branch order of polyCoefs accounts for commutator, therefore the
polyCoefs branches can be applied in range(Nphases) order
independent of commutator direction.
. commutator: 'up', first output appears when the switch is at the
branch p = 0 and then yields a new output for every branch,
therefore p = 0 first in polyCoefs.
polyCoefs = [[ 0, 4, 8], # p = 0, H0(z)
[ 1, 5, 9], # p = 1, H1(z)
[ 2, 6, 10], # p = 2, H2(z)
[ 3, 7, 11]]) # p = 3, H3(z),
Nphases = 4 rows (axis=0)
Ntaps = 3 columns (axis=1)
If Ncoefs < Ntaps * Nphases, then pad polyCoefs with zeros.
. commutator: 'down', new FIR sum output starts being aggregated when
switch is at branch p = Nphases-1, therefore p = Nphases - 1 first
in polyCoefs
polyCoefs = [[ 3, 7, 11], # p = 3, H3(z)
[ 2, 6, 10], # p = 2, H2(z)
[ 1, 5, 9], # p = 1, H1(z)
[ 0, 4, 8]]) # p = 0, H0(z)
Nphases = 4 rows (axis=0)
Ntaps = 3 columns (axis=1)
"""
polyCoefs = np.zeros(self.Nphases * self.Ntaps)
polyCoefs[0 : self.Ncoefs] = self.coefs
......@@ -116,12 +156,12 @@ class PolyPhaseFirFilterStructure:
"""Filter block of inData per poly phase
Input:
. inData: block of Nphases input data, time n increments, so inData[0] is
oldest data and inData[Nphases-1] is newest data
. inData: block of Nphases input data, time index n increments, so
inData[0] is oldest data and inData[Nphases-1] is newest data.
Return:
. outData: block of Nphases filtered output data with outData[0] is
oldest data and outData[Nphases-1] is newest data, so time n
increments, like with inData
oldest data and outData[Nphases-1] is newest data, so time index n
increments, like with inData.
"""
# Shift delay memory by one block
self.polyDelays = np.roll(self.polyDelays, 1, axis=1)
......@@ -129,18 +169,61 @@ class PolyPhaseFirFilterStructure:
self.polyDelays[:, 0] = inData
# Apply FIR coefs per element
zData = self.polyDelays * self.polyCoefs
# Sum per poly phase
# Sum FIR taps per poly phase
outData = np.sum(zData, axis=1)
return outData
def poly_structure_size_for_downsampling_whole_x(Lx, Ndown):
"""Polyphase structure size for downsampling a signal x with Lx samples
The polyphase structure has Ndown branches and branch size suitable to use lfilter() once per polyphase branch
for whole signal x. Instead of using pfs.polyDelays per pfs.Ntaps.
. Prepend x with Ndown - 1 zeros to have first down sampled sample at m = 0 start at n = 0.
. Skip any remaining last samples from x, that are not enough yield a new output FIR sum.
Input:
. Lx: len(x)
. Ndown: downsample rate
Return:
. Nx: Total number of samples from x, including prepended Ndown - 1 zeros.
. Nxp: Total number of samples used from x per polyphase branch, is Ny the number of samples that will be in
downsampled output y[m] = x[m D] for m = 0, 1, 2, ..., Nxp - 1
"""
# Numer of samples per polyphase
Nxp = (Ndown - 1 + Lx) // Ndown
# Used number of samples from x
Nx = 1 + Ndown * (Nxp - 1)
return Nx, Nxp
def poly_structure_data_for_downsampling_whole_x(x, Ndown):
"""Polyphase data structure for whole signal x
. see poly_structure_size_for_downsampling_whole_x()
Input:
. x: input samples
Return:
. polyX: polyphase data structure with size (Ndown, Nxp) for Ndown branches and Nxp samples from x per branch.
"""
Lx = len(x)
Nx, Nxp = poly_structure_size_for_downsampling_whole_x(Lx, Ndown)
polyX = np.zeros(Ndown * Nxp)
polyX[Ndown - 1] = x[0] # prepend x with Ndown - 1 zeros
polyX[Ndown:] = x[1 : Nx]
polyX = polyX.reshape(Nxp, Ndown).T
return polyX, Nx, Nxp
def upsample(x, Nup, coefs, verify=False, verbosity=1): # interpolate
"""Upsample x by factor I = Nup
"""Upsample x by factor U = Nup and LPF
Input:
. x: Input signal x[n]
. Nup: upsample (interpolation) factor
. coefs: FIR filter coefficients for antialiasing LPF
. coefs: prototype FIR filter coefficients for anti aliasing LPF
. verify: when True then verify that output y is the same when calculated directly or when calculated using the
polyphase implementation.
- verbosity: when > 0 print() status, else no print()
......@@ -149,34 +232,34 @@ def upsample(x, Nup, coefs, verify=False, verbosity=1): # interpolate
Assumptions:
. x[n] = 0 for n < 0
. m = n * I, so first upsampled output sample starts at m = 0 based on n = 0
. m = n * U, so first upsampled output sample starts at m = 0 based on n = 0
. no y[m] for m < 0
. insert I - 1 zeros after each x[n] to get v[m], so len(v) = len(y) = I * len(x)
. use coefs as anti aliasing filter, must be LPF with BW < fNyquist / I
. insert U - 1 zeros after each x[n] to get v[m], so len(v) = len(y) = U * len(x)
. use coefs as anti aliasing filter, must be LPF with BW < fNyquist / U
. len(coefs) typically is multiple of Nup. If shorter, then the coefs are extended with zeros.
Remarks:
. The input sample period is ts and the output sample period of the upsampled (= interpolated signal) is tsUp =
ts / I
. The group delay is (Ncoefs - 1) / 2 * tsUp = (Ncoefs - 1) / I / 2 * ts. With odd Ncoefs and symmetrical coefs
ts / U
. The group delay is (Ncoefs - 1) / 2 * tsUp = (Ncoefs - 1) / U / 2 * ts. With odd Ncoefs and symmetrical coefs
to have linear phase, the group delay is an integer number of tsUp periods, so:
- when (Ncoefs - 1) % 2 == 0, then the group delay is an integer number of tsUp periods
- when (Ncoefs - 1) % (2 * I) == 0, then the group delay is an integer number of ts periods
- when (Ncoefs - 1) % (2 * U) == 0, then the group delay is an integer number of ts periods
Procedure:
x[n]: 0 1 2 3 --> time n with unit Ts
v[m] = x[m / I], for m = 0, +-I, +-2I, ...
v[m] = x[m / U], for m = 0, +-U, +-2U, ...
= 0, else
v[m]: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 --> time m with unit Ts / I
v[m]: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 --> time m with unit Ts / U
| | | |
h[k]: 11 10 9 8 7 6 5 4 3 2 1 0 | | | --> coefs for m = 12
11 10 9 8 7 6 5 4 3 2 1 0 | | --> coefs for m = 13
11 10 9 8 7 6 5 4 3 2 1 0 | --> coefs for m = 14
11 10 9 8 7 6 5 4 3 2 1 0 --> coefs for m = 15
y[m]: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 --> time m with unit Ts / I
y[m]: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 --> time m with unit Ts / U
Calculate y[0, 1, 2, 3] at x[0]
y[0] = h[0] x[0]
......@@ -202,15 +285,15 @@ def upsample(x, Nup, coefs, verify=False, verbosity=1): # interpolate
y[14] = h[2] x[3] + h[6] x[2] + h[10] x[1]
y[15] = h[3] x[3] + h[7] x[2] + h[11] x[1]
==> Calculate y[n * I + [0, ..., I - 1]] at x[n]
==> Calculate y[n * U + [0, ..., U - 1]] at x[n]
y[n * I + 0] = lfilter(h[0, 4, 8], [1], x) # p = 0
y[n * I + 1] = lfilter(h[1, 5, 9], [1], x) # p = 1
y[n * I + 2] = lfilter(h[2, 6,10], [1], x) # p = 2
y[n * I + 3] = lfilter(h[3, 7,11], [1], x) # p = 3
y[n * U + 0] = lfilter(h[0, 4, 8], [1], x) # p = 0
y[n * U + 1] = lfilter(h[1, 5, 9], [1], x) # p = 1
y[n * U + 2] = lfilter(h[2, 6,10], [1], x) # p = 2
y[n * U + 3] = lfilter(h[3, 7,11], [1], x) # p = 3
"""
Nx = len(x)
a = [1.0]
a = [1.0] # FIR b = coefs, a = 1
# Polyphase implementation to avoid multiply by zero values
# coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
......@@ -222,14 +305,17 @@ def upsample(x, Nup, coefs, verify=False, verbosity=1): # interpolate
polyCoefs = pfs.polyCoefs
# Poly phases for data
# . Use polyY for whole data signal x with Nx samples, instead of pfs.polyDelays for pfs.Ntaps, to be able to use
# signal.lfilter()
# . Use polyX for whole data signal x with Nx samples, instead of pfs.polyDelays per pfs.Ntaps, to be able to use
# signal.lfilter() per branch for whole data signal x.
polyY = np.zeros((Nup, Nx))
# Filter x per polyphase, index p = 0, 1, .., Nup - 1 is the commutator in [Figure 10.13, 10.23 in LYONS, 3.25
# HARRIS]
# Filter x per polyphase, index p = 0, 1, .., Nup - 1, order in polyCoefs already accounts for the commutator
# direction [LYONS Fig 10.13, 10.23, 10.23 seems to have commutator arrows in wrong direction], [CROCHIERE Fig
# 3.25], [HARRIS 7.11, 7.12]
for p in range(Nup):
polyY[p] = signal.lfilter(polyCoefs[p], a, x)
# Output Nup samples for every input sample
y = polyY.T.reshape(1, Nup * Nx)[0]
if verify:
......@@ -245,20 +331,20 @@ def upsample(x, Nup, coefs, verify=False, verbosity=1): # interpolate
if verbosity:
print('> upsample():')
print('. Nup = ' + str(Nup))
print('. Nx = ' + str(Nx))
print('. len(y) = ' + str(len(y)))
print('. Nup =', str(Nup))
print('. Nx =', str(Nx))
print('. len(y) =', str(len(y)))
print('')
return y
def downsample(x, Ndown, coefs, verify=False, verbosity=1): # decimate
"""Downsample x by factor D = Ndown up
"""LPF and downsample x by factor D = Ndown
Input:
. x: Input signal x[n]
. Ndown: downsample (decimation) factor
. coefs: FIR filter coefficients for antialiasing LPF
. coefs: prototype FIR filter coefficients for anti aliasing LPF
. verify: when True then verify that output y is the same when calculated directly or when calculated using the
polyphase implementation.
- verbosity: when > 0 print() status, else no print()
......@@ -267,7 +353,7 @@ def downsample(x, Ndown, coefs, verify=False, verbosity=1): # decimate
Assumptions:
. x[n] = 0 for n < 0
. m = n // I, so first downsampled output sample starts at m = 0 based on n = 0
. m = n // D, so first downsampled output sample starts at m = 0 based on n = 0
. no y[m] for m < 0
. skip D - 1 values zeros after each x[n] to get y[m], so len(y) = len(x) // D
. use coefs as anti aliasing filter, must be LPF with BW < fNyquist / D
......@@ -329,12 +415,7 @@ def downsample(x, Ndown, coefs, verify=False, verbosity=1): # decimate
v[n - 1] = lfilter(h[1, 5, 9], [1], polyX[2])
v[n - 0] = lfilter(h[0, 4, 8], [1], polyX[3])
"""
# Number of samples from x per poly phase in polyX, prepended with Ndown - 1 zeros. Prepend Ndown - 1 zeros
# to have first down sampled sample at m = 0 start at n = 0.
Nxp = (Ndown - 1 + len(x)) // Ndown
# Used number of samples from x
Nx = 1 + Ndown * (Nxp - 1)
a = [1.0]
a = [1.0] # FIR b = coefs, a = 1
# Polyphase implementation to avoid calculating values that are removed
# coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
......@@ -345,19 +426,16 @@ def downsample(x, Ndown, coefs, verify=False, verbosity=1): # decimate
pfs = PolyPhaseFirFilterStructure(Ndown, coefs, commutator='down')
polyCoefs = pfs.polyCoefs
# Poly phases for data
# . Use polyX for whole data signal x with Nx samples, instead of pfs.polyDelays for pfs.Ntaps, to be able to use
# signal.lfilter() and to be able to prepend Ndown - 1 zeros
polyX = np.zeros(Ndown * Nxp)
polyX[Ndown - 1] = x[0] # prepend x with Ndown - 1 zeros
polyX[Ndown:] = x[1 : Nx]
polyX = polyX.reshape(Nxp, Ndown).T
polyY = np.zeros((Ndown, Nxp))
# Poly phases for whole data signal x, prepended with Ndown - 1 zeros
polyX, Nx, Nxp = poly_structure_data_for_downsampling_whole_x(x, Ndown) # size (Ndown, Nxp), Ny = Nxp
# Filter x per polyphase, index p = 0, 1, .., Ndown - 1 is the commutator in [Figure 10.14, 10.22 in LYONS, 3.29
# HARRIS]
# Filter x per polyphase, index p = 0, 1, .., Ndown - 1, the row order in polyCoefs already accounts for the
# commutator direction. [HARRIS Fig 6.12, 6.13], [LYONS Fig 10.14, 10.22], [CROCHIERE Fig 3.29]
polyY = np.zeros((Ndown, Nxp))
for p in range(Ndown):
polyY[p] = signal.lfilter(polyCoefs[p], a, polyX[p])
# Sum the branch outputs to get single downsampled output value
y = np.sum(polyY, axis=0)
if verify:
......@@ -373,17 +451,75 @@ def downsample(x, Ndown, coefs, verify=False, verbosity=1): # decimate
if verbosity:
print('> downsample():')
print('. len(x) = ' + str(len(x)))
print('. Ndown = ' + str(Ndown))
print('. Nx = ' + str(Nx))
print('. Nxp = ' + str(Nxp))
print('. len(y) = ' + str(len(y))) # = Nxp
print('. len(x) =', str(len(x)))
print('. Ndown =', str(Ndown))
print('. Nx =', str(Nx))
print('. Nxp =', str(Nxp))
print('. len(y) =', str(len(y))) # = Nxp
print('')
return y
def downsample_bpf(x, Ndown, k, Ndft, coefs, verbosity=1):
"""BPF at bin k of Ndft and downsample x by factor D = Ndown [HARRIS Fig 6.14]
Implement downsampling down converter for one bin [HARRIS Fig 6.14].
. see downsample()
Input:
. x: Input signal x[n]
. Ndown: downsample (decimation) factor
. k: Index of BPF center frequency w_k = 2 pi k / Ndft
. Ndft: number of frequency bins, DFT size
. coefs: prototype FIR filter coefficients for anti aliasing BPF
- verbosity: when > 0 print() status, else no print()
Return:
. y: Downsampled and down converted output signal y[m], m = n // D for bin
k. Complex baseband signal.
"""
a = [1.0] # FIR b = coefs, a = 1
# Polyphase implementation to avoid calculating values that are removed
# coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
# polyCoefs = [[ 3, 7, 11],
# [ 2, 6, 10],
# [ 1, 5, 9],
# [ 0, 4, 8]])
pfs = PolyPhaseFirFilterStructure(Ndown, coefs, commutator='down')
polyCoefs = pfs.polyCoefs
# Poly phases for whole data signal x, prepended with Ndown - 1 zeros
polyX, Nx, Nxp = poly_structure_data_for_downsampling_whole_x(x, Ndown) # size (Ndown, Nxp), Ny = Nxp
# Filter x per polyphase, order in polyCoefs accounts for commutator [HARRIS Fig 6.12, 6.13]
polyY = np.zeros((Ndown, Nxp))
for p in range(Ndown):
polyY[p] = signal.lfilter(polyCoefs[p], a, polyX[p])
# Phase rotate per polyphase, due to delay line at branch inputs [HARRIS Eq 6.8]
polyYC = np.zeros((Ndown, Nxp), dtype='cfloat')
for p in range(Ndown):
polyYC[p] = polyY[p] * np.exp(1j * 2 * np.pi * (Ndown - 1 - p) * k / Ndown)
# Sum the branch outputs to get single downsampled and downconverted output value
y = np.sum(polyYC, axis=0)
if verbosity:
print('> downsample_bpf():')
print('. len(x) =', str(len(x)))
print('. Ndown =', str(Ndown))
print('. Nx =', str(Nx))
print('. Nxp =', str(Nxp))
print('. len(y) =', str(len(y))) # = Nxp
print('. Ndft =', str(Ndft))
print('. k =', str(k))
print('')
return y
def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate and decimate by Nup / Ndown
"""Resample x by factor I / D = Nup / Ndown
"""Resample x by factor U / D = Nup / Ndown
x[n] --> Nup --> v[m] --> LPF --> w[m] --> Ndown --> y[k]
......@@ -401,7 +537,7 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate an
The resampling is done by first upsampling and then downsampling, because then only one shareed LPF is needed.
For upsampling an LPF is always needed, because it has to construct the inserted Nup - 1 zero values. For
downsampling the LPF of the upsampling already has restricted the input band width (BW), provided that the LPF
has BW < fNyquist / I and BW < fNyquist / D. The downsampling therefore does not need an LPF and reduces to
has BW < fNyquist / U and BW < fNyquist / D. The downsampling therefore does not need an LPF and reduces to
selecting appropriate output phase of the upsampler. The upsampler then only needs to calculate the output phase
that will be used by the downsampler, instead of calculating all output phases.
......@@ -416,13 +552,12 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate an
Return:
. y: Resampled output signal y[m]
Assumptions:
. x[n] = 0 for n < 0
. k = m // D = (n * I) // D, so first resampled output sample starts at k = 0 based on n = 0
. k = m // D = (n * U) // D, so first resampled output sample starts at k = 0 based on n = 0
. no y[k] for k < 0
. insert I - 1 zeros after each x[n] to get v[m], so len(v) = len(y) = I * len(x)
. use coefs as anti aliasing filter, must be LPF with BW < fNyquist / I and BW < fNyquist / D
. insert U - 1 zeros after each x[n] to get v[m], so len(v) = len(y) = U * len(x)
. use coefs as anti aliasing filter, must be LPF with BW < fNyquist / U and BW < fNyquist / D
. len(coefs) typically is multiple of Nup. If shorter, then the coefs are extended with zeros.
Remarks:
......@@ -430,17 +565,16 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate an
"""
Nx = len(x)
Nv = Nup * Nx
a = [1.0]
Nyp = (Nv + Ndown - 1) // Ndown # Number of samples from v, per poly phase in polyY
Ny = 1 + Ndown * (Nyp - 1) # Used number of samples from v
# Upsampling with LPF
w = upsample(x, Nup, coefs, verify=False)
# Downsampling selector
# This model is not efificent, because the upsample() has calculated all Nup phases, whereas it only needs to
# This model is not efficient, because the upsample() has calculated all Nup phases, whereas it only needs to
# calculate the phase that is selected by the downsampling. However, it is considered not worth the effort to
# make the model more efficient, because then only parts of lfilter() in upsample() need to be calulated, so
# then lfilter() needs to be implemented manually and possibly in a for loop {section 3.3.4 HARRIS].
# then lfilter() needs to be implemented manually and possibly in a for loop [section 3.3.4 CROCHIERE].
downSelect = np.arange(0, Nyp) * Ndown
y = w[downSelect]
......@@ -451,6 +585,7 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate an
v[0] = x
v = v.T.reshape(1, Nup * Nx)[0] # upsample
# . LPF
a = [1.0] # FIR b = coefs, a = 1
w = signal.lfilter(coefs, a, v)
# . Downsampling with calculating values that are removed
yVerify = np.zeros(Ndown * Nyp)
......@@ -465,11 +600,11 @@ def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1): # interpolate an
if verbosity:
print('> resample():')
print('. len(x) = ' + str(len(x)))
print('. Nx = ' + str(Nx))
print('. len(v) = ' + str(len(v)))
print('. Ny = ' + str(Ny))
print('. Nyp = ' + str(Nyp))
print('. len(y) = ' + str(len(y)))
print('. len(x) =', str(len(x)))
print('. Nx =', str(Nx))
print('. len(v) =', str(len(v)))
print('. Ny =', str(Ny))
print('. Nyp =', str(Nyp))
print('. len(y) =', str(len(y)))
print('')
return y
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment