added digital filter to ADC_Result

src/damaris/data/ADC_Result.py

@@ -9,6 +9,7 @@ import numpy
 import sys
 import datetime
 import tables
+from scipy.signal import filtfilt, remez
 #############################################################################
 #                                                                           #
 # Name: Class ADC_Result                                                    #
@@ -147,6 +148,97 @@ class ADC_Result(Resultable, Drawable, DamarisFFT, Signalpath):
         return self.is_clipped
 
 
+    def digital_filter(self, cutoff, numtaps=None):
+        """
+        Apply a zero-phase lowpass FIR filter to all channels.
+
+        Uses scipy.signal.filtfilt (forward-backward filtering) with a Parks-McClellan
+        FIR filter (remez). The forward-backward pass guarantees zero phase distortion
+        and doubles the effective filter order.
+
+        The stopband is automatically set to 0.1 * Nyquist beyond the cutoff. If
+        numtaps is not provided, it is estimated from the transition width to yield
+        approximately 40-50 dB stopband attenuation (typically ~40 taps).
+
+        Parameters
+        ----------
+        cutoff : float
+            Passband edge frequency in Hz. Frequencies above this value are attenuated.
+            Must be a positive number below the Nyquist frequency.
+        numtaps : int, optional
+            Number of FIR taps. If omitted, computed automatically from the transition
+            width.
+
+        Returns
+        -------
+        ADC_Result
+            A new ADC_Result with the filtered data (float64 per channel).
+            The original object is not modified.
+
+        Raises
+        ------
+        RuntimeError
+            If the result contains no data.
+        ValueError
+            If cutoff is invalid for the current sampling rate.
+
+        Examples
+        --------
+        >>> # Lowpass at 1 MHz (taps auto-computed, ~40)
+        ... filtered = adc.digital_filter(cutoff=1e6)
+        >>> # Lowpass at 500 kHz with fixed 63 taps
+        ... filtered = adc.digital_filter(cutoff=500e3, numtaps=63)
+        """
+        if not self.contains_data():
+            raise RuntimeError("digital_filter: no data present")
+
+        nyquist = self.sampling_rate / 2.0
+        cutoff = float(cutoff)
+
+        if cutoff <= 0 or cutoff >= nyquist:
+            raise ValueError(
+                f"cutoff {cutoff} Hz is outside valid range (0, {nyquist})"
+            )
+
+        # Transition width: 0.1 * Nyquist
+        trans_width = 0.1 * nyquist
+        stopband = min(cutoff + trans_width, nyquist * 0.99)
+
+        # Auto-estimate taps if not given
+        if numtaps is None:
+            numtaps = max(21, int(round(8.0 / (stopband - cutoff) * nyquist)))
+        if numtaps % 2 == 0:
+            numtaps += 1
+
+        coeffs = remez(numtaps, [0, cutoff, stopband, nyquist], [1, 0], fs=self.sampling_rate)
+
+        self.lock.acquire()
+        try:
+            tmp_y = []
+            for i in range(self.get_number_of_channels()):
+                data = numpy.asarray(self.y[i], dtype='float64')
+                filtered = filtfilt(
+                    coeffs, [1.0], data,
+                    padtype='odd',
+                    padlen=(numtaps - 1) * 3
+                )
+                tmp_y.append(filtered)
+
+            r = ADC_Result(
+                x=self.x[:],
+                y=tmp_y,
+                index=self.index[:],
+                sampl_freq=self.sampling_rate,
+                desc=self.description,
+                job_id=self.job_id,
+                job_date=self.job_date,
+            )
+        finally:
+            self.lock.release()
+
+        return r
+
+
     def add_sample_space(self, samples):
         "Adds space for n samples, where n can also be negative (deletes space). New space is filled up with \"0\""