Wander TDEV and MTIE Analyzers

src/vse_sync_pp/analyzers/__init__.py

@@ -19,5 +19,13 @@
         ts2phc.TimeErrorAnalyzer,
         phc2sys.TimeErrorAnalyzer,
         pmc.ClockStateAnalyzer,
+        gnss.TimeDeviationAnalyzer,
+        ppsdpll.TimeDeviationAnalyzer,
+        ts2phc.TimeDeviationAnalyzer,
+        phc2sys.TimeDeviationAnalyzer,
+        gnss.MaxTimeIntervalErrorAnalyzer,
+        ppsdpll.MaxTimeIntervalErrorAnalyzer,
+        ts2phc.MaxTimeIntervalErrorAnalyzer,
+        phc2sys.MaxTimeIntervalErrorAnalyzer,
     )
 }

src/vse_sync_pp/analyzers/analyzer.py

@@ -6,6 +6,11 @@
 from pandas import DataFrame
 from datetime import (datetime, timezone)
 
+import allantools
+import numpy as np
+
+from scipy import signal as scipy_signal
+
 from ..requirements import REQUIREMENTS
 
 
@@ -265,3 +270,208 @@ def explain(self, data):
             'duration': data.iloc[-1].timestamp - data.iloc[0].timestamp,
             'terror': self._statistics(data.terror, 'ns'),
         }
