symbolTimingSync.m

function [ xI ] = symbolTimingSync(TED, intpl, L, mfIn, mfOut, K1, K2, ...
    const, Ksym, rollOff, rcDelay, debug_s, debug_r)
% Symbol Timing Loop
% ---------------------
%
% Implements closed-loop symbol timing recovery with a configurable timing
% error detector (TED) and a configurable interpolator. The feedback
% control loop uses a proportional-plus-integrator (PI) controller and a
% modulo-1 counter to control the interpolator. Meanwhile, the TED can be
% configured from five alternative implementations:
%     - Maximum-likelihood TED (MLTED);
%     - Early-late TED (ELTED);
%     - Zero-crossing TED (ZCTED);
%     - Gardner TED (GTED);
%     - Mueller-Muller TED (MMTED).
% The interpolator can be chosen from four implementations:
%     - Polyphase;
%     - Linear;
%     - Quadratic;
%     - Cubic.
%
% When using a polyphase interpolator, the loop simultaneously synchronizes
% the symbol timing and implements the matched filter (MF). In this case,
% the loop processes the MF input sequence directly, and there is no need
% for an external MF block. In contrast, when using any other interpolation
% method (linear, quadratic, or cubic), the loop processes the MF output
% and produces interpolated values based on groups of MF output samples.
% For instance, the linear interpolator projects a line between a pair of
% MF output samples and produces an interpolant along this line. Thus, when
% using the linear, quadratic, or cubic interpolators, this symbol
% synchronizer loop must be preceded by a dedicated MF block.
%
% In any case, to support both pre-MF and post-MF interpolation approaches,
% this function takes both the MF input and MF output sequences as input
% arguments (i.e., the mfIn and mfOut arguments). The mfOut argument can be
% empty when using the polyphase interpolator, while the mfIn argument can
% be empty when using the other interpolation methods. The only exception
% is if using the MLTED. The MLTED uses the so-called derivative matched
% filter (dMF), which takes the MF input sequence and produces the
% differentiated MF output. Thus, when using the MLTED, even with the
% linear, quadratic, or cubic interpolators, the MF input sequence must be
% provided on the "mfIn" input argument.
%
% Note the reason why the matched filtering is executed jointly with symbol
% synchronization when using the polyphase interpolator is merely because
% it is possible to do so. The cascaded combination of the Tx root raised
% cosine (RRC) filter and the matched RRC filter results in a raised cosine
% filter, which, in turn, is an Lth-band filter (also known as Nyquist
% filter) that is adequate for interpolation (see [2]). Hence, the two
% tasks (interpolation and matched filtering) can be achieved in one go,
% which saves the need and computational cost of an extra dedicated MF
% block. If it were not for this approach, the computational cost of the
% polyphase interpolator would typically be higher than the other
% interpolation methods. In contrast, by using the polyphase interpolator
% jointly as the MF, its computational cost becomes nearly zero, since it
% only implements the indispensable MF computations.
%
% Furthermore, note all TED schemes except the GTED compute the symbol
% timing error using symbol decisions. Hence, the reference constellation
% and scaling factor must be provided through input arguments 'const' and
% 'Ksym'. Such TED schemes could also leverage prior knowledge and compute
% the timing error using known symbols instead of decisions. However, this
% function does not offer the data-aided alternative. Instead, it only
% implements the decision-directed flavor of each TED. Meanwhile, the
% GTED is the only scheme purely based on the raw input samples, which is
% not decision-directed nor data-aided. Hence, the 'const' and 'Ksym'
% arguments are irrelevant when using the GTED.
%
% Input Arguments:
% TED     -> TED scheme ('MLTED', 'ELTED', 'ZCTED', 'GTED', or 'MMTED').
% intpl   -> Interpolator: 0) Polyphase; 1) Linear; 2) Quadratic; 3) Cubic.
% L       -> Oversampling factor.
% mfIn    -> MF input sequence sampled at L samples/symbol.
% mfOut   -> MF output sequence sampled at L samples/symbol.
% K1      -> PI controller's proportional gain.
% K2      -> PI controller's integrator gain.
% const   -> Symbol constellation.
% Ksym    -> Symbol scaling factor to be undone prior to slicing.
% rollOff -> Matched filter's rolloff factor.
% rcDelay -> Raised cosine filter delay (double the MF RRC delay).
% debug_s -> Show static debug plots after the loop processing.
% debug_r -> Open scopes for real-time monitoring over the loop iterations.
%
% References:
%   [1] Michael Rice, Digital Communications - A Discrete-Time Approach.
%   New York: Prentice Hall, 2008.
%   [2] Milić, Ljiljana. Multirate Filtering for Digital Signal Processing:
%   MATLAB Applications. Information Science Reference, 2009.

