Coverage for src / tracekit / filtering / design.py: 99%
160 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-11 23:04 +0000
1"""Filter design functions for TraceKit.
3Provides high-level filter design API with support for Butterworth,
4Chebyshev, Bessel, and Elliptic filter types. Supports automatic
5order calculation from specifications.
8Example:
9 >>> from tracekit.filtering.design import LowPassFilter, design_filter
10 >>> # Simple filter creation
11 >>> lpf = LowPassFilter(cutoff=1e6, sample_rate=10e6, order=4)
12 >>> # Spec-based design
13 >>> filt = design_filter_spec(
14 ... passband=1e6, stopband=2e6,
15 ... passband_ripple=1.0, stopband_atten=40.0,
16 ... sample_rate=10e6
17 ... )
18"""
20from __future__ import annotations
22from typing import Any, Literal
24import numpy as np
25from scipy import signal
27from tracekit.core.exceptions import AnalysisError
28from tracekit.filtering.base import IIRFilter
30FilterType = Literal["butterworth", "chebyshev1", "chebyshev2", "bessel", "elliptic"]
31BandType = Literal["lowpass", "highpass", "bandpass", "bandstop"]
34def design_filter(
35 filter_type: FilterType,
36 cutoff: float | tuple[float, float],
37 sample_rate: float,
38 order: int,
39 btype: BandType = "lowpass",
40 *,
41 ripple_db: float = 1.0,
42 stopband_atten_db: float = 40.0,
43 analog: bool = False,
44 output: Literal["sos", "ba"] = "sos",
45) -> IIRFilter:
46 """Design an IIR filter with specified parameters.
48 Args:
49 filter_type: Type of filter ("butterworth", "chebyshev1", "chebyshev2",
50 "bessel", or "elliptic").
51 cutoff: Cutoff frequency in Hz. For bandpass/bandstop, tuple of (low, high).
52 sample_rate: Sample rate in Hz.
53 order: Filter order.
54 btype: Band type ("lowpass", "highpass", "bandpass", "bandstop").
55 ripple_db: Passband ripple in dB (for Chebyshev and Elliptic).
56 stopband_atten_db: Stopband attenuation in dB (for Chebyshev2 and Elliptic).
57 analog: If True, design analog filter (s-domain). Default is digital (z-domain).
58 output: Output format - "sos" for second-order sections (recommended),
59 "ba" for transfer function polynomials.
61 Returns:
62 IIRFilter object with designed coefficients.
64 Raises:
65 AnalysisError: If cutoff frequency is invalid or filter design fails.
67 Example:
68 >>> lpf = design_filter("butterworth", 1e6, 10e6, order=4)
69 >>> filtered = lpf.apply(trace)
71 References:
72 scipy.signal.iirfilter, butter, cheby1, cheby2, ellip, bessel
73 """
74 # Normalize cutoff frequency
75 if isinstance(cutoff, tuple):
76 Wn = cutoff if analog else (cutoff[0] / (sample_rate / 2), cutoff[1] / (sample_rate / 2))
77 else:
78 Wn = cutoff if analog else cutoff / (sample_rate / 2) # type: ignore[assignment]
80 # Validate normalized frequency
81 if not analog:
82 if isinstance(Wn, tuple):
83 if not (0 < Wn[0] < 1 and 0 < Wn[1] < 1):
84 raise AnalysisError(
85 f"Normalized cutoff must be in (0, 1), got {Wn}. "
86 f"Cutoff {cutoff} Hz must be less than Nyquist {sample_rate / 2} Hz."
87 )
88 elif not 0 < Wn < 1: # type: ignore[unreachable]
89 raise AnalysisError(
90 f"Normalized cutoff must be in (0, 1), got {Wn}. "
91 f"Cutoff {cutoff} Hz must be less than Nyquist {sample_rate / 2} Hz."
92 )
94 # Design filter
95 ftype_map = {
96 "butterworth": "butter",
97 "chebyshev1": "cheby1",
98 "chebyshev2": "cheby2",
99 "bessel": "bessel",
100 "elliptic": "ellip",
101 }
102 ftype = ftype_map.get(filter_type)
103 if ftype is None:
104 raise AnalysisError(f"Unknown filter type: {filter_type}")
106 try:
107 if filter_type == "butterworth":
108 if output == "sos":
109 sos = signal.butter(order, Wn, btype=btype, analog=analog, output="sos")
110 return IIRFilter(sample_rate=sample_rate, sos=sos)
111 else:
112 b, a = signal.butter(order, Wn, btype=btype, analog=analog, output="ba")
113 return IIRFilter(sample_rate=sample_rate, ba=(b, a))
115 elif filter_type == "chebyshev1":
116 if output == "sos":
117 sos = signal.cheby1(order, ripple_db, Wn, btype=btype, analog=analog, output="sos")
118 return IIRFilter(sample_rate=sample_rate, sos=sos)
119 else:
120 b, a = signal.cheby1(order, ripple_db, Wn, btype=btype, analog=analog, output="ba")
121 return IIRFilter(sample_rate=sample_rate, ba=(b, a))
123 elif filter_type == "chebyshev2":
124 if output == "sos":
125 sos = signal.cheby2(
126 order,
127 stopband_atten_db,
128 Wn,
129 btype=btype,
130 analog=analog,
131 output="sos",
132 )
133 return IIRFilter(sample_rate=sample_rate, sos=sos)
134 else:
135 b, a = signal.cheby2(
136 order,
137 stopband_atten_db,
138 Wn,
139 btype=btype,
140 analog=analog,
141 output="ba",
142 )
143 return IIRFilter(sample_rate=sample_rate, ba=(b, a))
145 elif filter_type == "bessel":
146 if output == "sos":
147 sos = signal.bessel(
148 order, Wn, btype=btype, analog=analog, output="sos", norm="phase"
149 )
150 return IIRFilter(sample_rate=sample_rate, sos=sos)
151 else:
152 b, a = signal.bessel(
153 order, Wn, btype=btype, analog=analog, output="ba", norm="phase"
154 )
155 return IIRFilter(sample_rate=sample_rate, ba=(b, a))
157 elif filter_type == "elliptic": 157 ↛ 182line 157 didn't jump to line 182 because the condition on line 157 was always true
158 if output == "sos":
159 sos = signal.ellip(
160 order,
161 ripple_db,
162 stopband_atten_db,
163 Wn,
164 btype=btype,
165 analog=analog,
166 output="sos",
167 )
168 return IIRFilter(sample_rate=sample_rate, sos=sos)
169 else:
170 b, a = signal.ellip(
171 order,
172 ripple_db,
173 stopband_atten_db,
174 Wn,
175 btype=btype,
176 analog=analog,
177 output="ba",
178 )
179 return IIRFilter(sample_rate=sample_rate, ba=(b, a))
181 else:
182 raise AnalysisError(f"Unsupported filter type: {filter_type}")
184 except Exception as e:
185 raise AnalysisError(f"Filter design failed: {e}") from e
188def design_filter_spec(
189 passband: float | tuple[float, float],
190 stopband: float | tuple[float, float],
191 sample_rate: float,
192 passband_ripple: float = 1.0,
193 stopband_atten: float = 40.0,
194 *,
195 filter_type: FilterType = "elliptic",
196 analog: bool = False,
197) -> IIRFilter:
198 """Design filter from passband/stopband specifications.
200 Automatically computes the minimum filter order required to meet
201 the specifications.
203 Args:
204 passband: Passband edge frequency in Hz. Tuple for bandpass/bandstop.
205 stopband: Stopband edge frequency in Hz. Tuple for bandpass/bandstop.
206 sample_rate: Sample rate in Hz.
207 passband_ripple: Maximum passband ripple in dB.
208 stopband_atten: Minimum stopband attenuation in dB.
209 filter_type: Filter type to design.
210 analog: If True, design analog filter.
212 Returns:
213 IIRFilter object with minimum-order design.
215 Raises:
216 AnalysisError: If filter order cannot be determined.
218 Example:
219 >>> # Design a filter with 1MHz passband, 2MHz stopband, 40dB rejection
220 >>> filt = design_filter_spec(
221 ... passband=1e6, stopband=2e6,
222 ... passband_ripple=1.0, stopband_atten=40.0,
223 ... sample_rate=10e6
224 ... )
225 """
226 # Normalize frequencies
227 if isinstance(passband, tuple):
228 wp = (passband[0] / (sample_rate / 2), passband[1] / (sample_rate / 2))
229 ws = (stopband[0] / (sample_rate / 2), stopband[1] / (sample_rate / 2)) # type: ignore[index]
230 else:
231 wp = passband / (sample_rate / 2) # type: ignore[assignment]
232 ws = stopband / (sample_rate / 2) # type: ignore[assignment, operator]
234 # Determine band type
235 if isinstance(passband, tuple):
236 # Bandpass or bandstop
237 if passband[0] < stopband[0]: # type: ignore[index]
238 btype: BandType = "bandstop"
239 else:
240 btype = "bandpass"
241 # Lowpass or highpass
242 elif passband < stopband: # type: ignore[operator]
243 btype = "lowpass"
244 else:
245 btype = "highpass"
247 # Compute minimum order
248 try:
249 if filter_type == "butterworth":
250 order, Wn = signal.buttord(wp, ws, passband_ripple, stopband_atten, analog=analog)
251 elif filter_type == "chebyshev1":
252 order, Wn = signal.cheb1ord(wp, ws, passband_ripple, stopband_atten, analog=analog)
253 elif filter_type == "chebyshev2":
254 order, Wn = signal.cheb2ord(wp, ws, passband_ripple, stopband_atten, analog=analog)
255 elif filter_type == "elliptic":
256 order, Wn = signal.ellipord(wp, ws, passband_ripple, stopband_atten, analog=analog)
257 else:
258 # Bessel doesn't have an ord function, estimate based on Butterworth
259 order, Wn = signal.buttord(wp, ws, passband_ripple, stopband_atten, analog=analog)
261 except Exception as e:
262 raise AnalysisError(f"Could not determine filter order: {e}") from e
264 # Design with computed order
265 cutoff = (
266 tuple(w * sample_rate / 2 for w in Wn)
267 if isinstance(Wn, np.ndarray)
268 else Wn * sample_rate / 2
269 )
271 return design_filter(
272 filter_type=filter_type,
273 cutoff=cutoff,
274 sample_rate=sample_rate,
275 order=int(order),
276 btype=btype,
277 ripple_db=passband_ripple,
278 stopband_atten_db=stopband_atten,
279 analog=analog,
280 )
283# =============================================================================
284# Convenience Filter Classes
285# =============================================================================
288class LowPassFilter(IIRFilter):
289 """Low-pass Butterworth filter.
291 Convenient class for creating low-pass filters with sensible defaults.
293 Example:
294 >>> lpf = LowPassFilter(cutoff=1e6, sample_rate=10e6, order=4)
295 >>> filtered = lpf.apply(trace)
296 """
298 def __init__(
299 self,
300 cutoff: float,
301 sample_rate: float,
302 order: int = 4,
303 *,
304 filter_type: FilterType = "butterworth",
305 ripple_db: float = 1.0,
306 stopband_atten_db: float = 40.0,
307 ) -> None:
308 """Initialize low-pass filter.
310 Args:
311 cutoff: Cutoff frequency in Hz (-3dB point for Butterworth).
312 sample_rate: Sample rate in Hz.
313 order: Filter order.
314 filter_type: Type of filter to use.
315 ripple_db: Passband ripple for Chebyshev/Elliptic.
316 stopband_atten_db: Stopband attenuation for Chebyshev2/Elliptic.
317 """
318 filt = design_filter(
319 filter_type=filter_type,
320 cutoff=cutoff,
321 sample_rate=sample_rate,
322 order=order,
323 btype="lowpass",
324 ripple_db=ripple_db,
325 stopband_atten_db=stopband_atten_db,
326 )
327 super().__init__(sample_rate=sample_rate, sos=filt.sos)
328 self._cutoff = cutoff
329 self._filter_type = filter_type
331 @property
332 def cutoff(self) -> float:
333 """Cutoff frequency in Hz."""
334 return self._cutoff
337class HighPassFilter(IIRFilter):
338 """High-pass Butterworth filter.
340 Example:
341 >>> hpf = HighPassFilter(cutoff=1e3, sample_rate=100e3, order=4)
342 >>> filtered = hpf.apply(trace)
343 """
345 def __init__(
346 self,
347 cutoff: float,
348 sample_rate: float,
349 order: int = 4,
350 *,
351 filter_type: FilterType = "butterworth",
352 ripple_db: float = 1.0,
353 stopband_atten_db: float = 40.0,
354 ) -> None:
355 """Initialize high-pass filter."""
356 filt = design_filter(
357 filter_type=filter_type,
358 cutoff=cutoff,
359 sample_rate=sample_rate,
360 order=order,
361 btype="highpass",
362 ripple_db=ripple_db,
363 stopband_atten_db=stopband_atten_db,
364 )
365 super().__init__(sample_rate=sample_rate, sos=filt.sos)
366 self._cutoff = cutoff
368 @property
369 def cutoff(self) -> float:
370 """Cutoff frequency in Hz."""
371 return self._cutoff
374class BandPassFilter(IIRFilter):
375 """Band-pass filter.
377 Example:
378 >>> bpf = BandPassFilter(low=1e3, high=10e3, sample_rate=100e3, order=4)
379 >>> filtered = bpf.apply(trace)
380 """
382 def __init__(
383 self,
384 low: float,
385 high: float,
386 sample_rate: float,
387 order: int = 4,
388 *,
389 filter_type: FilterType = "butterworth",
390 ripple_db: float = 1.0,
391 stopband_atten_db: float = 40.0,
392 ) -> None:
393 """Initialize band-pass filter."""
394 filt = design_filter(
395 filter_type=filter_type,
396 cutoff=(low, high),
397 sample_rate=sample_rate,
398 order=order,
399 btype="bandpass",
400 ripple_db=ripple_db,
401 stopband_atten_db=stopband_atten_db,
402 )
403 super().__init__(sample_rate=sample_rate, sos=filt.sos)
404 self._low = low
405 self._high = high
407 @property
408 def passband(self) -> tuple[float, float]:
409 """Passband frequencies (low, high) in Hz."""
410 return (self._low, self._high)
413class BandStopFilter(IIRFilter):
414 """Band-stop (notch) filter.
416 Example:
417 >>> bsf = BandStopFilter(low=50, high=60, sample_rate=1000, order=4)
418 >>> filtered = bsf.apply(trace) # Remove 50-60 Hz interference
419 """
421 def __init__(
422 self,
423 low: float,
424 high: float,
425 sample_rate: float,
426 order: int = 4,
427 *,
428 filter_type: FilterType = "butterworth",
429 ripple_db: float = 1.0,
430 stopband_atten_db: float = 40.0,
431 ) -> None:
432 """Initialize band-stop filter."""
433 filt = design_filter(
434 filter_type=filter_type,
435 cutoff=(low, high),
436 sample_rate=sample_rate,
437 order=order,
438 btype="bandstop",
439 ripple_db=ripple_db,
440 stopband_atten_db=stopband_atten_db,
441 )
442 super().__init__(sample_rate=sample_rate, sos=filt.sos)
443 self._low = low
444 self._high = high
446 @property
447 def stopband(self) -> tuple[float, float]:
448 """Stopband frequencies (low, high) in Hz."""
449 return (self._low, self._high)
452# =============================================================================
453# Filter Type Classes
454# =============================================================================
457class ButterworthFilter(IIRFilter):
458 """Butterworth filter with maximally flat passband.
460 Example:
461 >>> filt = ButterworthFilter(cutoff=1e6, sample_rate=10e6, order=4, btype="lowpass")
462 """
464 def __init__(
465 self,
466 cutoff: float | tuple[float, float],
467 sample_rate: float,
468 order: int = 4,
469 btype: BandType = "lowpass",
470 ) -> None:
471 filt = design_filter("butterworth", cutoff, sample_rate, order, btype)
472 super().__init__(sample_rate=sample_rate, sos=filt.sos)
475class ChebyshevType1Filter(IIRFilter):
476 """Chebyshev Type I filter with passband ripple.
478 Example:
479 >>> filt = ChebyshevType1Filter(cutoff=1e6, sample_rate=10e6, order=4, ripple_db=0.5)
480 """
482 def __init__(
483 self,
484 cutoff: float | tuple[float, float],
485 sample_rate: float,
486 order: int = 4,
487 btype: BandType = "lowpass",
488 ripple_db: float = 1.0,
489 ) -> None:
490 filt = design_filter("chebyshev1", cutoff, sample_rate, order, btype, ripple_db=ripple_db)
491 super().__init__(sample_rate=sample_rate, sos=filt.sos)
494class ChebyshevType2Filter(IIRFilter):
495 """Chebyshev Type II filter with stopband ripple.
497 Example:
498 >>> filt = ChebyshevType2Filter(cutoff=1e6, sample_rate=10e6, order=4, stopband_atten_db=40)
499 """
501 def __init__(
502 self,
503 cutoff: float | tuple[float, float],
504 sample_rate: float,
505 order: int = 4,
506 btype: BandType = "lowpass",
507 stopband_atten_db: float = 40.0,
508 ) -> None:
509 filt = design_filter(
510 "chebyshev2",
511 cutoff,
512 sample_rate,
513 order,
514 btype,
515 stopband_atten_db=stopband_atten_db,
516 )
517 super().__init__(sample_rate=sample_rate, sos=filt.sos)
520class BesselFilter(IIRFilter):
521 """Bessel filter with maximally flat group delay.
523 Best for preserving waveform shape during filtering.
525 Example:
526 >>> filt = BesselFilter(cutoff=1e6, sample_rate=10e6, order=4)
527 """
529 def __init__(
530 self,
531 cutoff: float | tuple[float, float],
532 sample_rate: float,
533 order: int = 4,
534 btype: BandType = "lowpass",
535 ) -> None:
536 filt = design_filter("bessel", cutoff, sample_rate, order, btype)
537 super().__init__(sample_rate=sample_rate, sos=filt.sos)
540class EllipticFilter(IIRFilter):
541 """Elliptic (Cauer) filter with equiripple passband and stopband.
543 Provides the sharpest transition band for a given order.
545 Example:
546 >>> filt = EllipticFilter(cutoff=1e6, sample_rate=10e6, order=4,
547 ... ripple_db=0.5, stopband_atten_db=60)
548 """
550 def __init__(
551 self,
552 cutoff: float | tuple[float, float],
553 sample_rate: float,
554 order: int = 4,
555 btype: BandType = "lowpass",
556 ripple_db: float = 1.0,
557 stopband_atten_db: float = 40.0,
558 ) -> None:
559 filt = design_filter(
560 "elliptic",
561 cutoff,
562 sample_rate,
563 order,
564 btype,
565 ripple_db=ripple_db,
566 stopband_atten_db=stopband_atten_db,
567 )
568 super().__init__(sample_rate=sample_rate, sos=filt.sos)
571def suggest_filter_type(
572 transition_bandwidth: float,
573 passband_ripple_db: float,
574 stopband_atten_db: float,
575) -> FilterType:
576 """Suggest best filter type based on requirements.
578 Recommends filter type based on design tradeoffs:
579 - Butterworth: Maximally flat passband, moderate rolloff
580 - Chebyshev1: Faster rolloff, passband ripple
581 - Chebyshev2: Faster rolloff, stopband ripple
582 - Elliptic: Sharpest rolloff, both passband and stopband ripple
583 - Bessel: Linear phase, slowest rolloff (for waveform preservation)
585 Args:
586 transition_bandwidth: Normalized transition bandwidth (stopband - passband) / sample_rate.
587 passband_ripple_db: Acceptable passband ripple in dB.
588 stopband_atten_db: Required stopband attenuation in dB.
590 Returns:
591 Recommended filter type.
593 Example:
594 >>> # Sharp transition, can tolerate some ripple
595 >>> ftype = suggest_filter_type(
596 ... transition_bandwidth=0.1,
597 ... passband_ripple_db=0.5,
598 ... stopband_atten_db=60.0
599 ... )
600 >>> print(ftype) # 'elliptic'
602 References:
603 API-020: Filter Design Auto-Order
604 """
605 # For very sharp transitions with ripple tolerance, use elliptic
606 if transition_bandwidth < 0.15 and passband_ripple_db >= 0.1:
607 return "elliptic"
609 # For moderate sharpness with low passband ripple, use Chebyshev2
610 if transition_bandwidth < 0.2 and passband_ripple_db < 0.1:
611 return "chebyshev2"
613 # For moderate sharpness with ripple tolerance, use Chebyshev1
614 if transition_bandwidth < 0.3 and passband_ripple_db >= 0.1:
615 return "chebyshev1"
617 # For phase linearity (waveform preservation), use Bessel
618 if stopband_atten_db < 40.0:
619 return "bessel"
621 # Default to Butterworth for balanced performance
622 return "butterworth"
625def auto_design_filter(
626 passband: float | tuple[float, float],
627 stopband: float | tuple[float, float],
628 sample_rate: float,
629 *,
630 passband_ripple_db: float = 1.0,
631 stopband_atten_db: float = 40.0,
632 suggest_type: bool = True,
633) -> tuple[IIRFilter, dict[str, Any]]:
634 """Automatically design optimal filter from specifications.
636 Automatically computes filter order and optionally suggests the best
637 filter type based on transition band and ripple requirements.
639 Args:
640 passband: Passband edge frequency in Hz. Tuple for bandpass/bandstop.
641 stopband: Stopband edge frequency in Hz. Tuple for bandpass/bandstop.
642 sample_rate: Sample rate in Hz.
643 passband_ripple_db: Maximum passband ripple in dB (default: 1.0).
644 stopband_atten_db: Minimum stopband attenuation in dB (default: 40.0).
645 suggest_type: If True, automatically suggest filter type (default: True).
647 Returns:
648 Tuple of (IIRFilter, design_info_dict).
649 design_info_dict contains: filter_type, order, cutoff, transition_bandwidth.
651 Example:
652 >>> # Automatic filter design with type suggestion
653 >>> filt, info = auto_design_filter(
654 ... passband=1e6,
655 ... stopband=1.5e6,
656 ... sample_rate=10e6,
657 ... stopband_atten_db=60.0
658 ... )
659 >>> print(f"Designed {info['filter_type']} filter with order {info['order']}")
660 >>> filtered = filt.apply(trace)
662 References:
663 API-020: Filter Design Auto-Order
664 """
665 # Compute transition bandwidth
666 if isinstance(passband, tuple):
667 # Bandpass/bandstop - use average
668 transition_bw = (
669 abs(stopband[0] - passband[0]) + abs(stopband[1] - passband[1]) # type: ignore[index]
670 ) / 2.0
671 else:
672 transition_bw = abs(stopband - passband) # type: ignore[operator]
674 normalized_transition = transition_bw / sample_rate
676 # Suggest filter type if requested
677 if suggest_type:
678 filter_type = suggest_filter_type(
679 transition_bandwidth=normalized_transition,
680 passband_ripple_db=passband_ripple_db,
681 stopband_atten_db=stopband_atten_db,
682 )
683 else:
684 filter_type = "butterworth"
686 # Design filter with auto-order computation
687 filt = design_filter_spec(
688 passband=passband,
689 stopband=stopband,
690 sample_rate=sample_rate,
691 passband_ripple=passband_ripple_db,
692 stopband_atten=stopband_atten_db,
693 filter_type=filter_type,
694 )
696 # Extract design info
697 design_info = {
698 "filter_type": filter_type,
699 "order": filt.sos.shape[0] * 2 if filt.sos is not None else 0,
700 "cutoff": passband,
701 "transition_bandwidth": transition_bw,
702 "passband_ripple_db": passband_ripple_db,
703 "stopband_atten_db": stopband_atten_db,
704 }
706 return filt, design_info
709__all__ = [
710 "BandPassFilter",
711 "BandStopFilter",
712 "BandType",
713 "BesselFilter",
714 "ButterworthFilter",
715 "ChebyshevType1Filter",
716 "ChebyshevType2Filter",
717 "EllipticFilter",
718 "FilterType",
719 "HighPassFilter",
720 "LowPassFilter",
721 "auto_design_filter",
722 "design_filter",
723 "design_filter_spec",
724 "suggest_filter_type",
725]