+
+
+def calculate_limit(accuracy, limit_percentage, tau):
+    """Calculate upper limit based on tau
+
+    `accuracy` is the list of functions to calculate upper limits
+    `limit_percentage` is the unaccuracy percentage
+    `tau` is the observation window interval
+
+    Return the upper limit value based on `tau`
+    """
+    for (low, high), f in accuracy.items():
+        if ((low is None or tau > low) and (tau <= high)):
+            return f(tau) * (limit_percentage / 100)
+
+
+def out_of_range(taus, samples, accuracy, limit):
+    """Check if the input samples are out of range.
+
+    `taus` list of observation windows intervals
+    `samples` are input samples
+    `accuracy` contains the list of upper bound limit functions
+    `limit` is the percentage to apply the upper limit
+
+    Return `True` if any value in `samples` is out of range
+    """
+    for tau, sample in zip(taus, samples):
+        mask = calculate_limit(accuracy, limit, tau)
+        if mask <= sample:
+            return True
+    return False
+
+
+def calculate_filter(input_signal, transient, sample_rate):
+    """Calculate digital low-pass filter from `input_signal`
+
+    scipy_signal.butter input arguments:
+        order of the filter: 1
+        critical Frequency in half cycles per sample:
+        (for digital filters is normalized from 0 to 1 where 1 is Nyquist frequency)
+          cutoff Frequency: 0.1Hz
+          sample_rate in samples per second
+        btype is band type or type of filter
+        analog=False since it is always a digital filter
+    scipy_signal.butter return arguments:
+        `numerator` coefficient vector and `denominator coefficient vector of the butterworth digital filter
+        """
+    numerator, denominator = scipy_signal.butter(1, 0.1 / (sample_rate / 2), btype="low", analog=False, output="ba")
+    lpf_signal = scipy_signal.filtfilt(numerator, denominator, input_signal.terror)
+    lpf_signal = lpf_signal[transient:len(lpf_signal)]
+    return lpf_signal
+
+
+class TimeIntervalErrorAnalyzerBase(Analyzer):
+    """Analyze Time Interval Error (also referred to as Wander).
+
+    Derived classes calculate specific Time Interval Error metric focused on measuring
+    the change of Time Error.
+    """
+    locked = frozenset()
+
+    def __init__(self, config):
+        super().__init__(config)
+        # samples in the initial transient period are ignored
+        self._transient = config.parameter('transient-period/s')
+        # minimum test duration for a valid test
+        self._duration_min = config.parameter('min-test-duration/s')
+        # limit initial tau observation windows from 1 to 10k taus
+        taus_below_10k = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90,
+                                  100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000,
+                                  4000, 5000, 6000, 7000, 8000, 9000, 10000])
+        # observation window upper limit to 100k samples in 5k increments
+        taus_above_10k = np.arange(15000, 100000, 5000)
+        # `taus_list` contains range limit of obervation window intervals for which to compute the statistic
+        self._taus_list = np.concatenate((taus_below_10k, taus_above_10k))
+        self._rate = None
+        self._lpf_signal = None
+
+    def prepare(self, rows):
+        idx = 0
+        try:
+            tstart = rows[0].timestamp + self._transient
+        except IndexError:
+            pass
+        else:
+            while idx < len(rows):
+                if tstart <= rows[idx].timestamp:
+                    break
+                idx += 1
+        return super().prepare(rows[idx:])
+
+    @staticmethod
+    def calculate_rate(data):
+        # calculate sample rate using 100 samples
+        cumdelta = 0
+
+        for i in range(1, 100 if len(data) > 100 else len(data)):
+            cumdelta = cumdelta + data.iloc[i].timestamp - data.iloc[i - 1].timestamp
+        return round((1 / (cumdelta / 100)))
+
+    def _test_common(self, data):
+        if len(data) == 0:
+            return ("error", "no data")
+        if frozenset(data.state.unique()).difference(self.locked):
+            return (False, "loss of lock")
+        if data.iloc[-1].timestamp - data.iloc[0].timestamp < self._duration_min:
+            return (False, "short test duration")
+        if len(data) - 1 < self._duration_min:
+            return (False, "short test samples")
+        if self._rate is None:
+            self._rate = self.calculate_rate(data)
+        if self._lpf_signal is None:
+            self._lpf_signal = calculate_filter(data, self._transient, self._rate)
+        return None
+
+    def _explain_common(self, data):
+        if len(data) == 0:
+            return {}
+        if self._rate is None:
+            self._rate = self.calculate_rate(data)
+        if self._lpf_signal is None:
+            self._lpf_signal = calculate_filter(data, self._transient, self._rate)
+        return None
+
+
+class TimeDeviationAnalyzerBase(TimeIntervalErrorAnalyzerBase):
+    """Analyze Time Deviation (TDEV).
+
+    Derived classes must override class attribute `locked`, specifying a
+    frozenset of values representing locked states.
+    """
+    def __init__(self, config):
+        super().__init__(config)
+        # required system time deviation output
+        self._accuracy = config.requirement('time-deviation-in-locked-mode/ns')
+        # limit of inaccuracy at observation point
+        self._limit = config.parameter('time-deviation-limit/%')
+        # list of observation windows intervals to calculate TDEV
+        # `_taus` is a subset of `taus_list`
+        self._taus = None
+        # TDEV samples
+        self._samples = None
+
+    def test(self, data):
+        result = self._test_common(data)
+        if result is None:
+            if self._samples is None:
+                self._taus, self._samples, errors, ns = allantools.tdev(self._lpf_signal, rate=self._rate, data_type="phase", taus=self._taus_list) # noqa
+            if out_of_range(self._taus, self._samples, self._accuracy, self._limit):
+                return (False, "unacceptable time deviation")
+            return (True, None)
+        return result
+
+    def explain(self, data):
+        analysis = self._explain_common(data)
+        if analysis is None:
+            if self._samples is None:
+                self._taus, self._samples, errors, ns = allantools.tdev(self._lpf_signal, rate=self._rate, data_type="phase", taus=self._taus_list) # noqa
+            return {
+                'timestamp': self._timestamp_from_dec(data.iloc[0].timestamp),
+                'duration': data.iloc[-1].timestamp - data.iloc[0].timestamp,
+                'tdev': self._statistics(self._samples, 'ns'),
+            }
+        return analysis
+
+
+class MaxTimeIntervalErrorAnalyzerBase(TimeIntervalErrorAnalyzerBase):
+    """Analyze Maximum Time Interval Error (MTIE).
+
+    Derived classes must override class attribute `locked`, specifying a
+    frozenset of values representing locked states.
+    """
+    def __init__(self, config):
+        super().__init__(config)
+        # required system maximum time interval error output in us
+        self._accuracy = config.requirement('maximum-time-interval-error-in-locked-mode/us')
+        # limit of inaccuracy at observation point
+        self._limit = config.parameter('maximum-time-interval-error-limit/%')
+        # list of observation windows intervals to calculate MTIE
+        # `_taus` will be a subset of `taus_list`
+        self._taus = None
+        # MTIE samples
+        self._samples = None
+
+    def test(self, data):
+        result = self._test_common(data)
+        if result is None:
+            if self._samples is None:
+                self._taus, self._samples, errors, ns = allantools.mtie(self._lpf_signal, rate=self._rate, data_type="phase", taus=self._taus_list) # noqa
+            if out_of_range(self._taus, self._samples, self._accuracy, self._limit):
+                return (False, "unacceptable mtie")
+            return (True, None)
+        return result
+
+    def explain(self, data):
+        analysis = self._explain_common(data)
+        if analysis is None:
+            if self._samples is None:
+                self._taus, self._samples, errors, ns = allantools.mtie(self._lpf_signal, rate=self._rate, data_type="phase", taus=self._taus_list) # noqa
+            return {
+                'timestamp': self._timestamp_from_dec(data.iloc[0].timestamp),
+                'duration': data.iloc[-1].timestamp - data.iloc[0].timestamp,
+                'mtie': self._statistics(self._samples, 'ns'),
+            }
+        return analysis