if (nargin < 12)
    debug_s = 0;
end

if (nargin < 13)
    debug_r = 0;
end

% Midpoint between consecutive symbols
%
% Some of the TED schemes (ELTED, ZCTED, and GTED) rely on the interpolants
% located halfway between two consecutive output interpolants (or output
% symbols). For instance, the ZCTED computes the timing error using the
% zero-crossing value obtained from interpolation, referred to as the
% "zero-crossing interpolant". When processing the k-th strobe, the loop
% estimates the zero-crossing interpolant by applying the offset "mu(k)" on
% the sample located at the basepoint index "m(k) - L/2".
%
% However, the problem is that the midpoint offset at "m(k) - L/2" only
% works if L is even. If L is odd, we can take "ceil(L/2)" and compensate
% for the discrepancy using the fractional timing offset mu. For instance,
% for L=3, the basepoint index for the zero-crossing interpolant would be
% located at "m(k) - 2", and the fractional timing offset can be adjusted
% to "mu(k) + 0.5". Similarly, the ELTED's "early" value would be computed
% using the sample at "m(k) + 2" as the basepoint index and "mu(k) - 0.5"
% as the fractional symbol timing offset.
%
% Finally, note that "mu(k) +-0.5" may not fall within the [0, 1) range.
% For instance, if mu(k)=0, then "mu(k) - 0.5" would be negative. In this
% case, we can move the basepoint index and readjust mu(k). For example, if
% "mu(k) - 0.5" is negative, we can move the the basepoint index to the
% preceding sample and use "mu(k) - 0.5 + 1" instead. See the adjustment on
% the "interpolate()" function.
midpointOffset = ceil(L / 2);
muOffset = midpointOffset - L/2; % 0.5 if L is odd, 0 if L is even

% Modulation order
M = numel(const);

% Make sure the input vectors are column vectors
if (size(mfIn, 1) == 1)
    mfIn = mfIn(:);
end
if (size(mfOut, 1) == 1)
    mfOut = mfOut(:);
end

% Create an alias for the input vector to be used. As explained above, use
% the MF input with the polyphase interpolator and the MF output otherwise.
if (intpl == 0)
    inVec = mfIn;
else
    inVec = mfOut;
end

%% Optional System Objects for Step-by-step Debugging of the Loop

% Constellation Diagram
if (debug_r)
    hScope = comm.ConstellationDiagram(...
        'SymbolsToDisplaySource', 'Property',...
        'SamplesPerSymbol', 1, ...
        'MeasurementInterval', 256, ...
        'ReferenceConstellation', ...
        Ksym * const);
    hScope.XLimits = [-1.5 1.5]*max(real(const));
    hScope.YLimits = [-1.5 1.5]*max(imag(const));
end

% Time scope used to debug the fractional error
if (debug_r)
    hTScopeCounter = dsp.TimeScope(...
        'Title', 'Fractional Inverval', ...
        'NumInputPorts', 1, ...
        'ShowGrid', 1, ...
        'ShowLegend', 1, ...
        'BufferLength', 1e5, ...
        'TimeSpanOverrunAction', 'Scroll', ...
        'TimeSpan', 1e4, ...
        'TimeUnits', 'None', ...
        'YLimits', [-1 1]);
end

