fix smoothing

scripts/transients_with_blackout.py

@@ -9,13 +9,17 @@ from ringfit.fit import *
 path = "/home/hiro/Documents/org/roam/code/fitting_ringdown/data/09_05_24/Nicely_hybridised_2 2024,05,09, 15h57min00sec/"
 scan = ScanData.from_dir(path, truncation=[0, 50])
 
-fig, ax = plot_scan(scan, smoothe_output=50, normalize=True, laser=False, steps=True)
+# %% interactive
+fig = plt.figure("interactive")
+fig.clf()
+ax, ax_trans = fig.subplots(2, 1)
+plot_scan(scan, smoothe_output=10**2, normalize=True, laser=False, steps=True, ax=ax)
+
 time, output, _ = scan.for_step(10)
 t, o, params, cov, scaled = fit_transient(time, output, window_size=100)
 
-plt.figure()
+ax_trans.plot(t, o)
-plt.plot(t, o)
+ax_trans.plot(t, transient_model(t, *params))
-plt.plot(t, transient_model(t, *params))
+ax_trans.set_title(
-plt.title(
+    f"Transient 2, τ_γ={1/scaled[1] * 10**6}μs ω/2π={scaled[1] / (2*np.pi * 10**3)}kHz"
-    f"Transient 2, γ={scaled[1] / 10**3}kHz ω/2π={scaled[1] / (2*np.pi * 10**3)}kHz"
 )

scripts/transients_without_amplification.py

@@ -6,25 +6,30 @@ from ringfit.data import *
 from ringfit.plotting import *
 from ringfit.fit import *
 
-path = (
+path = "/home/hiro/Documents/org/roam/code/fitting_ringdown/data/08_05_24/characterization_first"
-    "/home/hiro/Documents/org/roam/code/fitting_ringdown/data/08_05_24/nice_transient_2"
+scan = ScanData.from_dir(path)
-)
-scan = ScanData.from_dir(path, truncation=[0, 50])
 
 # %% interactive
-STEPS = [2, 3, 5]
+STEPS = [25, 27, 28]
 fig = plt.figure("interactive")
-ax, *axs = fig.subplots(1, len(STEPS))
+fig.clf()
-plot_scan(scan, smoothe_output=500, normalize=True, laser=True, steps=True, ax=ax)
+ax, *axs = fig.subplots(1, len(STEPS) + 1)
 
+plot_scan(
+    scan, smoothe_output=10 ** (-6), normalize=True, laser=True, steps=True, ax=ax
+)
+
+# %% plot steps
 for ax, STEP in zip(axs, STEPS):
-    time, output, _ = scan.for_step(step=STEP)
+    data = scan.for_step(step=STEP).smoothed(10**-6)
-    t, o, params, cov, scaled = fit_transient(time, output, window_size=100)
 
+    t, o, params, cov, scaled = fit_transient(data.time, data.output)
+
+    ax.cla()
     ax.plot(t, o)
     ax.plot(t, transient_model(t, *params))
     ax.set_title(
-        f"Transient 2, γ={scaled[1] / 10**3:.2f}kHz ({cov[1] / 10**3:.2f}kHz)\n ω/2π={scaled[0] / (2*np.pi * 10**3):.5f}kHz\n step={STEP}"
+        f"Transient {STEP}, γ={1/scaled[1] * 10**6:.2f}μs ({cov[1]/scaled[1]**2 * 10**6:.2f}μs)\n ω/2π={scaled[0] / (2*np.pi * 10**3):.5f}kHz\n step={STEP}"
     )
 
     freq_unit = params[1] / scaled[1]

scripts/weird_oscillations.py

@@ -19,7 +19,7 @@ fig.clf()
 (ax1, *axs) = fig.subplots(nrows=1, ncols=len(STEPS) + 1)
 
 # %% Plot scan
-plot_scan(scan, smoothe_output=100, normalize=True, laser=False, steps=True, ax=ax1)
+plot_scan(scan, smoothe_output=50**2, normalize=True, laser=False, steps=True, ax=ax1)
 
 
 # %% Plot Frequency Fits
@@ -30,7 +30,7 @@ def fit_frequency(step, ax):
     end = int(0.8 * l)
     time = time[begin:end]
     output = output[begin:end]
-    output = smoothe_signal(output, 0.05)
+    output = smoothe_signal(output, 50**2)
     output = shift_and_normalize(output)
     output -= output.mean()
 

src/ringfit/data.py

@@ -6,6 +6,7 @@ import numpy as np
 from pathlib import Path
 import functools
 from . import utils
+import gc
 
 
 class ScanData:
@@ -14,8 +15,8 @@ class ScanData:
         laser: np.ndarray,
         output: np.ndarray,
         time: np.ndarray,
-        truncation: [float, float] = [0, 100],
+        truncation: tuple[float, float] = (0.0, 100.0),
-        sparcity: float = 1,
+        max_frequency: float = np.inf,
     ):
         """
         A class to hold the data from an oscilloscope scan where the
@@ -36,13 +37,23 @@ class ScanData:
 
         length = len(laser)
         begin = int(truncation[0] * length / 100)
-        end = int(truncation[1] * length / 100)
+        end = int(truncation[1] * length / 100) + 1
 
-        every = int(1 / sparcity)
+        self._laser = laser[begin:end]
+        self._output = output[begin:end]
+        self._time = time[begin:end]
 
-        self._laser = laser[begin:end:every]
+        step = time[2] - time[1]
-        self._output = output[begin:end:every]
+
-        self._time = time[begin:end:every]
+        if 1 / step > max_frequency:
+            new_step = 1 / max_frequency
+            index_step = int(new_step / step)
+
+            self._laser = self._laser[::index_step]
+            self._output = self._output[::index_step]
+            self._time = self._time[::index_step]
+
+        gc.collect()
 
     @classmethod
     def from_dir(cls, directory: str | Path, **kwargs):
@@ -79,6 +90,38 @@ class ScanData:
 
         return self._time
 
+    @property
+    def timestep(self):
+        """The time between each sample."""
+
+        return self._time[2] - self._time[1]
+
+    def __len__(self):
+        """The number of samples in the data."""
+        return len(self._laser)
+
+    def smoothed(self, window: float):
+        """
+        Return a smoothed version of the data where the signals are
+        smoothened using a window of size ``window`` (in time units).
+        """
+
+        step = self.timestep
+        return ScanData(
+            self._laser,
+            utils.smoothe_signal(self._output.copy(), window, step),
+            self._time,
+        )
+
+    def sparsified(self, max_frequency: float):
+        """Return a sparsified version of the data where the frequency
+        is limited to ``max_frequency``.
+        """
+
+        return ScanData(
+            self._laser, self._output, self._time, max_frequency=max_frequency
+        )
+
     @functools.cache
     def laser_steps(self, *args, **kwargs) -> np.ndarray:
         """Find the indices of the laser signal ``laser`` where the