src/vse_sync_pp/analyzers/gnss.py

@@ -3,6 +3,8 @@
 """Analyze GNSS log messages"""
 
 from .analyzer import TimeErrorAnalyzerBase
+from .analyzer import TimeDeviationAnalyzerBase
+from .analyzer import MaxTimeIntervalErrorAnalyzerBase
 
 
 class TimeErrorAnalyzer(TimeErrorAnalyzerBase):
@@ -17,3 +19,19 @@ class TimeErrorAnalyzer(TimeErrorAnalyzerBase):
     # 4 = GPS + dead reckoning combined
     # 5 = time only fix
     locked = frozenset({3, 4, 5})
+
+
+class TimeDeviationAnalyzer(TimeDeviationAnalyzerBase):
+    """Analyze time deviation"""
+    id_ = 'gnss/time-deviation'
+    parser = 'gnss/time-error'
+    # see 'state' values in `TimeErrorAnalyzer` comments
+    locked = frozenset({3, 4, 5})
+
+
+class MaxTimeIntervalErrorAnalyzer(MaxTimeIntervalErrorAnalyzerBase):
+    """Analyze time deviation"""
+    id_ = 'gnss/mtie'
+    parser = 'gnss/time-error'
+    # see 'state' values in `TimeErrorAnalyzer` comments
+    locked = frozenset({3, 4, 5})

src/vse_sync_pp/analyzers/phc2sys.py

@@ -3,10 +3,26 @@
 """Analyze phc2sys log messages"""
 
 from .analyzer import TimeErrorAnalyzerBase
+from .analyzer import TimeDeviationAnalyzerBase
+from .analyzer import MaxTimeIntervalErrorAnalyzerBase
 
 
 class TimeErrorAnalyzer(TimeErrorAnalyzerBase):
     """Analyze time error"""
     id_ = 'phc2sys/time-error'
     parser = id_
     locked = frozenset({'s2'})
+
+
+class TimeDeviationAnalyzer(TimeDeviationAnalyzerBase):
+    """Analyze time deviation"""
+    id_ = 'phc2sys/time-deviation'
+    parser = 'phc2sys/time-error'
+    locked = frozenset({'s2'})
+
+
+class MaxTimeIntervalErrorAnalyzer(MaxTimeIntervalErrorAnalyzerBase):
+    """Analyze max time interval error"""
+    id_ = 'phc2sys/mtie'
+    parser = 'phc2sys/time-error'
+    locked = frozenset({'s2'})

src/vse_sync_pp/analyzers/ppsdpll.py

@@ -3,6 +3,8 @@
 """Analyze ppsdpll log messages"""
 
 from .analyzer import TimeErrorAnalyzerBase