%% Derivative Matched Filter (dMF) - used with the MLTED only
% As explained above, the polyphase interpolator implements the matched
% filtering concurrently with interpolation, so there is no need for a
% dedicated MF block. When using the polyphase interpolator with an MLTED,
% the same applies for the dMF part. Namely, an additional polyphase filter
% is adopted just for the differential matched filtering. The polyphase dMF
% obtains the interpolants jointly with dMF filtering, so there is no need
% to use a dedicated dMF block. In contrast, the other interpolators
% require a dedicated dMF block and process the dMF output directly. For
% such interpolators, implement the dedicated dMF block right here:
if (intpl ~= 0 && strcmp(TED, 'MLTED'))
    mf = rcosdesign(rollOff, rcDelay, L);
    dmf = derivativeMf(mf, L);
    dMfOut = filter(dmf, 1, mfIn);
end

%% Interpolator Design

% Polyphase filter bank
if (intpl == 0)
    % Interpolation factor
    %
    % Note the polyphase filter's interpolation factor is completely
    % independent from the receiver's oversampling factor L. For instance,
    % the receiver may be running with L=2 and the polyphase filter may
    % still apply a significantly higher interpolation factor. Furthermore,
    % note that there is no performance penalty in using a high
    % interpolation factor, aside from using more memory to store the
    % polyphase filter bank. In the end, a single subfilter is used per
    % strobe anyway. A high interpolation factor (e.g., 128) is preferable
    % when the receiver oversampling is low (e.g., L=2).
    polyInterpFactor = 128;

    % Polyphase MF realization
    %
    % Note the RRC filter is not an Lth-band filter, so it is not strictly
    % adequate for a polyphase interpolation filter. However, its cascaded
    % combination with the RRC filter used on the Tx side for pulse shaping
    % yields a raised cosine filter, which is a proper Lth-band filter.
    %
    % The RRC filter is normally designed based on the receiver's
    % oversampling factor L. However, the following RRC is also an
    % interpolator, and the interpolator aims to "divide" each sampling
    % period into polyInterpFactor instants (or phases). Hence, design the
    % filter with a combined oversampling factor of "L * polyInterpFactor".
    % Later on, after applying the polyphase decomposition, each resulting
    % subfilter will end up having an oversampling of only L.
    interpMf = sqrt(polyInterpFactor) * ...
        rcosdesign(rollOff, rcDelay, L * polyInterpFactor);
    polyMf = polyDecomp(interpMf, polyInterpFactor);
    assert(size(polyMf, 1) == polyInterpFactor);

    % Polyphase dMF
    %
    % Each subfilter of the polyphase MF is a phase-offset RRC on its own,
    % equivalent to a phase-offset version of the filter produced by
    % "rcosdesign(rollOff, rcDelay, L)". Correspondingly, to polyphase dMF
    % shall contain the differentiated rows of the polyphase MF.
    polyDMf = zeros(size(polyMf));
    for i = 1:polyInterpFactor
        polyDMf(i, :) = derivativeMf(polyMf(i, :), L);
    end

    % To facilitate the convolution computation using inner products, flip
    % all subfilters (rows of polyMf and polyDMf) from left to right.
    polyMf = fliplr(polyMf);
    polyDMf = fliplr(polyDMf);
else
    polyMf = []; % dummy variable
end

% Quadratic and cubic interpolators
%
% Define the matrix bl_i with the Farrow coefficients to multiply the
% samples surrounding the desired interpolant. Each column of b_mtx holds
% b_l(i) for a fixed l (exponent of mu(k) in 8.76) and for i (neighbor
% sample index) from -2 to 1. After the fliplr operations, the first column
% becomes the one associated with l=0 and the last with l=2. Each of those
% columns are filters to process the samples from x(mk-1) to x(mk+2), i.e.,
% the sample before the basepoint index (x(mk-1)), the sample at the
% basepoint index (x(mk)), and two samples ahead (x(mk+1) and x(mk+2)).
% Before the flipud, the first row of bl_i would have the taps for i=-2,
% which would multiply x(mk+2). For convenience, however, the flipping
% ensures the first row has the taps for i=+1, which multiply x(mk-1). This
% order facilitates the dot product used later.
if (intpl == 2)
    % Farrow coefficients for alpha=0.5 (see Table 8.4.1)
    alpha = 0.5;
    b_mtx = flipud(fliplr(...
        [+alpha,      -alpha, 0; ...
         -alpha, (1 + alpha), 0; ...
         -alpha, (alpha - 1), 1; ...
         +alpha,      -alpha, 0]));
