multirate.py

#! /usr/bin/env python3
###############################################################################
#
# Copyright 2024
# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

###############################################################################

# Author: Eric Kooistra
# 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 sys import exit
from scipy import signal
from .utilities import c_rtol, c_atol, ceil_div

c_firA = [1.0]  # FIR b = coefs, a = 1


def zeros(shape, cmplx=False):
    if cmplx:
        return np.zeros(shape, dtype='cfloat')
    else:
        return np.zeros(shape)


def ones(shape, cmplx=False):
    return zeros(shape, cmplx) + 1


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 = zeros(len(x) * U, np.iscomplexobj(x))
    xU[0::U] = x
    return xU


###############################################################################
# Polyphase filter structure for input block processing
###############################################################################

class PolyPhaseFirFilterStructure:
    """Polyphase FIR filter structure (PFS) per block of data

    Purpose of PFS implementation is to avoid multiply by zero values in case
    of upsampling and to avoid calculating samples that will be discarded in
    case of downsampling. The output result of the PFS FIR is identical to
    using the FIR filter. The spectral aliasing and spectral replication are
    due to the sample rate change, not to the implementation structure.

    The FIR data can be real to process real data, or complex to process
    equivalent baseband data. The FIR coefficients are real.

    The PFS is is suitable for downsampling and for upsampling:
    - Upsampling uses the PFS as Transposed Direct-Form FIR filter, where the
      commutator selects the polyphases from phase 0 to Nphases - 1 to yield
      output samples for every input sample [HARRIS Fig 7.11, 7.12, CROCHIERE
      Fig 3.25].
    - Downsampling uses the PFS as a Direct-Form FIR filter, where the
      commutator selects the polyphases from phase Nphases - 1 downto 0 to sum
      samples from the past for one output sample [HARRIS Fig 6.12, 6.13,
      CROCHIERE Fig 3.29].
    - In both cases the commutator rotates counter clockwise, because the PFS
      uses a delay line of z^(-1) called type 1 structure [CROCHIERE 3.3.3,
      VAIDYANATHAN 4.3].

    Input:
    . coefs: Real FIR coefficients b[0 : Nphases * Ntaps - 1].
    . Nphases: number of polyphases, is number of rows (axis=0) in the
               polyphase structure.
    . cmplx: Define PFS delays for real (when False) or complex (when True)
             FIR data.
    Derived:
    . Ntaps: number of taps per polyphase, is number of columns (axis=1) in
      the polyphase structure.
    . polyCoefs: FIR coefficients storage with Nphases rows and Ntaps columns.
    . polyDelays: FIR taps storage with Nphases rows and Ntaps columns.

    Storage of FIR coefficients and FIR data:
      Stream of input samples x[n], n >= 0 for increasing time. Show index n
      also as value so x[n] = n:

                 oldest sample                             newest sample
        x         v                                                  v
        samples  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...]
        blocks   [         0,          1,            2,              3]
                                                                     ^
                                                           newest block

      The delayLine index corresponds to the FIR coefficient index k of
      b[k] = impulse response h[k]. For Ncoefs = 12, with Nphases = 4 rows and
      Ntaps = 3 columns:

        polyCoefs = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11] = b[k]
                                 row = polyphase
        polyCoefs = [[0, 4, 8],    0
                     [1, 5, 9],    1
                     [2, 6,10],    2
                     [3, 7,11]]    3
             column:  0, 1, 2

      Shift x[n] samples into delayLine from left, show sample index n in (n)
      and delayLine index k in [k]. Shift in per sample or per block.

        shift in -->(15,14,13,12,11,10, 9, 8, 7, 6, 5, 4) --> shift out
        delayLine = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11]

      For polyDelays shift in from left-top (n = 15 is newest input sample)
      and shift out from right-bottom (n = 4 is oldest sample).

                       v              v
        polyDelays = [[0, 4, 8],   ((15,11, 7),
                      [1, 5, 9],    (14,10, 6),
                      [2, 6,10],    (13, 9, 5),
                      [3, 7,11]]    (12, 8, 4))
                             v              v

      Output sample y[n] = sum(polyCoefs * polyDelays).

      For downsampling by Ndown = Nphases shift blocks of Ndown samples in
      from left in polyDelays, put newest block[3] = x[12,13,14,15] in first
      column [0] of polyDelays.

      Store input samples in polyDelays structure, use mapping to access as
      delayLine:
        . map_to_delay_line()
        . map_to_poly_delays()
    """
    def __init__(self, Nphases, coefs, cmplx=False):
        self.coefs = coefs
        self.Ncoefs = len(coefs)
        self.Nphases = Nphases  # branches, rows
        self.Ntaps = ceil_div(self.Ncoefs, Nphases)  # taps, columns
        self.cmplx = cmplx
        self.init_poly_coeffs()
        self.reset_poly_delays()

    def init_poly_coeffs(self):
        """Polyphase structure for FIR filter coefficients coefs in Nphases.

        Input:
        . coefs = prototype FIR filter coefficients (= b[k] = impulse response)
        Return:
        . If Ncoefs < Nphases * Ntaps, then pad coefs with zeros.
        . E.g. coefs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
               coefs.reshape((4, 3)) = [[0, 1, 2],
                                        [3, 4, 5],
                                        [6, 7, 8],
                                        [9,10,11])
        . Therefore for Nphases = 4 rows and Ntaps = 3 columns use:
          polyCoefs = coefs.reshape((3, 4)).T
                    = [[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)
        """
        polyCoefs = np.zeros(self.Nphases * self.Ntaps)  # real
        polyCoefs[0 : self.Ncoefs] = self.coefs
        self.polyCoefs = polyCoefs.reshape((self.Ntaps, self.Nphases)).T

    def reset_poly_delays(self):
        self.polyDelays = zeros((self.Nphases, self.Ntaps), self.cmplx)

    def map_to_delay_line(self):
        delayLine = self.polyDelays.T.reshape(-1)
        return delayLine

    def map_to_poly_delays(self, delayLine):
        self.polyDelays = delayLine.reshape((self.Ntaps, self.Nphases)).T

    def shift_in_data(self, inData, flipped):
        """Shift block of data into the polyDelays structure.

        View polyDelays as delay line if L < Nphases. Shift in from the left
        in the delay line, so last (= newest) sample in inData will appear at
        most left of the delay line. If L = Nphases, then it is possible to
        shift the data directly into the first (= left) column of polyDelays.

        Input:
        . inData: Block of one or more input samples with index as time index
            n in inData[n], so oldest sample at index 0 and newest sample at
            index -1.
        . flipped: False then inData order is inData[n] and still needs to be
            flipped. If True then the inData is already flipped.
        """
        L = len(inData)
        xData = inData if flipped else np.flip(inData)
        if L < self.Nphases:
            delayLine = self.map_to_delay_line()
            # Equivalent code:
            # delayLine = np.concatenate((inData, delayLine[L:]))
            delayLine = np.roll(delayLine, L)
            delayLine[:L] = xData
            self.map_to_poly_delays(delayLine)
        else:
            # Equivalent code for L == Nphases: Shift in inData block directly
            # at column 0
            self.polyDelays = np.roll(self.polyDelays, 1, axis=1)
            self.polyDelays[:, 0] = xData

    def filter_block(self, inData, flipped):
        """Filter block of inData per polyphase.

        Input:
        . inData: block of input data with length <= Nphases, time index n in
            inData[n] increments, so inData[0] is oldest data and inData[-1] is
            newest data. The inData block size is = 1 for upsampling or > 1
            and <= Nphases for downsampling.
        . flipped: False then inData order is inData[n] with newest sample last
            and still needs to be flipped. If True then the inData is already
            flipped in time, so with newest sample first.
        Return:
        . pfsData: block of polyphase FIR filtered output data for Nphases, with
            pfsData[p] and p = 0:Nphases-1 from top to bottom.
        """
        # Shift in one block of input data (1 <= len(inData) <= Nphases)
        self.shift_in_data(inData, flipped)
        # Apply FIR coefs per delay element
        zData = self.polyDelays * self.polyCoefs
        # Sum FIR taps per polyphase
        pfsData = np.sum(zData, axis=1)
        # Output block of Nphases filtered output data
        return pfsData

    def filter_up(self, inSample):
        """Filter same input sample at each polyphase, to upsample it by factor
        U = Nphases.

        Input:
        . inSample: input sample that is applied to each polyphase
        Return:
        . pfsData: block of polyphase FIR filtered output data for Nphases,
            with pfsData[p] and p = 0:Nphases-1 from top to bottom. The
            pfsData[0] is the first and thus oldest and pfsData[-1] is the last
            and thus newest interpolated data.
        """
        inWires = inSample * ones(self.Nphases, np.iscomplexobj(inSample))
        pfsData = self.filter_block(inWires, flipped=True)
        return pfsData


