Sync

softlora.sync

Synchronization for a Nyquist-rate LoRa receiver.

This module implements the three-stage low-complexity synchronization algorithm of Xhonneux et al., "A Low-Complexity LoRa Synchronization Algorithm Robust to Sampling Time Offsets", IEEE Internet of Things Journal, 2021 (arXiv:1912.11344v2). Algorithm 1 of the paper is transcribed line-by-line; every estimator carries its equation number from the paper.

The three stages are:

1. stage1_detect_and_cfo -- preamble detection on three consecutive upchirps and estimation of the fractional CFO (Eq. 14).
2. stage2_preliminary_sto -- preliminary estimation of the fractional STO from the three upchirps following detection (Eq. 20).
3. stage3_final_sync -- definitive estimation of the integer CFO and STO (Eq. 17, 18) and of the fractional STO (Eq. 20), using the last full preamble upchirp (U7) and the first SFD downchirp (D1), followed by the frequency and time corrections.

The paper's "realign the receiver by lambda_tilde_STO" step of stage 2 is folded into the final time correction: the net timing correction applied to the received signal is L_STO + lambda_sto (the lambda_tilde_STO and -lambda_tilde_STO terms of the paper's final realignment cancel), which yields the same physical payload position while keeping every estimator exact. The paper assumes a standard 8-upchirp preamble (N_preamble_up=8, N_netid=2), in which case D1 is window l+8; the general formula l + (N_preamble_up - N_detect) + N_netid + 1 is used so that preambles with a different number of upchirps (e.g. the 10-upchirp OTA capture) are handled identically.

apply_freq_correction(signal, L_CFO, lambda_cfo, N)

Remove the estimated carrier frequency offset from signal.

The total CFO in bins is L_CFO + lambda_cfo; the signal is shifted in frequency by -(L_CFO + lambda_cfo) * BW / N.

Source code in softlora/sync.py

def apply_freq_correction(signal, L_CFO, lambda_cfo, N):
    """Remove the estimated carrier frequency offset from ``signal``.

    The total CFO in bins is ``L_CFO + lambda_cfo``; the signal is shifted in
    frequency by ``-(L_CFO + lambda_cfo) * BW / N``.
    """
    n = np.arange(len(signal))
    return signal * np.exp(-1j * 2 * np.pi * (L_CFO + lambda_cfo) * n / N)

apply_time_shift(signal, shift_samples)

Delay signal by shift_samples samples (positive = later).

The integer part is applied by dropping/prepending samples and the fractional part by a zero-padded FFT phase ramp. The padding absorbs the circular wrap of the FFT delay so the interior of signal is delayed exactly (a plain FFT phase ramp is circular and would alias the wrap back into the signal, which badly distorts wideband chirps). A negative shift advances the signal.

Source code in softlora/sync.py

def apply_time_shift(signal, shift_samples):
    """Delay ``signal`` by ``shift_samples`` samples (positive = later).

    The integer part is applied by dropping/prepending samples and the
    fractional part by a zero-padded FFT phase ramp.  The padding absorbs the
    circular wrap of the FFT delay so the interior of ``signal`` is delayed
    exactly (a plain FFT phase ramp is circular and would alias the wrap back
    into the signal, which badly distorts wideband chirps).  A negative shift
    advances the signal.
    """
    if shift_samples == 0:
        return signal
    int_shift = int(np.floor(shift_samples))
    frac_shift = shift_samples - int_shift

    if int_shift > 0:
        signal = np.concatenate([np.zeros(int_shift, dtype=signal.dtype), signal])
    elif int_shift < 0:
        signal = signal[-int_shift:]

    if abs(frac_shift) > 1e-10 and len(signal) > 0:
        pad = 128
        x = np.concatenate([
            np.zeros(pad, dtype=complex),
            np.asarray(signal, dtype=complex),
            np.zeros(pad, dtype=complex),
        ])
        X = np.fft.fft(x)
        freqs = np.fft.fftfreq(len(x))
        x2 = np.fft.ifft(X * np.exp(-1j * 2 * np.pi * freqs * frac_shift))
        signal = x2[pad:pad + len(signal)]

    return signal

cfo_bin_products(Y_cur, Y_prev, i, N)

Sum of the five DFT bin products centered on i (Eq. 14).

z = sum_{p=-2}^{2} Y_cur[(i+p) mod N] * conj(Y_prev[(i+p) mod N]). For two consecutive unmodulated upchirps the phase of z equals 2*pi*lambda_CFO, the per-symbol phase accumulation of the fractional CFO.