elseif (intpl == 3)
    % Table 8.4.2
    b_mtx = flipud(fliplr(...
        [+1/6,    0, -1/6, 0; ...
         -1/2, +1/2,   +1, 0; ...
         +1/2,   -1, -1/2, 1; ...
         -1/6, +1/2, -1/3, 0]));
else
    b_mtx = []; % dummy variable
end

%% Timing Recovery Loop

% Constants
nSamples = length(inVec);
nSymbols = ceil(nSamples / L);

% Preallocate
xI = zeros(nSymbols, 1); % Output interpolants
mu = zeros(nSymbols, 1); % Fractional symbol timing offset estimate
v  = zeros(nSamples, 1); % PI output
e  = zeros(nSamples, 1); % Error detected by the TED

% Initialize
k      = 0; % interpolant/symbol index
strobe = 0; % strobe signal
cnt    = 1; % modulo-1 counter
vi     = 0; % PI filter integrator

% NOTE: by starting cnt with value 1, the first strobe is asserted when
% "n=L+1" and takes effect on iteration "n=L+2", while setting the
% basepoint index to "m_k=L+1". Furthermore, because the counter step is
% "W=1/L" before the first strobe and "cnt=0" when the first strobe is
% asserted, the first fractional interval estimate from Eq. (8.89) is
% "mu=0". Consequently, the first interpolant tends to be closer (or equal
% to) the sample at the basepoint index m_k, namely to "x(L+1)". This is
% not strictly necessary, but is important to understand, e.g., when
% evaluating the loop in unit tests. Also, this strategy is useful to
% ensure there is enough "memory" when the time comes to compute the first
% interpolant, as the interpolation equations use samples from the past.
%
% Furthermore, note that some TED schemes (ZCTED, GTED, and MMTED) use the
% previous output interpolant or its decision in the error computation.
% Since the first strobe takes effect on iteration "n=L+2", with a
% basepoint at "n=L+1", let the raw input sample at "n=1" be considered as
% the "last output interpolant" when computing the timing error on the
% first strobe, even though "inVec(1)" is not strictly an output
% interpolant. Again, this is an arbitrary choice, which is noteworthy when
% testing the loop. By picking inVec(1) as the starting "last_xI", we can
% compute the timing error right from the first strobe (k=1). Otherwise, we
% would need to compute the error conditionally on "k > 1".
%
% Lastly, note the above "cnt=1" initialization applies only to the first
% iteration. In all other iterations, the counter is always within [0, 1).
last_xI = inVec(1);

% End the loop with enough margin for the computations
%
% Do not process the last L samples when using the ELTED due to the
% look-ahead scheme used to compute the "early" interpolant. When using the
% quadratic or cubic interpolators (which use "x(m_k + 2)"), leave one
% extra sample in the end. In all other cases, process all samples.
if (strcmp(TED, 'ELTED'))
    n_end = nSamples - L;
elseif (intpl > 1)
    n_end = nSamples - 1;
else
    n_end = nSamples;
end