###############################################################################
# Up, down, resample low pass input signal
###############################################################################

def polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros):
    """Polyphase data structure for downsampling whole signal x.

    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 Nzeros = 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 to yield a
      new output FIR sum.

    Input:
    . x: all input samples x[n] for n = 0 : Lx - 1
    . Ndown: downsample rate and number of polyphase branches
    . Nzeros: number of zero samples to prepend for x
    Return:
    . polyX: polyphase data structure with size (Ndown, Nxp) for Ndown branches
        and Nxp samples from x per branch.
    . Nx: Total number of samples from x, including prepended Ndown - 1 zeros.
    . Nxp: Total number of samples used from x per polyphase branch, is the
        number of samples Ny, that will be in downsampled output y[m] = x[m D],
        for m = 0, 1, 2, ..., Nxp - 1.
    """
    Lx = len(x)
    Nxp = (Nzeros - 1 + Lx) // Ndown  # Number of samples per polyphase
    Nx = 1 + Ndown * (Nxp - 1)  # Used number of samples from x

    # Load x into polyX with Ndown rows = polyphases
    # . prepend x with Ndown - 1 zeros
    # . skip any last remaining samples from x, that are not enough yield a new
    #   output FIR sum.
    polyX = zeros(Ndown * Nxp, np.iscomplexobj(x))
    polyX[Ndown - 1] = x[0]
    polyX[Ndown:] = x[1 : Nx]
    # . Store data in time order per branch, so with oldest data left, to match
    #   use with lfilter(). Note this differs from order of tap data in
    #   pfs.polyDelays where newest data is left, because the block data is
    #   flipped by shift_in_data(), to match use in pfs.filter_block(), so that
    #   it can use pfs.polyDelays * pfs.polyCoefs to filter the block.
    # . The newest sample x[-1] has to be in top branch of polyX and the oldest
    #   sample x[0] in bottom branch, to match branch order of 0:Ndown-1 from
    #   top to bottom in the pfs.polyCoefs, therefore do flipud(polyX). This
    #   flipud() is similar as the flip() per block in shift_in_data().
    #
    #         x[0] is oldest                              x[-1] is newest
    #           |                                            |
    #   x[n] = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15]
    #          (0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15), x[n] = n
    #                                                              newest
    #   poly                 poly        v                             |
    #   Coefs = [[0, 4, 8],  Delays = ((15,11, 7),  polyX = ((3, 7,11,15),
    #            [1, 5, 9],            (14,10, 6),           (2, 6,10,14),
    #            [2, 6,10],            (13, 9, 5),           (1, 5, 9,13),
    #            [3, 7,11]]            (12, 8, 4))           (0, 4, 8,12))
    #                                          v              |
    #                                                       oldest
    polyX = np.flipud(polyX.reshape(Nxp, Ndown).T)
    return polyX, Nx, Nxp