Source code in softlora/sync.py

def cfo_bin_products(Y_cur, Y_prev, i, N):
    """Sum of the five DFT bin products centered on ``i`` (Eq. 14).

    ``z = sum_{p=-2}^{2} Y_cur[(i+p) mod N] * conj(Y_prev[(i+p) mod N])``.
    For two consecutive unmodulated upchirps the phase of ``z`` equals
    ``2*pi*lambda_CFO``, the per-symbol phase accumulation of the fractional
    CFO.
    """
    z = 0j
    for p in range(-2, 3):
        z += Y_cur[(i + p) % N] * np.conjugate(Y_prev[(i + p) % N])
    return z

est_integer_offsets(s_up, s_down, N)

Integer CFO and STO from the upchirp/downchirp peak bins (Eq. 17, 18).

.. math::

\hat L_{CFO} = \tfrac12 \Gamma_N[(\hat s_{up} + \hat s_{down}) \bmod N]
\hat L_{STO} = (\hat s_{up} - \hat L_{CFO}) \bmod N

with Gamma_N the wrap to [-N/2, N/2) and hat L_STO left in [0, N) as in the paper.

Source code in softlora/sync.py

def est_integer_offsets(s_up, s_down, N):
    """Integer CFO and STO from the upchirp/downchirp peak bins (Eq. 17, 18).

    .. math::

        \\hat L_{CFO} = \\tfrac12 \\Gamma_N[(\\hat s_{up} + \\hat s_{down}) \\bmod N]
        \\hat L_{STO} = (\\hat s_{up} - \\hat L_{CFO}) \\bmod N

    with ``Gamma_N`` the wrap to ``[-N/2, N/2)`` and ``hat L_STO`` left in
    ``[0, N)`` as in the paper.
    """
    sum_mod = (s_up + s_down) % N
    L_CFO = Gamma_N(sum_mod, N) // 2
    L_STO = (s_up - L_CFO) % N
    return L_CFO, L_STO

est_lambda_cfo(z_sum)

Fractional CFO in bins from the summed bin products (Eq. 14).

lambda_hat_CFO = angle(z) / (2*pi), in (-0.5, 0.5].

Source code in softlora/sync.py

def est_lambda_cfo(z_sum):
    """Fractional CFO in bins from the summed bin products (Eq. 14).

    ``lambda_hat_CFO = angle(z) / (2*pi)``, in ``(-0.5, 0.5]``.
    """
    return np.angle(z_sum) / (2 * np.pi)

est_lambda_sto(Y, i, M_hat, N)

Fractional STO in samples from three DFT bins (Eq. 20).

.. math::

\hat\lambda_{STO} = -\Re\left(
    \frac{e^{-j2\pi\hat M/N}Y_{i+1} - e^{j2\pi\hat M/N}Y_{i-1}}
          {2Y_i - e^{-j2\pi\hat M/N}Y_{i+1} - e^{j2\pi\hat M/N}Y_{i-1}}
\right)

M_hat is an estimate of M = N - floor(L_STO + lambda_STO). With M_hat = 0 the estimator reduces to the Jacobsen estimator of Eq. 19.

Source code in softlora/sync.py

def est_lambda_sto(Y, i, M_hat, N):
    """Fractional STO in samples from three DFT bins (Eq. 20).

    .. math::

        \\hat\\lambda_{STO} = -\\Re\\left(
            \\frac{e^{-j2\\pi\\hat M/N}Y_{i+1} - e^{j2\\pi\\hat M/N}Y_{i-1}}
                  {2Y_i - e^{-j2\\pi\\hat M/N}Y_{i+1} - e^{j2\\pi\\hat M/N}Y_{i-1}}
        \\right)

    ``M_hat`` is an estimate of ``M = N - floor(L_STO + lambda_STO)``.  With
    ``M_hat = 0`` the estimator reduces to the Jacobsen estimator of Eq. 19.
    """
    A = np.exp(-1j * 2 * np.pi * M_hat / N) * Y[(i + 1) % N]
    B = np.exp(1j * 2 * np.pi * M_hat / N) * Y[(i - 1) % N]
    denom = 2 * Y[i] - A - B
    if abs(denom) < 1e-12:
        return 0.0
    return -np.real((A - B) / denom)

stage1_detect_and_cfo(iq, downchirp, N, N_detect=3)

Stage 1: preamble detection and fractional CFO estimation.

