Coverage for src / tracekit / analyzers / power / ripple.py: 92%
99 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-11 23:04 +0000
1"""Ripple measurement for TraceKit.
3Provides AC ripple analysis for DC power supply outputs.
6Example:
7 >>> from tracekit.analyzers.power.ripple import ripple, ripple_statistics
8 >>> r_pp, r_rms = ripple(dc_output_trace)
9 >>> print(f"Ripple: {r_pp*1e3:.2f} mV pp, {r_rms*1e3:.2f} mV rms")
10"""
12from __future__ import annotations
14import numpy as np
15from scipy import signal
17from tracekit.core.exceptions import AnalysisError
18from tracekit.core.types import WaveformTrace
21def ripple(
22 trace: WaveformTrace,
23 *,
24 dc_coupling: bool = False,
25) -> tuple[float, float]:
26 """Measure AC ripple on a DC signal.
28 Isolates the AC component from the DC offset and measures
29 peak-to-peak and RMS ripple.
31 Args:
32 trace: DC voltage/current waveform with AC ripple.
33 dc_coupling: If True, include DC component in measurement.
34 If False (default), remove DC for pure AC ripple.
36 Returns:
37 Tuple of (ripple_pp, ripple_rms) in signal units.
39 Example:
40 >>> vpp, vrms = ripple(output_voltage)
41 >>> print(f"Ripple: {vpp*1e3:.2f} mV pp, {vrms*1e3:.2f} mV rms")
43 References:
44 IEC 61000-4-7 (power quality)
45 """
46 data = trace.data
48 if dc_coupling:
49 ac_component = data
50 else:
51 # Remove DC (mean)
52 ac_component = data - np.mean(data)
54 ripple_pp = float(np.max(ac_component) - np.min(ac_component))
55 ripple_rms = float(np.sqrt(np.mean(ac_component**2)))
57 return ripple_pp, ripple_rms
60def ripple_percentage(
61 trace: WaveformTrace,
62) -> tuple[float, float]:
63 """Measure ripple as percentage of DC level.
65 Args:
66 trace: DC voltage/current waveform with AC ripple.
68 Returns:
69 Tuple of (ripple_pp_percent, ripple_rms_percent).
71 Example:
72 >>> pp_pct, rms_pct = ripple_percentage(output_voltage)
73 >>> print(f"Ripple: {pp_pct:.2f}% pp, {rms_pct:.2f}% rms")
74 """
75 dc_level = float(np.mean(trace.data))
77 if dc_level == 0: 77 ↛ 78line 77 didn't jump to line 78 because the condition on line 77 was never true
78 return np.nan, np.nan
80 r_pp, r_rms = ripple(trace)
82 return (r_pp / dc_level * 100, r_rms / dc_level * 100)
85def ripple_frequency(
86 trace: WaveformTrace,
87 *,
88 min_frequency: float | None = None,
89 max_frequency: float | None = None,
90) -> float:
91 """Find dominant ripple frequency.
93 Args:
94 trace: DC voltage waveform with AC ripple.
95 min_frequency: Minimum frequency to consider (Hz).
96 max_frequency: Maximum frequency to consider (Hz).
98 Returns:
99 Dominant ripple frequency in Hz.
101 Example:
102 >>> f_ripple = ripple_frequency(output_voltage)
103 >>> print(f"Ripple frequency: {f_ripple/1e3:.2f} kHz")
104 """
105 data = trace.data
106 sample_rate = trace.metadata.sample_rate
108 # Remove DC
109 ac_data = data - np.mean(data)
111 # FFT
112 n = len(ac_data)
113 fft_result = np.abs(np.fft.rfft(ac_data))
114 freqs = np.fft.rfftfreq(n, 1 / sample_rate)
116 # Apply frequency limits
117 freq_mask = np.ones(len(freqs), dtype=bool)
118 if min_frequency is not None:
119 freq_mask &= freqs >= min_frequency
120 if max_frequency is not None:
121 freq_mask &= freqs <= max_frequency
123 # Exclude DC
124 freq_mask[0] = False
126 if not np.any(freq_mask): 126 ↛ 127line 126 didn't jump to line 127 because the condition on line 126 was never true
127 return 0.0
129 # Find peak
130 masked_fft = fft_result.copy()
131 masked_fft[~freq_mask] = 0
132 peak_idx = np.argmax(masked_fft)
134 return float(freqs[peak_idx])
137def ripple_harmonics(
138 trace: WaveformTrace,
139 fundamental_freq: float | None = None,
140 n_harmonics: int = 10,
141) -> dict[int, float]:
142 """Analyze ripple harmonics.
144 Args:
145 trace: DC voltage waveform with AC ripple.
146 fundamental_freq: Fundamental ripple frequency. If None, auto-detect.
147 n_harmonics: Number of harmonics to analyze.
149 Returns:
150 Dictionary mapping harmonic number to amplitude.
152 Example:
153 >>> harmonics = ripple_harmonics(output_voltage)
154 >>> for h, amp in harmonics.items():
155 ... print(f"H{h}: {amp*1e3:.2f} mV")
156 """
157 data = trace.data
158 sample_rate = trace.metadata.sample_rate
160 # Remove DC
161 ac_data = data - np.mean(data)
163 # Find fundamental if not provided
164 if fundamental_freq is None:
165 fundamental_freq = ripple_frequency(trace)
167 if fundamental_freq <= 0:
168 return {}
170 # FFT
171 n = len(ac_data)
172 fft_result = np.abs(np.fft.rfft(ac_data)) * 2 / n # Scale for amplitude
173 freqs = np.fft.rfftfreq(n, 1 / sample_rate)
175 harmonics = {}
176 for h in range(1, n_harmonics + 1):
177 target_freq = h * fundamental_freq
178 # Find closest bin
179 idx = np.argmin(np.abs(freqs - target_freq))
180 if idx < len(fft_result): 180 ↛ 176line 180 didn't jump to line 176 because the condition on line 180 was always true
181 harmonics[h] = float(fft_result[idx])
183 return harmonics
186def ripple_statistics(
187 trace: WaveformTrace,
188) -> dict[str, float]:
189 """Calculate comprehensive ripple statistics.
191 Args:
192 trace: DC voltage waveform with AC ripple.
194 Returns:
195 Dictionary with:
196 - dc_level: DC (mean) level
197 - ripple_pp: Peak-to-peak ripple
198 - ripple_rms: RMS ripple
199 - ripple_pp_percent: Peak-to-peak as % of DC
200 - ripple_rms_percent: RMS as % of DC
201 - ripple_frequency: Dominant ripple frequency
202 - crest_factor: Ripple peak / ripple RMS
204 Example:
205 >>> stats = ripple_statistics(output_voltage)
206 >>> print(f"DC: {stats['dc_level']:.2f} V")
207 >>> print(f"Ripple: {stats['ripple_pp']*1e3:.2f} mV pp")
208 """
209 data = trace.data
210 dc_level = float(np.mean(data))
211 ac_data = data - dc_level
213 r_pp = float(np.max(ac_data) - np.min(ac_data))
214 r_rms = float(np.sqrt(np.mean(ac_data**2)))
215 r_peak = float(np.max(np.abs(ac_data)))
217 crest_factor = r_peak / r_rms if r_rms > 0 else 0.0
219 if dc_level != 0: 219 ↛ 223line 219 didn't jump to line 223 because the condition on line 219 was always true
220 r_pp_pct = r_pp / dc_level * 100
221 r_rms_pct = r_rms / dc_level * 100
222 else:
223 r_pp_pct = np.nan
224 r_rms_pct = np.nan
226 return {
227 "dc_level": dc_level,
228 "ripple_pp": r_pp,
229 "ripple_rms": r_rms,
230 "ripple_pp_percent": r_pp_pct,
231 "ripple_rms_percent": r_rms_pct,
232 "ripple_frequency": ripple_frequency(trace),
233 "crest_factor": crest_factor,
234 }
237def extract_ripple(
238 trace: WaveformTrace,
239 *,
240 high_pass_freq: float | None = None,
241 filter_order: int = 4,
242) -> WaveformTrace:
243 """Extract AC ripple component from DC signal.
245 Args:
246 trace: DC voltage waveform with AC ripple.
247 high_pass_freq: High-pass filter cutoff. If None, uses simple DC removal.
248 filter_order: Order of the Butterworth high-pass filter (default: 4).
249 Higher orders give sharper cutoff but more phase distortion.
251 Returns:
252 Waveform trace containing only the AC ripple component.
254 Raises:
255 AnalysisError: If high_pass_freq exceeds Nyquist frequency
257 Example:
258 >>> ac_ripple = extract_ripple(output_voltage)
259 >>> # Now analyze or plot just the ripple
260 >>> # With custom filter order for sharper cutoff:
261 >>> ac_ripple = extract_ripple(output_voltage, high_pass_freq=10, filter_order=6)
262 """
263 data = trace.data
264 sample_rate = trace.metadata.sample_rate
266 if high_pass_freq is not None:
267 # Use high-pass filter
268 nyquist = sample_rate / 2
269 if high_pass_freq >= nyquist:
270 raise AnalysisError(
271 f"High-pass frequency {high_pass_freq} Hz must be less than "
272 f"Nyquist frequency {nyquist} Hz"
273 )
275 sos = signal.butter(filter_order, high_pass_freq / nyquist, btype="high", output="sos")
276 ac_data = signal.sosfiltfilt(sos, data)
277 else:
278 # Simple DC removal
279 ac_data = data - np.mean(data)
281 return WaveformTrace(
282 data=ac_data.astype(np.float64),
283 metadata=trace.metadata,
284 )
287def ripple_envelope(
288 trace: WaveformTrace,
289 *,
290 method: str = "hilbert",
291 peak_window_size: int | None = None,
292) -> WaveformTrace:
293 """Extract ripple envelope (for amplitude modulation analysis).
295 Args:
296 trace: DC voltage waveform with AC ripple.
297 method: Envelope detection method ("hilbert" or "peak").
298 peak_window_size: Window size for peak envelope detection (samples).
299 If None, defaults to a size that covers approximately one ripple period.
300 Only used when method="peak".
302 Returns:
303 Waveform trace containing the ripple envelope.
305 Raises:
306 AnalysisError: If unknown envelope method specified
308 Example:
309 >>> envelope = ripple_envelope(output_voltage)
310 >>> # Analyze envelope for beat frequencies, etc.
311 >>> # With custom peak window size:
312 >>> envelope = ripple_envelope(output_voltage, method="peak", peak_window_size=200)
313 """
314 # First extract AC component
315 ac_trace = extract_ripple(trace)
316 ac_data = ac_trace.data
318 if method == "hilbert":
319 analytic_signal = signal.hilbert(ac_data)
320 envelope = np.abs(analytic_signal)
321 elif method == "peak":
322 # Simple peak detection
323 from scipy.ndimage import maximum_filter1d
325 # Determine window size
326 if peak_window_size is None: 326 ↛ 338line 326 didn't jump to line 338 because the condition on line 326 was always true
327 # Default: scale to signal frequency if possible
328 # Try to detect ripple frequency and use ~1 period
329 ripple_freq = ripple_frequency(trace)
330 sample_rate = trace.metadata.sample_rate
331 if ripple_freq > 0: 331 ↛ 336line 331 didn't jump to line 336 because the condition on line 331 was always true
332 # Use approximately one period of the ripple
333 peak_window_size = max(10, int(sample_rate / ripple_freq))
334 else:
335 # Fallback: use 1% of signal length, min 10, max 1000
336 peak_window_size = max(10, min(1000, len(ac_data) // 100))
338 envelope = maximum_filter1d(np.abs(ac_data), size=peak_window_size)
339 else:
340 raise AnalysisError(f"Unknown envelope method: {method}")
342 return WaveformTrace(
343 data=envelope.astype(np.float64),
344 metadata=trace.metadata,
345 )
348__all__ = [
349 "extract_ripple",
350 "ripple",
351 "ripple_envelope",
352 "ripple_frequency",
353 "ripple_harmonics",
354 "ripple_percentage",
355 "ripple_statistics",
356]