% Start with enough history samples for the interpolator.
%
% As mentioned earlier, the first strobe only takes effect on iteration
% "n=L+2", and the first basepoint index is "m_k=L+1". Hence, the first
% interpolation can access up to L samples from the past. This amount is
% generally sufficient for the linear, quadratic, and cubic interpolators,
% which at maximum access the sample preceding the basepoint index (in this
% case, index "n = m_k - 1 = L"). Thus, the loop can be started right from
% "n=1". Furthermore, this approach works even with the ELTED, ZCTED and
% GTED schemes, which need to compute the zero-crossing (or late)
% interpolants. The zero-crossing interpolant is computed based on the
% basepoint index at "m_k - L/2", so the interpolator only uses up to index
% "m_k - L/2 - 1", which is guaranteed to be available for L >= 2. For
% instance, if L=2, the first basepoint index is "m_k=3", the zero-crossing
% basepoint index is 2, and the interpolator accesses index 1 to compute
% the zero-crossing interpolant.
%
% The only scenario where there may not be enough history samples for the
% interpolation is if using the polyphase interpolator. The sample history
% required by the polyphase interpolator depends on the filter length
% adopted on each polyphase branch. Say, if the polyphase branch filter has
% length N, it processes the samples from index "m_k - N + 1" to m_k. This
% range only works if "m_k - N + 1 >= 1", namely if "m_k >= N". And since
% the first basepoint index occurs after L iterations from the start, the
% loop must start at index "N - L". Furthermore, when using the ELTED,
% ZCTED, or GTED, all of which compute the zero-crossing interpolant using
% basepoint index "m_k - ceil(L/2)", the starting index must be offset by
% another "ceil(L/2)" samples.
if (intpl == 0)
    poly_branch_len = size(polyMf, 2);
    n_start = max(1, poly_branch_len - L);
    if (strcmp(TED, 'ELTED') || strcmp(TED, 'ZCTED') || ...
        strcmp(TED, 'GTED'))
        n_start = n_start + ceil(L/2);
    end
else
    n_start = 1;
end

