Coverage for src / tracekit / analyzers / jitter / spectrum.py: 100%
60 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-11 23:04 +0000
1"""Jitter spectrum analysis.
3This module provides FFT-based analysis of TIE data to identify
4periodic jitter sources and their frequencies.
7Example:
8 >>> from tracekit.analyzers.jitter.spectrum import jitter_spectrum
9 >>> result = jitter_spectrum(tie_data, sample_rate=1e9)
10 >>> print(f"Dominant frequency: {result.dominant_frequency / 1e3:.1f} kHz")
12References:
13 IEEE 2414-2020: Standard for Jitter and Phase Noise
14"""
16from __future__ import annotations
18from dataclasses import dataclass
19from typing import TYPE_CHECKING
21import numpy as np
22from scipy.signal import find_peaks
24if TYPE_CHECKING:
25 from numpy.typing import NDArray
28@dataclass
29class JitterSpectrumResult:
30 """Result of jitter spectrum analysis.
32 Attributes:
33 frequencies: Frequency array in Hz.
34 magnitude: Magnitude spectrum in seconds.
35 magnitude_db: Magnitude spectrum in dB (relative to 1 ps).
36 dominant_frequency: Frequency of largest component in Hz.
37 dominant_magnitude: Magnitude at dominant frequency in seconds.
38 noise_floor: Estimated noise floor in seconds.
39 peaks: List of (frequency, magnitude) tuples for detected peaks.
40 """
42 frequencies: NDArray[np.float64]
43 magnitude: NDArray[np.float64]
44 magnitude_db: NDArray[np.float64]
45 dominant_frequency: float | None
46 dominant_magnitude: float | None
47 noise_floor: float
48 peaks: list[tuple[float, float]]
51def jitter_spectrum(
52 tie_data: NDArray[np.float64],
53 sample_rate: float,
54 *,
55 window: str = "hann",
56 detrend: bool = True,
57 n_peaks: int = 10,
58) -> JitterSpectrumResult:
59 """Compute FFT of TIE data to identify jitter frequency components.
61 Identifies periodic jitter sources by frequency, useful for
62 debugging EMI-induced jitter and power supply noise.
64 Args:
65 tie_data: Time Interval Error data in seconds.
66 sample_rate: Sample rate of TIE data (edges per second).
67 window: Window function ("hann", "hamming", "blackman", "none").
68 detrend: Remove linear trend from data before FFT.
69 n_peaks: Number of peaks to identify.
71 Returns:
72 JitterSpectrumResult with frequency analysis.
74 Example:
75 >>> result = jitter_spectrum(tie_data, sample_rate=1e9)
76 >>> for freq, mag in result.peaks:
77 ... print(f"{freq/1e3:.1f} kHz: {mag*1e12:.2f} ps")
79 References:
80 IEEE 2414-2020 Section 6.8
81 """
82 valid_data = tie_data[~np.isnan(tie_data)]
83 n = len(valid_data)
85 if n < 16:
86 return JitterSpectrumResult(
87 frequencies=np.array([]),
88 magnitude=np.array([]),
89 magnitude_db=np.array([]),
90 dominant_frequency=None,
91 dominant_magnitude=None,
92 noise_floor=0.0,
93 peaks=[],
94 )
96 # Detrend if requested
97 if detrend:
98 # Remove linear trend
99 x = np.arange(n)
100 slope, intercept = np.polyfit(x, valid_data, 1)
101 data_detrended = valid_data - (slope * x + intercept)
102 else:
103 data_detrended = valid_data - np.mean(valid_data)
105 # Apply window
106 if window == "hann":
107 win = np.hanning(n)
108 elif window == "hamming":
109 win = np.hamming(n)
110 elif window == "blackman":
111 win = np.blackman(n)
112 else:
113 win = np.ones(n)
115 # Compensate for window power loss
116 window_factor = np.sqrt(np.mean(win**2))
117 data_windowed = data_detrended * win
119 # Zero-pad to next power of 2
120 nfft = int(2 ** np.ceil(np.log2(n)))
122 # Compute FFT
123 spectrum = np.fft.rfft(data_windowed, n=nfft)
124 frequencies = np.fft.rfftfreq(nfft, d=1.0 / sample_rate)
126 # Calculate magnitude spectrum
127 # Scale for proper amplitude: 2/N for single-sided, compensate for window
128 magnitude = np.abs(spectrum) * 2 / n / window_factor
130 # Convert to dB (relative to 1 ps = 1e-12 s)
131 reference = 1e-12 # 1 ps
132 magnitude_db = 20 * np.log10(magnitude / reference + 1e-20)
134 # Estimate noise floor (median of spectrum)
135 noise_floor = float(np.median(magnitude))
137 # Find peaks
138 peaks = identify_periodic_components(
139 frequencies,
140 magnitude,
141 n_peaks=n_peaks,
142 min_height=noise_floor * 3,
143 )
145 # Get dominant component
146 if len(peaks) > 0:
147 dominant_frequency = peaks[0][0]
148 dominant_magnitude = peaks[0][1]
149 else:
150 dominant_frequency = None
151 dominant_magnitude = None
153 return JitterSpectrumResult(
154 frequencies=frequencies,
155 magnitude=magnitude,
156 magnitude_db=magnitude_db,
157 dominant_frequency=dominant_frequency,
158 dominant_magnitude=dominant_magnitude,
159 noise_floor=noise_floor,
160 peaks=peaks,
161 )
164def identify_periodic_components(
165 frequencies: NDArray[np.float64],
166 magnitude: NDArray[np.float64],
167 *,
168 n_peaks: int = 10,
169 min_height: float | None = None,
170 min_distance: int = 3,
171) -> list[tuple[float, float]]:
172 """Identify periodic components from jitter spectrum.
174 Finds peaks in the magnitude spectrum that indicate periodic
175 jitter sources.
177 Args:
178 frequencies: Frequency array in Hz.
179 magnitude: Magnitude spectrum.
180 n_peaks: Maximum number of peaks to return.
181 min_height: Minimum peak height (default: 3x median).
182 min_distance: Minimum distance between peaks in bins.
184 Returns:
185 List of (frequency, magnitude) tuples, sorted by magnitude.
186 """
187 if len(magnitude) < 3:
188 return []
190 min_height_value: float
191 min_height_value = float(np.median(magnitude) * 3) if min_height is None else min_height
193 # Find peaks
194 peak_indices, _properties = find_peaks(
195 magnitude,
196 height=min_height_value,
197 distance=min_distance,
198 )
200 if len(peak_indices) == 0:
201 return []
203 # Get peak heights
204 peak_heights = magnitude[peak_indices]
206 # Sort by magnitude (descending)
207 sorted_order = np.argsort(peak_heights)[::-1]
208 sorted_indices = peak_indices[sorted_order][:n_peaks]
210 # Build result list
211 peaks = [(float(frequencies[idx]), float(magnitude[idx])) for idx in sorted_indices]
213 return peaks
216__all__ = [
217 "JitterSpectrumResult",
218 "identify_periodic_components",
219 "jitter_spectrum",
220]