Slides N-sample windows over iq, dechirps each with the reference downchirp and computes z^l (Eq. 14) for every window that has a predecessor. A preamble is declared when N_detect consecutive windows peak within one FFT bin. The fractional CFO is then the phase of the sum of the last two bin-product vectors (Eq. 14):

lambda_hat_CFO = angle(z^l + z^{l-1}) / (2*pi).

Parameters:

Name	Type	Description	Default
iq	ndarray	Complex baseband signal at fs = BW (N samples per symbol).	required
downchirp	ndarray	N-sample reference downchirp (conjugate of the upchirp).	required
N	int	Samples per symbol (2**sf).	required
N_detect	int	Number of consecutive matching windows required.	3

Returns:

Name	Type	Description
l	int or None	Index of the last window of the detected preamble run (None if no preamble was found).
symbols	list of int	Argmax bin of every processed window.
fft_results	list of ndarray	Dechirped window FFTs of every processed window.
lambda_cfo	float or None	Fractional CFO in bins (units of BW/N), or None if undetected.

Source code in softlora/sync.py

def stage1_detect_and_cfo(iq, downchirp, N, N_detect=3):
    """Stage 1: preamble detection and fractional CFO estimation.

    Slides N-sample windows over ``iq``, dechirps each with the reference
    downchirp and computes ``z^l`` (Eq. 14) for every window that has a
    predecessor.  A preamble is declared when ``N_detect`` consecutive windows
    peak within one FFT bin.  The fractional CFO is then the phase of the sum
    of the last two bin-product vectors (Eq. 14):

    ``lambda_hat_CFO = angle(z^l + z^{l-1}) / (2*pi)``.

    Parameters
    ----------
    iq : ndarray
        Complex baseband signal at fs = BW (N samples per symbol).
    downchirp : ndarray
        N-sample reference downchirp (conjugate of the upchirp).
    N : int
        Samples per symbol (2**sf).
    N_detect : int
        Number of consecutive matching windows required.

    Returns
    -------
    l : int or None
        Index of the last window of the detected preamble run (None if no
        preamble was found).
    symbols : list of int
        Argmax bin of every processed window.
    fft_results : list of ndarray
        Dechirped window FFTs of every processed window.
    lambda_cfo : float or None
        Fractional CFO in bins (units of BW/N), or None if undetected.
    """
    num_windows = len(iq) // N
    symbols = []
    fft_results = []
    z_list = []

    l = None
    for w in range(num_windows):
        window = iq[w * N:(w + 1) * N]
        if len(window) < N:
            break
        Y = np.fft.fft(window * downchirp)
        s = int(np.argmax(np.abs(Y)))
        symbols.append(s)
        fft_results.append(Y)
        if w >= 1:
            z_list.append(cfo_bin_products(Y, fft_results[w - 1], s, N))

        if w >= N_detect - 1:
            vals = symbols[w - N_detect + 1:w + 1]
            if all(abs(v - vals[0]) <= 1 for v in vals):
                l = w
                break

    if l is None:
        return None, symbols, fft_results, None

    lambda_cfo = est_lambda_cfo(z_list[l - 1] + z_list[l - 2])
    return l, symbols, fft_results, lambda_cfo

stage2_preliminary_sto(iq, downchirp, N, l, lambda_cfo, N_detect=3)

Stage 2: first correction of the fractional STO.

Dechirps the N_detect windows following the detected preamble with the fractional CFO removed (e^{-j2*pi*lambda_hat_CFO * n/N} with n the absolute sample index) and averages their spectra. s_tilde_up is the peak bin of the first of these windows (Algorithm 1 step 12), giving M_tilde = N - s_tilde_up; the preliminary fractional STO is then estimated on the averaged spectrum with Eq. 20.

Parameters:

Name	Type	Description	Default
iq	ndarray	Complex baseband signal at fs = BW.	required
downchirp	ndarray	N-sample reference downchirp.	required
N	int	Samples per symbol.	required
l	int	Last window index of the detected preamble (from stage 1).	required
lambda_cfo	float	Fractional CFO in bins (from stage 1).	required
N_detect	int	Number of windows to average.	3

Returns:

Name	Type	Description
s_tilde_up	int	Peak bin of the first post-detection window (used as the stage-2 approximation of L_STO).
lambda_sto_prelim	float	Preliminary fractional STO in samples (Eq. 20 with M_tilde).
Y_avg	ndarray	CFO-corrected, averaged spectrum of the post-detection windows (used by stage 3 to compute the definitive fractional STO).