for n = n_start:n_end
    if strobe == 1
        % Interpolation
        xI(k) = interpolate(intpl, inVec, m_k, mu(k), b_mtx, polyMf);

        % Timing Error Detector:
        a_hat_k = Ksym * slice(xI(k) / Ksym, M); % Data Symbol Estimate
        switch (TED)
            case 'MLTED' % Maximum Likelihood TED
                % dMF interpolant
                if (intpl == 0)
                    xdotI = interpolate(intpl, mfIn, m_k, mu(k), ...
                        b_mtx, polyDMf);
                else
                    xdotI = interpolate(intpl, dMfOut, m_k, mu(k), b_mtx);
                end
                % Decision-directed version of Eq. (8.98), i.e., Eq. (8.27)
                % adapted to complex symbols:
                e(n) = real(a_hat_k) * real(xdotI) + ...
                    imag(a_hat_k) * imag(xdotI);
            case 'ELTED' % Early-late TED
                % Early and late interpolants
                early_idx = m_k + midpointOffset;
                late_idx = m_k - midpointOffset;
                early_mu = mu(k) - muOffset;
                late_mu = mu(k) + muOffset;
                x_early = interpolate(intpl, inVec, early_idx, ...
                    early_mu, b_mtx, polyMf);
                x_late = interpolate(intpl, inVec, late_idx, ...
                    late_mu, b_mtx, polyMf);
                % Decision-directed version of (8.99), i.e., (8.34)
                % adapted to complex symbols:
                e(n) = real(a_hat_k) * (real(x_early) - real(x_late)) + ...
                    imag(a_hat_k) * (imag(x_early) - imag(x_late));
            case 'ZCTED' % Zero-crossing TED
                % Estimate of the previous data symbol
                a_hat_prev = Ksym * slice(last_xI / Ksym, M);

                % Zero-crossing interpolant
                zc_idx = m_k - midpointOffset;
                zc_mu = mu(k) + muOffset;
                x_zc = interpolate(intpl, inVec, zc_idx, zc_mu, ...
                    b_mtx, polyMf);

                % Decision-directed version of (8.100), i.e., (8.37)
                % adapted to complex symbols:
                e(n) = real(x_zc) * ...
                    (real(a_hat_prev) - real(a_hat_k)) + ...
                    imag(x_zc) * (imag(a_hat_prev) - imag(a_hat_k));
            case 'GTED' % Gardner TED
                % Zero-crossing interpolant, same as used by the ZCTED
                zc_idx = m_k - midpointOffset;
                zc_mu = mu(k) + muOffset;
                x_zc = interpolate(intpl, inVec, zc_idx, zc_mu, ...
                    b_mtx, polyMf);

                % Equation (8.101):
                e(n) = real(x_zc) * (real(last_xI) - real(xI(k))) ...
                    + imag(x_zc) * (imag(last_xI) - imag(xI(k)));
            case 'MMTED' % Mueller and Müller TED
                % Estimate of the previous data symbol
                a_hat_prev = Ksym * slice(last_xI / Ksym, M);

                % Decision-directed version of (8.102), i.e., (8.49)
                % adapted to complex symbols:
                e(n) = real(a_hat_prev) * real(xI(k)) - ...
                    real(a_hat_k) * real(last_xI) + ...
                    imag(a_hat_prev) * imag(xI(k)) - ...
                    imag(a_hat_k) * imag(last_xI);
        end

        % Update the "last output interpolant" for the next strobe
        last_xI = xI(k);

        % Real-time debugging scopes
        if (debug_r)
            step(hScope, xI(k))
            step(hTScopeCounter, mu(k));
        end
    else
        % Make the error null on the iterations without a strobe. This is
        % equivalent to upsampling the TED output.
        e(n) = 0;
    end

    % Loop Filter
    vp   = K1 * e(n);        % Proportional
    vi   = vi + (K2 * e(n)); % Integral
    v(n) = vp + vi;          % PI Output
    % NOTE: since e(n)=0 when strobe=0, the PI output can be simplified to
    % "v(n) = vi" on iterations without a strobe. It is only when strobe=1
    % that "vp != 0" and that vi (integrator output) changes. Importantly,
    % note the counter step W below changes briefly to "1/L + vp + vi" when
    % strobe=1 and, then, changes back to "1/L + vi" when strobe=0. In the
    % meantime, "vi" remains constant until the next strobe.

    % Adjust the step used by the modulo-1 counter (see below Eq. 8.86)
    W = 1/L + v(n);

    % Check whether the counter will underflow on the next cycle, i.e.,
    % whenever "cnt < W". When that happens, the strobe signal must
    % indicate the underflow occurrence and trigger updates on:
    %
    % - The basepoint index: set to the index right **before** the
    %   underflow. When strobe=1, it means an underflow will occur on the
    %   **next** cycle. Hence, the index before the underflow is exactly
    %   the current index.
    % - The estimate of the fractional symbol timing offset: the estimate
    %   is based on the counter value **before** the underflow (i.e., on
    %   the current cycle) and the current counter step, according to
    %   equation (8.89).
    strobe = cnt < W;
    if (strobe)
        k = k + 1; % Update the interpolant Index
        m_k = n; % Basepoint index (the index **before** the underflow)
        mu(k) = cnt / W; % Equation (8.89)
    end

    % Next modulo-1 counter value:
    cnt = mod(cnt - W, 1);
end

% Trim the output vector
if (strobe) % ended on a strobe (before filling the k-th interpolant)
    xI = xI(1:k-1);
else
    xI = xI(1:k);
end

%% Static Debug Plots
if (debug_s)
    figure
    plot(e)
    ylabel('Timing Error $e(t)$', 'Interpreter', 'latex')
    xlabel('Symbol $k$', 'Interpreter', 'latex')

    figure
    plot(v)
    title('PI Controller Output')
    ylabel('$v(n)$', 'Interpreter', 'latex')
    xlabel('Sample $n$', 'Interpreter', 'latex')

    figure
    plot(mu, '.')
    title('Fractional Error')
    ylabel('$\mu(k)$', 'Interpreter', 'latex')
    xlabel('Symbol $k$', 'Interpreter', 'latex')
end

end