@@ -133,8 +176,9 @@ class ScanData:
             raise ValueError("The step must be between 0 and the number of steps.")
 
         padded_steps = [0, *time_steps, len(self._output) - 1]
-        return (
+
-            self._time[padded_steps[step] : padded_steps[step + 1]],
+        return ScanData(
-            self._output[padded_steps[step] : padded_steps[step + 1]],
             self._laser[padded_steps[step] : padded_steps[step + 1]],
+            self._output[padded_steps[step] : padded_steps[step + 1]],
+            self._time[padded_steps[step] : padded_steps[step + 1]],
         )

src/ringfit/fit.py

@@ -15,24 +15,13 @@ def transient_model(t, Δω, γ, amplitude, phase):
     return np.imag(osci)
 
 
-def fit_transient(time: np.ndarray, transient: np.ndarray, window_size: int = 100):
+def fit_transient(time: np.ndarray, transient: np.ndarray):
     """
     Fit a transient signal ``transient`` over ``time`` to a damped
     oscillation model.
-
-    The smoothing window is calculated as the length of the transient
-    divided by the ``window_size``.
     """
 
-    # data_length = len(transient)
-    # begin, end = transient.argmax(), -int(data_length * 0.01)
-
-    # time = time[begin:end]
-    # output_data = transient[begin:end]
-
     output_data = transient
-    window = len(output_data) // window_size
-    output_data = uniform_filter1d(output_data, window)
     output_data -= output_data[0]
     output_data /= abs(output_data).max()
 

src/ringfit/plotting.py

@@ -31,13 +31,19 @@ def plot_scan(
     output=True,
     steps: bool | int = False,
     normalize=False,
-    smoothe_output: bool | int = False,
+    smoothe_output: bool | float = False,
     ax=None,
     **kwargs,
 ):
     if not (laser or output):
         raise ValueError("At least one of 'laser' or 'output' must be True.")
 
+    if smoothe_output:
+        if not isinstance(smoothe_output, float):
+            smoothe_output = 10 ** (-7)
+        data = data.smoothed(smoothe_output)
+        data = data.sparsified(2 / smoothe_output)
+
     time, output_data, laser_data = (
         (data.time, data.output, data.laser)
         if isinstance(steps, bool)
@@ -58,13 +64,6 @@ def plot_scan(
                 output_data.max() - output_data.min()
             )
 
-        if smoothe_output:
-            if not isinstance(smoothe_output, int):
-                smoothe_output = 60
-
-            window = len(output_data) // smoothe_output
-            output_data = uniform_filter1d(output_data, window)
-
         lines = ax.plot(time, output_data, **kwargs)
 
     if isinstance(steps, bool) and steps:

src/ringfit/utils.py

@@ -38,8 +38,11 @@ def shift_and_normalize(array: np.ndarray) -> np.ndarray:
     return shifted / abs(shifted).max()
 
 
-def smoothe_signal(signal: np.ndarray, window_size: float = 0.01) -> np.ndarray:
+def smoothe_signal(
+    signal: np.ndarray, window_size: float = 0.01, time_step: float = 1
+) -> np.ndarray:
     """Smoothe the signal ``signal`` using a uniform filter with a window
-    size of ``window_size * len(signal)``."""
+    size of ``window_size / time_step``."""
-    window = int(len(signal) * window_size)
+
+    window = int(window_size / time_step)
     return uniform_filter1d(signal, window)