Source code in softlora/sync.py

def stage2_preliminary_sto(iq, downchirp, N, l, lambda_cfo, N_detect=3):
    """Stage 2: first correction of the fractional STO.

    Dechirps the ``N_detect`` windows following the detected preamble with the
    fractional CFO removed (``e^{-j2*pi*lambda_hat_CFO * n/N}`` with ``n`` the
    absolute sample index) and averages their spectra.  ``s_tilde_up`` is the
    peak bin of the *first* of these windows (Algorithm 1 step 12), giving
    ``M_tilde = N - s_tilde_up``; the preliminary fractional STO is then
    estimated on the averaged spectrum with Eq. 20.

    Parameters
    ----------
    iq : ndarray
        Complex baseband signal at fs = BW.
    downchirp : ndarray
        N-sample reference downchirp.
    N : int
        Samples per symbol.
    l : int
        Last window index of the detected preamble (from stage 1).
    lambda_cfo : float
        Fractional CFO in bins (from stage 1).
    N_detect : int
        Number of windows to average.

    Returns
    -------
    s_tilde_up : int
        Peak bin of the first post-detection window (used as the stage-2
        approximation of ``L_STO``).
    lambda_sto_prelim : float
        Preliminary fractional STO in samples (Eq. 20 with ``M_tilde``).
    Y_avg : ndarray
        CFO-corrected, averaged spectrum of the post-detection windows (used
        by stage 3 to compute the definitive fractional STO).
    """
    n = np.arange(N)
    Y_avg = np.zeros(N, dtype=complex)
    for m in range(1, N_detect + 1):
        idx = l + m
        offset = idx * N
        window = iq[offset:offset + N]
        if len(window) < N:
            window = np.pad(window, (0, N - len(window)))
        phase_corr = np.exp(-1j * 2 * np.pi * lambda_cfo * (offset + n) / N)
        Y_avg += np.fft.fft(window * phase_corr * downchirp)

    offset1 = (l + 1) * N
    window1 = iq[offset1:offset1 + N]
    if len(window1) < N:
        window1 = np.pad(window1, (0, N - len(window1)))
    phase1 = np.exp(-1j * 2 * np.pi * lambda_cfo * (offset1 + n) / N)
    Y1 = np.fft.fft(window1 * phase1 * downchirp)

    s_tilde_up = int(np.argmax(np.abs(Y1)))
    M_tilde = N - s_tilde_up
    lambda_sto_prelim = est_lambda_sto(Y_avg, s_tilde_up, M_tilde, N)
    return s_tilde_up, lambda_sto_prelim, Y_avg

stage3_final_sync(iq, downchirp, N, l, lambda_cfo, Y_avg, N_preamble_up=8, N_netid=2, N_detect=3, down_shift=0)

Stage 3: integer CFO/STO and definitive fractional STO.

Demodulates the final full preamble upchirp (window l + up_shift: U7 = l+4 when the preamble has at least 7 upchirps, the last upchirp l + (N_preamble_up - N_detect) otherwise) and the first SFD downchirp D1 (window l + (N_preamble_up - N_detect) + N_netid + 1, equal to l+8 for the standard 8-upchirp preamble) with the fractional CFO removed, resolves the integer offsets with Eq. 17/18, and computes the definitive fractional STO with Eq. 20 using the averaged spectrum Y_avg stored by stage 2 and M_hat = N - L_STO.