%% Interpolation
function [xI] = interpolate(method, x, m_k, mu, b_mtx, poly_f)
% [xI] = interpolate(method, x, m_k, mu, b_mtx, poly_h) returns the
% interpolant xI obtained from the vector of samples x.
%
% Args:
%     method -> Interpolation method: polyphase (0), linear (1), quadratic
%               (2), or cubic (3).
%     x      -> Vector of samples based on which the interpolant shall be
%               computed, including the basepoint and surrounding samples.
%     m_k    -> Basepoint index, the index preceding the interpolant.
%     mu     -> Estimated fractional interval between the basepoint index
%               and the desired interpolant instant.
%     b_mtx  -> Matrix with the coefficients for the polynomial
%               interpolator used with method=2 or method=3.
%     poly_f -> Polyphase filter bank that should process the input samples
%               when using the polyphase interpolator (method=0).

    % Adjust the basepoint if mu falls out of the nominal [0,1) range. This
    % step is necessary only to support odd oversampling ratios, when a
    % +-0.5 offset is added to the original mu estimate. In contrast, with
    % an even oversampling ratio, mu is within [0,1) by definition.
    if (mu < 0)
        m_k = m_k - 1;
        mu = mu + 1;
    elseif (mu >= 1)
        m_k = m_k + 1;
        mu = mu - 1;
    end
    assert(mu >= 0 && mu < 1);

    switch (method)
    case 0 % Polyphase interpolator
        % Choose the polyphase subfilter using mu. Use the floor operator
        % to make sure the resulting branch (polyBranch) is always within
        % the acceptable range [1, polyInterpFactor], given that mu is
        % within [0, 1). Also, note it is perfectly feasible to use "round"
        % instead of "floor". In this case, it is only necessary to add an
        % extra subfilter to the filter bank. More specifically, to add a
        % shifted-by-one version of the first subfilter, namely the same
        % subfilter as the first but with a delay shorter by one sampling
        % period (see commit 0537d70). However, this extra complexity is
        % unnecessary, especially if the polyInterpFactor is high enough,
        % when the subfilter phases are already very close to each other.
        polyInterpFactor = size(poly_f, 1);
        polyBranch = floor(polyInterpFactor * mu) + 1;
        polySubfilt = poly_f(polyBranch, :);
        N = length(polySubfilt);
        xI = polySubfilt * x((m_k - N + 1) : m_k);
    case 1 % Linear Interpolator (See Eq. 8.61)
        xI = mu * x(m_k + 1) + (1 - mu) * x(m_k);
    case 2 % Quadratic Interpolator
        % Recursive computation based on Eq. 8.77
        v_l = x(m_k - 1 : m_k + 2).' * b_mtx;
        xI = (v_l(3) * mu + v_l(2)) * mu + v_l(1);
    case 3 % Cubic Interpolator
        % Recursive computation based on Eq. 8.78
        v_l = x(m_k - 1 : m_k + 2).' * b_mtx;
        xI = ((v_l(4) * mu + v_l(3)) * ...
            mu + v_l(2)) * mu + v_l(1);
    end
end

%% Function to map Rx symbols into constellation points
function [z] = slice(y, M)
if (isreal(y))
    % Move the real part of input signal; scale appropriately and round the
    % values to get ideal constellation index
    z_index = round( ((real(y) + (M-1)) ./ 2) );
    % clip the values that are outside the valid range
    z_index(z_index <= -1) = 0;
    z_index(z_index > (M-1)) = M-1;
    % Regenerate Symbol (slice)
    z = z_index*2 - (M-1);
else
    M_bar = sqrt(M);
    % Move the real part of input signal; scale appropriately and round the
    % values to get ideal constellation index
    z_index_re = round( ((real(y) + (M_bar - 1)) ./ 2) );
    % Move the imaginary part of input signal; scale appropriately and
    % round the values to get ideal constellation index
    z_index_im = round( ((imag(y) + (M_bar - 1)) ./ 2) );

    % clip the values that are outside the valid range
    z_index_re(z_index_re <= -1)       = 0;
    z_index_re(z_index_re > (M_bar-1)) = M_bar-1;
    z_index_im(z_index_im <= -1)       = 0;
    z_index_im(z_index_im > (M_bar-1)) = M_bar-1;

    % Regenerate Symbol (slice)
    z = (z_index_re*2 - (M_bar-1)) + 1j*(z_index_im*2 - (M_bar-1));
end
end