def polyphase_frontend(x, Nphases, coefs, sampling):
    """Calculate PFS FIR filter output of x for Nphases.

    The polyY[Nphases] output of this PFS frontend can be:
    . summed for dowsampling by factor D = Nphases
    . used as is for upsampling by factor U = Nphases
    . used with a range of Nphases LOs for a single channel downsampler
      downconverter, or upsampler upconverter
    . used with a IDFT for a analysis filterbank (PFB), or synthesis PFB.

    Input:
    . x: Input signal x[n]
    . Nphases: number of polyphase branches in PFS
    . coefs: prototype FIR filter coefficients for anti aliasing LPF. The
        len(coefs) typically is multiple of Nphases. If shorter, then the
        coefs are extended with zeros.
    . sampling: 'up' or 'down'
    Return:
    . polyY[Nphases]: Output of the FIR filtered branches in the PFS.
    . Nx = np.size(polyY), total number of samples from x
    . Nxp = np.size(polyY, axis=1), total number of samples from x per
        polyphase.
    """
    # Use polyphase FIR filter coefficients from pfs class.
    pfs = PolyPhaseFirFilterStructure(Nphases, coefs, np.iscomplexobj(x))

    # Define polyphases for whole data signal x with Nx samples, instead of
    # using polyDelays per Ntaps from pfs class, to be able to use
    # signal.lfilter() per polyphase branch for the whole data signal x.
    if sampling == 'up':
        Nx = len(x)
        Nxp = Nx
        Nup = Nphases
        # Filter whole x per polyphase, so Nxp = Nx, because the Nup branches
        # each use all x to yield Nup output samples per input sample from x.
        # The commutator index order for upsampling is p = 0, 1, .., Nup - 1,
        # so from top to bottom in the PFS.
        polyY = zeros((Nup, Nx), np.iscomplexobj(x))

        pCommutator = range(Nup)
        for p in pCommutator:
            polyY[p] = Nup * signal.lfilter(pfs.polyCoefs[p], c_firA, x)
    else:
        # 'down':
        # . Prepend x with Ndown - 1 zeros, to have y[m] for m = 0 start at
        #   n = 0 of x[n].
        # . Size of polyX is (Ndown, Nxp), and length of y is Ny = Nxp is
        #   length of each branch.
        Ndown = Nphases
        Nzeros = Ndown - 1
        polyX, Nx, Nxp = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros)
        # print(polyX[:, 0])
        # Filter Ndown parts of x per polyphase, because the FIR filter output
        # y will sum. The commutator index order for downsampling is p =
        # Ndown - 1,..., 1, 0, so from bottom to top in the PFS. However, the
        # commutator index order is only relevant when the branches are
        # calculated sequentially to reuse the same hardware, because for the
        # output y the branches are summed anyway.
        polyY = zeros((Ndown, Nxp), np.iscomplexobj(x))
        pCommutator = np.flip(np.arange(Ndown))
        for p in pCommutator:
            polyY[p] = signal.lfilter(pfs.polyCoefs[p], c_firA, polyX[p])
    return polyY, Nx, Nxp