down_shift selects which SFD downchirp window is demodulated (0 = D1, the paper's window; 1 = D2, used by the decoder to disambiguate an integer STO falling exactly on the N/2 wrap boundary).

Parameters:

Name	Type	Description	Default
iq	ndarray	Complex baseband signal at fs = BW.	required
downchirp	ndarray	N-sample reference downchirp.	required
N	int	Samples per symbol.	required
l	int	Last window index of the detected preamble (from stage 1).	required
lambda_cfo	float	Fractional CFO in bins (from stage 1).	required
Y_avg	ndarray	Averaged spectrum from stage 2.	required
N_preamble_up	int	Number of preamble upchirps.	8
N_netid	int	Number of network-identifier symbols before the SFD.	2
N_detect	int	Preamble detection window count.	3
down_shift	int	Offset (in symbols) from D1 of the downchirp window to demodulate.	0

Returns:

Name	Type	Description
L_CFO	int	Integer CFO in bins (Eq. 17).
L_STO	int	Integer STO in samples, in [0, N) (Eq. 18).
lambda_sto	float	Definitive fractional STO in samples (Eq. 20).
s_up	int	Peak bin of the final preamble upchirp (U7, or the last upchirp for preambles with fewer than 7 upchirps).
s_down	int	Peak bin of the SFD downchirp (D1 or D2).
M_hat	int	N - L_STO, the STO boundary estimate used in Eq. 20.

Source code in softlora/sync.py

def stage3_final_sync(iq, downchirp, N, l, lambda_cfo, Y_avg,
                      N_preamble_up=8, N_netid=2, N_detect=3,
                      down_shift=0):
    """Stage 3: integer CFO/STO and definitive fractional STO.

    Demodulates the final full preamble upchirp (window ``l + up_shift``:
    U7 = ``l+4`` when the preamble has at least 7 upchirps, the last upchirp
    ``l + (N_preamble_up - N_detect)`` otherwise) and the first SFD downchirp
    D1 (window ``l + (N_preamble_up - N_detect) + N_netid + 1``, equal to
    ``l+8`` for the standard 8-upchirp preamble) with the fractional CFO
    removed, resolves the integer offsets with Eq. 17/18, and computes the
    definitive fractional STO with Eq. 20 using the averaged spectrum
    ``Y_avg`` stored by stage 2 and ``M_hat = N - L_STO``.

    ``down_shift`` selects which SFD downchirp window is demodulated (0 = D1,
    the paper's window; 1 = D2, used by the decoder to disambiguate an integer
    STO falling exactly on the ``N/2`` wrap boundary).

    Parameters
    ----------
    iq : ndarray
        Complex baseband signal at fs = BW.
    downchirp : ndarray
        N-sample reference downchirp.
    N : int
        Samples per symbol.
    l : int
        Last window index of the detected preamble (from stage 1).
    lambda_cfo : float
        Fractional CFO in bins (from stage 1).
    Y_avg : ndarray
        Averaged spectrum from stage 2.
    N_preamble_up : int
        Number of preamble upchirps.
    N_netid : int
        Number of network-identifier symbols before the SFD.
    N_detect : int
        Preamble detection window count.
    down_shift : int
        Offset (in symbols) from D1 of the downchirp window to demodulate.

    Returns
    -------
    L_CFO : int
        Integer CFO in bins (Eq. 17).
    L_STO : int
        Integer STO in samples, in ``[0, N)`` (Eq. 18).
    lambda_sto : float
        Definitive fractional STO in samples (Eq. 20).
    s_up : int
        Peak bin of the final preamble upchirp (U7, or the last upchirp for
        preambles with fewer than 7 upchirps).
    s_down : int
        Peak bin of the SFD downchirp (D1 or D2).
    M_hat : int
        ``N - L_STO``, the STO boundary estimate used in Eq. 20.
    """
    n = np.arange(N)

    # U7 (window l+4) per Algorithm 1 when the preamble has at least 7
    # upchirps; otherwise the last full upchirp (window
    # l + N_preamble_up - N_detect) so preambles with fewer upchirps
    # (e.g. 6 upchirps + 2 net-id + 2.25 SFD) are handled identically.
    up_shift = min(4, N_preamble_up - N_detect)
    up_offset = (l + up_shift) * N
    window_up = iq[up_offset:up_offset + N]
    if len(window_up) < N:
        window_up = np.pad(window_up, (0, N - len(window_up)))
    phase_up = np.exp(-1j * 2 * np.pi * lambda_cfo * (up_offset + n) / N)
    Y_up = np.fft.fft(window_up * phase_up * downchirp)
    s_up = int(np.argmax(np.abs(Y_up)))

    d1 = (N_preamble_up - N_detect) + N_netid + 1 + down_shift
    down_offset = (l + d1) * N
    window_dn = iq[down_offset:down_offset + N]
    if len(window_dn) < N:
        window_dn = np.pad(window_dn, (0, N - len(window_dn)))
    phase_dn = np.exp(-1j * 2 * np.pi * lambda_cfo * (down_offset + n) / N)
    Y_dn = np.fft.fft(window_dn * phase_dn * np.conj(downchirp))
    s_down = int(np.argmax(np.abs(Y_dn)))

    L_CFO, L_STO = est_integer_offsets(s_up, s_down, N)
    M_hat = N - L_STO
    lambda_sto = est_lambda_sto(Y_avg, s_up, M_hat, N)

    return L_CFO, L_STO, lambda_sto, s_up, s_down, M_hat