|
| 1 | +# SPDX-License-Identifier: Apache-2.0 |
| 2 | +# Copyright (c) 2025-2026 Maurice Garcia |
| 3 | + |
| 4 | +from __future__ import annotations |
| 5 | + |
| 6 | +import math |
| 7 | +from dataclasses import dataclass |
| 8 | + |
| 9 | +import numpy as np |
| 10 | +from numpy.typing import NDArray |
| 11 | +from pydantic import BaseModel, Field |
| 12 | + |
| 13 | +from pypnm.lib.types import ( |
| 14 | + BandwidthHz, |
| 15 | + FloatSeries, |
| 16 | + Microseconds, |
| 17 | + PreEqAtdmaCoefficients, |
| 18 | +) |
| 19 | + |
| 20 | +from pypnm.lib.constants import DOCSIS_ROLL_OFF_FACTOR |
| 21 | + |
| 22 | +MIN_CHANNEL_WIDTH_HZ: BandwidthHz = BandwidthHz(0) |
| 23 | +MIN_TAPS_PER_SYMBOL: int = 0 |
| 24 | +MIN_ROLLOFF: float = 0.0 |
| 25 | +ONE: float = 1.0 |
| 26 | +TWO_PI: float = 2.0 * math.pi |
| 27 | +MICROSECONDS_PER_SECOND: float = 1_000_000.0 |
| 28 | + |
| 29 | + |
| 30 | +class GroupDelayModel(BaseModel): |
| 31 | + """Immutable ATDMA group delay results derived from pre-equalization taps. |
| 32 | +
|
| 33 | + Stores derived timing and delay series used for analysis and reporting. |
| 34 | + """ |
| 35 | + |
| 36 | + channel_width_hz: BandwidthHz = Field(..., description="ATDMA channel width in Hz.") |
| 37 | + rolloff: float = Field(..., description=f"RRC roll-off factor α (typical DOCSIS = {DOCSIS_ROLL_OFF_FACTOR}).") |
| 38 | + taps_per_symbol: int = Field(..., description="Taps per symbol from the pre-EQ header.") |
| 39 | + symbol_rate: float = Field(..., description="Derived symbol rate (sym/s): BW / (1 + rolloff).") |
| 40 | + symbol_time_us: Microseconds = Field(..., description="Derived symbol time in microseconds (µs): 1/symbol_rate.") |
| 41 | + sample_period_us: Microseconds = Field(..., description="Sample period in microseconds (µs): Tsym / taps_per_symbol.") |
| 42 | + fft_size: int = Field(..., description="FFT size used to evaluate the frequency response (N taps).") |
| 43 | + delay_samples: FloatSeries = Field(..., description="Group delay in samples per FFT bin (tap-period units).") |
| 44 | + delay_us: FloatSeries = Field(..., description="Group delay in microseconds per FFT bin.") |
| 45 | + model_config = {"frozen": True} |
| 46 | + |
| 47 | + |
| 48 | +@dataclass(frozen=True, slots=True) |
| 49 | +class GroupDelayCalculator: |
| 50 | + """Compute ATDMA group delay from upstream pre-equalization coefficients. |
| 51 | +
|
| 52 | + This calculator derives **group delay** (the negative slope of the unwrapped |
| 53 | + phase response) from a 24-tap ATDMA upstream FIR equalizer. The input taps are |
| 54 | + complex coefficients (real, imag) taken from `docsIf3CmStatusUsEqData.*` |
| 55 | + after decoding to signed integers (your existing `DocsEqualizerData` class |
| 56 | + already handles endianness + 12/16-bit interpretation and yields taps). |
| 57 | +
|
| 58 | + Conceptually, the equalizer taps represent a discrete-time FIR filter: |
| 59 | +
|
| 60 | + h[n] = re[n] + j·im[n] for n = 0..N-1 |
| 61 | +
|
| 62 | + The processing steps are: |
| 63 | +
|
| 64 | + 1) **Time → Frequency conversion** |
| 65 | + Compute the N-point FFT to obtain the complex frequency response: |
| 66 | +
|
| 67 | + H[k] = FFT{ h[n] } , k = 0..N-1 |
| 68 | +
|
| 69 | + 2) **Phase extraction and unwrap** |
| 70 | + Extract the phase angle of each bin and unwrap it to remove 2π discontinuities: |
| 71 | +
|
| 72 | + φ[k] = unwrap(angle(H[k])) |
| 73 | +
|
| 74 | + 3) **Group delay in samples** |
| 75 | + Group delay is defined as: |
| 76 | +
|
| 77 | + τ(ω) = - dφ(ω) / dω |
| 78 | +
|
| 79 | + With FFT bins, ω[k] = 2π·k/N. We approximate the derivative numerically, |
| 80 | + resulting in group delay measured in **tap-sample periods** (i.e., "samples"). |
| 81 | +
|
| 82 | + 4) **Convert delay from samples → microseconds** |
| 83 | + To express delay in time units, we need the tap sample period. |
| 84 | +
|
| 85 | + For DOCSIS ATDMA upstream channels, symbol rate is typically derived from |
| 86 | + channel width and roll-off (root-raised cosine shaping): |
| 87 | +
|
| 88 | + Rs = BW / (1 + α) |
| 89 | +
|
| 90 | + Then: |
| 91 | +
|
| 92 | + Tsym = 1 / Rs |
| 93 | + Tsamp = Tsym / taps_per_symbol |
| 94 | +
|
| 95 | + Finally: |
| 96 | +
|
| 97 | + delay_us[k] = delay_samples[k] · Tsamp(µs) |
| 98 | +
|
| 99 | + Notes and expectations: |
| 100 | +
|
| 101 | + - This class does **not** assume the main tap location is centered; it reports |
| 102 | + the group delay implied by the taps as provided. |
| 103 | + - The FFT size is set to **N = number of taps** by default. If you later want |
| 104 | + a smoother curve, you can zero-pad (e.g., 128/256 points) without changing |
| 105 | + the underlying physics—only the sampling density in frequency. |
| 106 | + - `taps_per_symbol` comes from the pre-EQ header byte (often 1). |
| 107 | + - `channel_width_hz` must be provided to compute absolute time units (µs). |
| 108 | + Without it, you can still compute delay in samples, but not in seconds. |
| 109 | +
|
| 110 | + Attributes: |
| 111 | + channel_width_hz: ATDMA upstream channel width in Hz (e.g., 1_600_000). |
| 112 | + taps_per_symbol: Tap sampling density per symbol from the pre-EQ header. |
| 113 | + Used to convert symbol time to tap-sample time. |
| 114 | + rolloff: DOCSIS shaping roll-off factor α. Typical default is 0.25. |
| 115 | +
|
| 116 | + Returns: |
| 117 | + A `GroupDelayModel` containing: |
| 118 | + - derived symbol rate/time and sample period |
| 119 | + - group delay arrays per FFT bin in samples and microseconds |
| 120 | + """ |
| 121 | + |
| 122 | + channel_width_hz: BandwidthHz |
| 123 | + taps_per_symbol: int |
| 124 | + rolloff: float = DOCSIS_ROLL_OFF_FACTOR |
| 125 | + |
| 126 | + def __post_init__(self) -> None: |
| 127 | + if int(self.channel_width_hz) <= MIN_CHANNEL_WIDTH_HZ: |
| 128 | + raise ValueError("channel_width_hz must be > 0.") |
| 129 | + if self.taps_per_symbol <= MIN_TAPS_PER_SYMBOL: |
| 130 | + raise ValueError("taps_per_symbol must be > 0.") |
| 131 | + if not math.isfinite(self.rolloff): |
| 132 | + raise ValueError("rolloff must be finite.") |
| 133 | + if self.rolloff < MIN_ROLLOFF: |
| 134 | + raise ValueError("rolloff must be >= 0.") |
| 135 | + |
| 136 | + @staticmethod |
| 137 | + def _to_complex_array(coefficients: list[PreEqAtdmaCoefficients]) -> NDArray[np.complex128]: |
| 138 | + taps: NDArray[np.complex128] = np.empty(len(coefficients), dtype=np.complex128) |
| 139 | + for i, (re, im) in enumerate(coefficients): |
| 140 | + taps[i] = complex(float(re), float(im)) |
| 141 | + return taps |
| 142 | + |
| 143 | + def symbol_rate(self) -> float: |
| 144 | + bw = float(int(self.channel_width_hz)) |
| 145 | + return bw / (ONE + self.rolloff) |
| 146 | + |
| 147 | + def symbol_time_us(self) -> Microseconds: |
| 148 | + sr = self.symbol_rate() |
| 149 | + ts = ONE / sr |
| 150 | + return Microseconds(ts * MICROSECONDS_PER_SECOND) |
| 151 | + |
| 152 | + def sample_period_us(self) -> Microseconds: |
| 153 | + tsym_us = float(self.symbol_time_us()) |
| 154 | + return Microseconds(tsym_us / float(self.taps_per_symbol)) |
| 155 | + |
| 156 | + def compute(self, coefficients: list[PreEqAtdmaCoefficients]) -> GroupDelayModel: |
| 157 | + if len(coefficients) == 0: |
| 158 | + raise ValueError("coefficients cannot be empty.") |
| 159 | + |
| 160 | + h_time = self._to_complex_array(coefficients) |
| 161 | + |
| 162 | + n = int(h_time.shape[0]) |
| 163 | + h_freq = np.fft.fft(h_time, n=n) |
| 164 | + |
| 165 | + phase = np.unwrap(np.angle(h_freq)) |
| 166 | + omega = TWO_PI * (np.arange(n, dtype=np.float64) / float(n)) |
| 167 | + |
| 168 | + dphi_domega = np.gradient(phase, omega) |
| 169 | + delay_samples = -dphi_domega |
| 170 | + |
| 171 | + tsamp_us = float(self.sample_period_us()) |
| 172 | + delay_us = delay_samples * tsamp_us |
| 173 | + |
| 174 | + delay_samples_list: FloatSeries = [float(x) for x in delay_samples.tolist()] |
| 175 | + delay_us_list: FloatSeries = [float(x) for x in delay_us.tolist()] |
| 176 | + |
| 177 | + sr = self.symbol_rate() |
| 178 | + tsym_us = float(self.symbol_time_us()) |
| 179 | + tsamp = float(self.sample_period_us()) |
| 180 | + |
| 181 | + return GroupDelayModel( |
| 182 | + channel_width_hz = BandwidthHz(int(self.channel_width_hz)), |
| 183 | + rolloff = float(self.rolloff), |
| 184 | + taps_per_symbol = int(self.taps_per_symbol), |
| 185 | + symbol_rate = float(sr), |
| 186 | + symbol_time_us = Microseconds(tsym_us), |
| 187 | + sample_period_us = Microseconds(tsamp), |
| 188 | + fft_size = int(n), |
| 189 | + delay_samples = delay_samples_list, |
| 190 | + delay_us = delay_us_list, |
| 191 | + ) |
0 commit comments