def upsample(x, Nup, coefs, verify=False, verbosity=1):  # interpolate
    """Upsample x by factor U = Nup and LPF.

    Input:
    . x: Input signal x[n]
    . Nup: upsample (interpolation) factor
    . 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()
    Return:
    . y: Upsampled output signal y[m]

    Assumptions:
    . x[n] = 0 for 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 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 / 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 * U) == 0, then the group delay is an integer
        number of ts periods.
    """
    # Polyphase FIR filter input x
    polyY, Nx, _ = polyphase_frontend(x, Nup, coefs, 'up')

    # Output Nup samples y[m] for every input sample x[n]
    y = polyY.T.reshape(1, Nup * Nx)[0]

    if verify:
        # Inefficient upsampling implementation with multiply by zero values.
        # Verify that x --> U --> LPF --> y yields identical result y as with
        # using the PFS: x --> PFS FIR --> y.
        xZeros = zeros((Nup, Nx), np.iscomplexobj(x))
        xZeros[0] = x
        xZeros = xZeros.T.reshape(1, Nup * Nx)[0]  # upsample
        yVerify = Nup * signal.lfilter(coefs, c_firA, xZeros)  # LPF
        print('> Verify upsample():')
        if np.allclose(y, yVerify, rtol=c_rtol, atol=c_atol):
            print('  . PASSED: correct upsample result')
        else:
            exit('ERROR: wrong upsample result')
    if verbosity:
        print('> Log upsample():')
        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
    """LPF and downsample x by factor D = Ndown.

    Input:
    . x: Input signal x[n]
    . Ndown: downsample (decimation) factor
    . 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()
    Return:
    . y: Downsampled output signal y[m]

    Assumptions:
    . x[n] = 0 for 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
    . len(coefs) typically is multiple of Ndown. If shorter, then the coefs
      are extended with zeros.

    Remarks:
    . The input sample period is ts and the output sample period of the
      downsampled (= decimated signal) is tsDown = ts * D
    . The group delay is (Ncoefs - 1) / 2 * ts. With odd Ncoefs and symmetrical
      coefs to have linear phase, the group delay is an integer number of ts
      periods, so:
      - when (Ncoefs - 1) %  2      == 0, then the group delay is an integer
        number of ts periods
      - when (Ncoefs - 1) % (2 * D) == 0, then the group delay is an integer
        number of tsDown periods
    """
    # Polyphase FIR filter input x
    polyY, Nx, Nxp = polyphase_frontend(x, Ndown, coefs, 'down')

    # FIR filter sum
    y = np.sum(polyY, axis=0)

    if verify:
        # Inefficient downsampling implementation with calculating values that
        # are removed. Verify that x --> LPF --> D --> y yields identical
        # result y as with using the PFS: x --> PFS FIR --> y.
        yVerify = zeros(Ndown * Nxp, np.iscomplexobj(x))
        yVerify[0 : Nx] = signal.lfilter(coefs, c_firA, x[0 : Nx])  # LPF
        yVerify = yVerify.reshape(Nxp, Ndown).T   # polyphases
        yVerify = yVerify[0]   # downsample by D
        print('> Verify downsample():')
        if np.allclose(y, yVerify, rtol=c_rtol, atol=c_atol):
            print('  . PASSED: correct downsample result')
        else:
            exit('ERROR: wrong downsample result')
    if verbosity:
        print('> Log 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('')
    return y


def resample(x, Nup, Ndown, coefs, verify=False, verbosity=1):  # interpolate and decimate by Nup / Ndown
    """Resample x by factor U / D = Nup / Ndown

    x[n] --> Nup --> v[m] --> LPF --> w[m] --> Ndown --> y[k]

    The w[m] is the result of upsampling x[n] by Nup. The LPF has narrowed the
    pass band to suit both Nup and Ndown. Therefore the downsampling by Ndown
    implies that only one in every Ndown samples of w[m] needs to be calculated
    to get y[k] = w[k * D]. Hence resampling by Nup / Ndown takes less
    processing than upsampling when Ndown > 1.

    Poly phases in up and downsampling:
    The upsampling structure outputs the interpolated samples via the separate
    output phases. The same input sequence x is used to calculate all Nup
    phases, whereby each phase uses a different set of coefficients from the
    LPF. The downsampling structure sums the Ndown phases to output the
    downsampled samples. Each phase uses a different set of coefficients from
    the LPF to filter Ndown delay phases of the input sequence x.

    Resampling is upsampling with downsampling by phase selection:
    The resampling is done by first upsampling and then downsampling, because
    then only one shared 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 / 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.

    Input:
    . x: Input signal x[n]
    . Nup: upsample (interpolation) factor
    . Ndown: downsample (decimation) factor
    . coefs: 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()
    Return:
    . y: Resampled output signal y[m]

    Assumptions:
    . x[n] = 0 for 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 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:
    . The group delay is the same as for the upsampler, because the
      downsampling has no additional filtering.
    """
    Nx = len(x)
    Nv = Nup * Nx
    # Nyp is number of samples from v, per polyphase in polyY
    Nyp = (Nv + Ndown - 1) // Ndown
    # Ny is used number of samples from v
    Ny = 1 + Ndown * (Nyp - 1)

    # Upsampling with LPF
    w = upsample(x, Nup, coefs, verify=False)
    # Downsampling selector
    # 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 [CROCHIERE 3.3.4].
    downSelect = np.arange(0, Nyp) * Ndown
    y = w[downSelect]

    if verify:
        # Inefficient implementation
        # . Upsampling with multiply by zero values
        v = np.zeros((Nup, Nx))   # polyphases
        v[0] = x
        v = v.T.reshape(1, Nup * Nx)[0]  # upsample
        # . LPF
        w = Nup * signal.lfilter(coefs, c_firA, v)
        # . Downsampling with calculating values that are removed
        yVerify = np.zeros(Ndown * Nyp)
        yVerify[0 : Ny] = w[0 : Ny]
        yVerify = yVerify.reshape(Nyp, Ndown).T   # polyphases
        yVerify = yVerify[0]   # downsample

        print('> Verify resample():')
        if np.allclose(y, yVerify, rtol=c_rtol, atol=c_atol):
            print('  . PASSED: correct resample result')
        else:
            exit('ERROR: wrong resample result')
    if verbosity:
        print('> Log 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('')
    return y


###############################################################################
# Single bandpass channel up and downsampling and up and downconversion
###############################################################################

def unit_circle_loops_phasor_arr(k, N, sign):
    """Return array of N phasors on k loops along the unit circle.

    Polyphase dependent phase offsets for bin k, when N = Ndown = Ndft [HARRIS
    Eq 6.8]. For k = 1 this yields the roots of unity.

    Input:
    . k: bin in range(N)
    . N: number of phasors
    . sign: +1 or -1 term in exp()
    Return:
    . pArr: exp(sign 2j pi k / N) for k in 0 : N - 1
    """
    pArr = np.array([np.exp(sign * 2j * np.pi * p * k / N) for p in range(N)])
    return pArr


def time_shift_phasor_arr(k, M, Ndft, Msamples, sign):
    """Return array of Msamples phasors in time to compensate for oversampling
    time shift.

    The time shift due to down or upsampling causes a frequency component of
    k * M / Ndft. With oversampling M < Ndft, and then after down or upsampling
    there remains a frequency offset that can be compensated by given by mArr
    [HARRIS Eq 9.3].

    Input:
    . k: Index of BPF center frequency w_k = 2 pi k / Ndft
    . M: downsample or upsample factor
    . Ndft: DFT size, number of bins and number of polyphases in PFS FIR filter
    . Msamples: Requested number of output samples in time
    . sign: +1 or -1 term in exp()
    Return:
    . mArr = exp(sign 2j pi k * M / Ndft * m),  for m in 0 : Msamples - 1
    """
    mArr = np.exp(sign * 2j * np.pi * k * M / Ndft * np.arange(Msamples))
    return mArr


def maximal_downsample_bpf(x, Ndown, k, coefs, verbosity=1):
    """BPF x at bin k in range(Ndown) and downsample x by factor D = Ndown and
    Ndown = Ndft.

    Implement maximal downsampling downconverter for one bin (= critically
    sampled) [HARRIS Fig 6.14].

    The BPF is centered at w_k = 2pi k / Ndft, where Ndft is number of
    frequency bins, is DFT size. The downsampling is maximal so Ndown = Ndft.
    The polyphase structure has Nphases = Ndown branches, so the input x
    data that shifts in remains in each branch. Therefore each branch can be
    FIR filtered independently for the whole input x using
    polyphase_frontend().

    Input:
    . x: Input signal x[n]
    . Ndown: downsample factor
    . k: Index of BPF center frequency w_k = 2 pi k / Ndown
    . coefs: prototype FIR filter coefficients for anti aliasing BPF
    - verbosity: when > 0 print() status, else no print()
    Return:
    . yc: Downsampled and downconverted output signal yc[m], m = n // D for
          bin k. Complex baseband signal.
    """
    # Polyphase FIR filter input x
    polyY, Nx, Nxp = polyphase_frontend(x, Ndown, coefs, 'down')

    # Phase rotate per polyphase for bin k, due to delay line at branch inputs
    # [HARRIS Eq 6.8]
    kPhasors = unit_circle_loops_phasor_arr(k, Ndown, 1)
    polyYc = np.zeros((Ndown, Nxp), dtype='cfloat')
    for p in range(Ndown):
        polyYc[p] = polyY[p] * kPhasors[p]  # row = row * scalar

    # Sum the branch outputs to get single downsampled and downconverted output
    # complex baseband value yc.
    yc = np.sum(polyYc, axis=0)

    if verbosity:
        print('> Log maximal_downsample_bpf():')
        print('  . len(x)  =', str(len(x)))
        print('  . Nx      =', str(Nx))
        print('  . Nxp     =', str(Nxp))
        print('  . len(yc) =', str(len(yc)))  # = Nxp
        print('  . Ndown   =', str(Ndown))
        print('  . k       =', str(k))
        print('')
    return yc


def maximal_upsample_bpf(xBase, Nup, k, coefs, verbosity=1):
    """BPF xBase at bin k in range(Ndft), and upsample xBase by factor U = Nup =
    Ndft.

    Implement maximal upsampling upconverter for one bin (= critically
    sampled) [HARRIS Fig. 7.16].

    The BPF is centered at w_k = 2pi k / Ndft, where Ndft is number of
    frequency bins, is DFT size. The upsampling is maximal so Nup = Ndft. The
    polyphase structure has Nphases = Nup branches, so the output yc data
    shifts out from the same branch for each block of Nup samples. Therefore
    each branch can be FIR filtered independently for the whole input xBase
    using polyphase_frontend().

    Input:
    . xBase: Input equivalent baseband signal
    . Nup: upsample factor
    . k: Index of BPF center frequency w_k = 2 pi k / Nup
    . coefs: prototype FIR filter coefficients for anti aliasing and
        interpolating BPF
    - verbosity: when > 0 print() status, else no print()
    Return:
    . yc: Upsampled and upconverted output signal yc[n], n = m U at
          intermediate frequency (IF) of bin k. Complex positive frequencies
          only signal.
    """
    # Polyphase FIR filter input xBase
    polyY, Nx, Nxp = polyphase_frontend(xBase, Nup, coefs, 'up')

    # Phase rotate per polyphase for bin k, due to delay line at branch inputs
    # [HARRIS Eq 7.8 = 6.8, Fig 7.16], can be applied after the FIR filter,
    # because the kPhasors per polyphase are constants.
    kPhasors = unit_circle_loops_phasor_arr(k, Nup, 1)
    polyYc = np.zeros((Nup, Nxp), dtype='cfloat')
    for p in range(Nup):
        polyYc[p] = polyY[p] * kPhasors[p]  # row = row * scalar

    # Output Nup samples yc[m] for every input sample xBase[n]
    yc = polyYc.T.reshape(1, Nup * Nx)[0]

    if verbosity:
        print('> Log maximal_upsample_bpf():')
        print('  . len(xBase) =', str(len(xBase)))
        print('  . Nx         =', str(Nx))
        print('  . Nxp        =', str(Nxp))
        print('  . len(yc)    =', str(len(yc)))  # = Nxp
        print('  . Nup        =', str(Nup))
        print('  . k          =', str(k))
        print('')
    return yc


def non_maximal_downsample_bpf(x, Ndown, k, Ndft, coefs, verbosity=1):
    """BPF x at bin k in range(Ndown), and downsample x by factor D = Ndown.

    Implement nonmaximal downsampling downconverter for one bin. The maximal
    downsampling downconverter for one bin has the kPhasors per polyphase
    branch of [HARRIS Eq. 6.8 and Fig 6.14]. For nonmaximal this needs to be
    extended with the tPhasors per polyphase branch of [HARRIS Eq. 9.3, to
    compensate for the oversampling time shift.

    The BPF is centered at w_k = 2pi k / Ndft, where Ndft is number of
    frequency bins, is DFT size. The polyphase FIR structure has Nphases = Ndft
    branches, to fit the requested number of bins. The polyphase FIR structure
    is maximally downsampled (= critically sampled) for Ndown = Ndft, but it
    can support any Ndown <= Ndft. The input data shifts in per Ndown samples,
    so it appears in different branches when Ndown < Ndft and a new block is
    shifted in. Therefore the input data cannot be FIR filtered per branch for
    the whole input x. Instead it needs to be FIR filtered per block of Ndown
    input samples from x, using pfs.polyDelays in pfs.filter_block().

    Input:
    . x: Input signal x[n]
    . Ndown: downsample factor
    . k: Index of BPF center frequency w_k = 2 pi k / Ndft
    . Ndft: DFT size, number of polyphases in PFS FIR filter
    . coefs: prototype LPF FIR filter coefficients for anti aliasing BPF
    - verbosity: when > 0 print() status, else no print()
    Return:
    . yc: Downsampled and downconverted output signal yc[m], m = n // D for
          bin k. Complex baseband signal.
    """
    # Prepend x with Ndown - 1 zeros, and represent x in Nblocks of Ndown
    # samples
    Nzeros = Ndown - 1
    xBlocks, Nx, Nblocks = polyphase_data_for_downsampling_whole_x(x, Ndown, Nzeros)
    # print(xBlocks[:, 0])

    # Prepare output
    yc = np.zeros(Nblocks, dtype='cfloat')

    # PFS with Ndft polyphases
    pfs = PolyPhaseFirFilterStructure(Ndft, coefs)
    kPhasors = unit_circle_loops_phasor_arr(k, Ndft, 1)  # [HARRIS Eq 6.8]

    # Oversampling time shift compensation via frequency dependent phase shift
    tPhasors = time_shift_phasor_arr(k, Ndown, Ndft, Nblocks, -1)  # [HARRIS Eq 9.3]

    for b in range(Nblocks):
        # Filter block
        inData = xBlocks[:, b]
        pfsData = pfs.filter_block(inData, flipped=True)
        # Phase rotate polyphases for bin k
        pfsBinData = pfsData * kPhasors  # [HARRIS Eq 6.8, 9.3]
        # Sum the polyphases to get single downsampled and downconverted output
        # value
        yc[b] = np.sum(pfsBinData) * tPhasors[b]

    if verbosity:
        print('> non_maximal_downsample_bpf():')
        print('  . len(x)   =', str(len(x)))
        print('  . Nx       =', str(Nx))
        print('  . Nblocks  =', str(Nblocks))
        print('  . len(yc)  =', str(len(yc)))  # = Nblocks
        print('  . Ndown    =', str(Ndown))
        print('  . Ndft     =', str(Ndft))
        print('  . k        =', str(k))
        print('')
    return yc


def non_maximal_upsample_bpf(xBase, Nup, k, Ndft, coefs, verbosity=1):
    """BPF xBase at bin k in range(Nup), and upsample xBase by factor U = Nup.

    Implement nonmaximal upsampling upconverter for one bin. The maximal
    upsampling upconverter for one bin has the kPhasors per polyphase
    branch similar of [HARRIS Fig 7.16]. For nonmaximal this needs to be
    extended with the tPhasors per polyphase branch similar as for down in
    [HARRIS Eq. 9.3], to compensate for the oversampling time shift.

    The BPF is centered at w_k = 2pi k / Ndft, where Ndft is number of
    frequency bins, is DFT size. The polyphase FIR structure has Nphases = Ndft
    branches, to fit the requested number of bins. The polyphase FIR structure
    is maximally upsampled (= critically sampled) for Nup = Ndft, but it
    can support any Nup <= Ndft. The output data shifts out per Nup samples,
    so it appears from different branches when Nup < Ndft and a new block is
    shifted out. Therefore the output data cannot be FIR filtered per branch
    for the whole input xBase. Instead it needs to be FIR filtered per block of
    Nup output samples, using pfs.polyDelays in pfs.filter_up().

    TODO
    . This code only runs ok for Ros = 1, 2 when Ndft = 16, but not for
      Ros = 4, so need to fix that. According to [CROCHIERE 7.2.7] the
      polyphase structure is only suitable for Ros is integer >= 1. For other
      Ros > 1 the weighted overlap-add (WOLA) structure is  suitable, so
      need to add WOLA.

    Input:
    . xBase: Input equivalent baseband signal xBase[m]
    . Nup: upsample factor
    . k: Index of BPF center frequency w_k = 2 pi k / Ndft
    . Ndft: DFT size, number of polyphases in PFS FIR filter
    . coefs: prototype LPF FIR filter coefficients for anti aliasing and
             interpolationBPF
    - verbosity: when > 0 print() status, else no print()
    Return:
    . yc: Upsampled and upconverted output signal yc[n], n = m U at
          intermediate frequency (IF) of bin k. Complex positive frequencies
          only signal.
    """
    Nblocks = len(xBase)

    # Prepare output
    polyYc = np.zeros((Nup, Nblocks), dtype='cfloat')

    # PFS with Ndft polyphases
    pfsUp = PolyPhaseFirFilterStructure(Nup, coefs, cmplx=True)
    kPhasors = unit_circle_loops_phasor_arr(k, Ndft, 1)  # [HARRIS Eq 7.8]

    # Oversampling time shift compensation via frequency dependent phase shift
    tPhasors = time_shift_phasor_arr(k, Nup, Ndft, Nblocks, -1)  # [HARRIS Eq 9.3]

    # Polyphase FIR filter input xBase
    for b in range(Nblocks):
        # Filter block
        inSample = xBase[b]
        pfsData = Nup * pfsUp.filter_up(inSample)
        # Phase rotate polyphases for bin k
        pfsBinData = pfsData * kPhasors[0:Nup]
        # Output Nup samples yc[n] for every input sample xBase[m]
        outBinData = pfsBinData * tPhasors[b]
        polyYc[:, b] = outBinData

    yc = polyYc.T.reshape(1, Nup * Nblocks)[0]

    if verbosity:
        print('> non_maximal_upsample_bpf():')
        print('  . len(xBase) =', str(len(xBase)))
        print('  . Nblocks    =', str(Nblocks))
        print('  . len(yc)    =', str(len(yc)))  # = Nblocks
        print('  . Nup        =', str(Nup))
        print('  . Ndft       =', str(Ndft))
        print('  . k          =', str(k))
        print('')
    return yc