+from .analyzer import TimeDeviationAnalyzerBase
+from .analyzer import MaxTimeIntervalErrorAnalyzerBase
 
 
 class TimeErrorAnalyzer(TimeErrorAnalyzerBase):
@@ -24,3 +26,29 @@ def prepare(self, rows):
         return super().prepare([
             r._replace(terror=float(r.terror)) for r in rows
         ])
+
+
+class TimeDeviationAnalyzer(TimeDeviationAnalyzerBase):
+    """Analyze DPLL time deviation"""
+    id_ = 'ppsdpll/time-deviation'
+    parser = 'dpll/time-error'
+    # see 'state' values in `TimeErrorAnalyzer` comments
+    locked = frozenset({2, 3})
+
+    def prepare(self, rows):
+        return super().prepare([
+            r._replace(terror=float(r.terror)) for r in rows
+        ])
+
+
+class MaxTimeIntervalErrorAnalyzer(MaxTimeIntervalErrorAnalyzerBase):
+    """Analyze DPLL max time interval error"""
+    id_ = 'ppsdpll/mtie'
+    parser = 'dpll/time-error'
+    # see 'state' values in `TimeErrorAnalyzer` comments
+    locked = frozenset({2, 3})
+
+    def prepare(self, rows):
+        return super().prepare([
+            r._replace(terror=float(r.terror)) for r in rows
+        ])

src/vse_sync_pp/analyzers/ts2phc.py

@@ -1,12 +1,27 @@
 ### SPDX-License-Identifier: GPL-2.0-or-later
 
 """Analyze ts2phc log messages"""
-
 from .analyzer import TimeErrorAnalyzerBase
+from .analyzer import TimeDeviationAnalyzerBase
+from .analyzer import MaxTimeIntervalErrorAnalyzerBase
 
 
 class TimeErrorAnalyzer(TimeErrorAnalyzerBase):
     """Analyze time error"""
     id_ = 'ts2phc/time-error'
     parser = id_
     locked = frozenset({'s2'})
+
+
+class TimeDeviationAnalyzer(TimeDeviationAnalyzerBase):
+    """Analyze time deviation"""
+    id_ = 'ts2phc/time-deviation'
+    parser = 'ts2phc/time-error'
+    locked = frozenset({'s2'})
+
+
+class MaxTimeIntervalErrorAnalyzer(MaxTimeIntervalErrorAnalyzerBase):
+    """Analyze max time interval error"""
+    id_ = 'ts2phc/mtie'
+    parser = 'ts2phc/time-error'
+    locked = frozenset({'s2'})

src/vse_sync_pp/requirements.py

@@ -1,15 +1,39 @@
 ### SPDX-License-Identifier: GPL-2.0-or-later
 
-"""Requirements specified in standards"""
+"""Requirements specified in ITU-T G.8272/Y.1367"""
 
 REQUIREMENTS = {
     'G.8272/PRTC-A': {
+        'maximum-time-interval-error-in-locked-mode/us': {
+            (None, 273): lambda t: 0.000275 * t + 0.025,
+            (274, None): lambda t: 0.10
+        },
+        'time-deviation-in-locked-mode/ns': {
+            (None, 100): lambda t: 3,
+            (101, 1000): lambda t: 0.03 * t,
+            (1001, 10000): lambda t: 30
+        },
         'time-error-in-locked-mode/ns': 100,
     },
     'G.8272/PRTC-B': {
+        'maximum-time-interval-error-in-locked-mode/us': {
+            (None, 54.5): lambda t: 0.000275 * t + 0.025,
+            (54.5, None): lambda t: 0.04
+        },
+        'time-deviation-in-locked-mode/ns': {
+            (None, 100): lambda t: 1,
+            (101, 500): lambda t: 0.01 * t,
+            (501, 100000): lambda t: 5
+        },
         'time-error-in-locked-mode/ns': 40,
     },
     'workload/RAN': {
         'time-error-in-locked-mode/ns': 100,
+        'time-deviation-in-locked-mode/ns': {
+            (None, 100000): lambda t: 100
+        },
+        'maximum-time-interval-error-in-locked-mode/us': {
+            (None, 100000): lambda t: 1
+        }
     },
 }