some experimenatal data analysis

scripts/experiments/007_generate_fake_analysis_data.py

@@ -67,7 +67,7 @@ def generate_phase_one_data():
     total_lifetimes = 20
     eom_off_lifetime = total_lifetimes * 1 / 2
     fluct_size = 0.05
-    SNR = 0.8
+    SNR = 0.5
 
     params, t, solution = make_params_and_solve(
         total_lifetimes,

scripts/experiments/008_analyze_11_07_transients.py

@@ -0,0 +1,90 @@
+from rabifun.system import *
+from rabifun.plots import *
+from rabifun.utilities import *
+from rabifun.analysis import *
+from ringfit.data import ScanData
+from ringfit.plotting import *
+import functools
+import gc
+
+# %% load data
+path = "../../data/11_07_24/second_signal"
+scan = ScanData.from_dir(path, extension="npz")
+
+
+# %% plot scan
+gc.collect()
+fig = plt.figure("interactive", constrained_layout=True, figsize=(20, 3 * 5))
+fig.clf()
+(ax, ax2, ax_signal, ax_window, ax_spectrum) = fig.subplot_mosaic("AB;CC;DE").values()
+plot_scan(scan, ax=ax_signal, smoothe_output=1e-8, linewidth=0.5)
+
+
+# %% window
+T_step = 0.0002
+t_peak = 0.032201
+N = 100 - 3
+t_scan_peak = t_peak - T_step * N + 4.07e-4  # T * N_steps
+win_length = 5e-05
+window = t_scan_peak - win_length / 2, t_scan_peak + win_length / 2 * 0.8
+
+ax_signal.axvline(t_peak, color="r", linestyle="--")
+ax_signal.axvline(t_scan_peak, color="r", linestyle="--")
+ax_signal.axvspan(
+    *window,
+    color="r",
+    linestyle="--",
+)
+mask = (scan.time > window[0]) & (scan.time < window[1])
+
+ax_window.clear()
+ax_window.plot(scan.time[mask], scan.output[mask], linewidth=0.1)
+
+freq, fft = fourier_transform(
+    scan.time, scan.output, window=window, low_cutoff=2e6, high_cutoff=90e6
+)
+ax.clear()
+ax.plot(freq, np.abs(fft) ** 2)
+
+ringdown_params = RingdownParams(
+    fω_shift=0,
+    mode_window=(0, 4),
+    fΩ_guess=13e6,
+    fδ_guess=0.2 * 13e6,
+    η_guess=0.2e6,
+    absolute_low_cutoff=2e6,
+)
+
+
+peak_info = find_peaks(scan, ringdown_params, window, prominence=0.01)
+peak_info = refine_peaks(peak_info, ringdown_params)
+
+
+plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params)
+print(peak_info.peak_freqs * 1e-6)
+
+
+Ω, ΔΩ, δ, Δδ, ladder = extract_Ω_δ(
+    peak_info,
+    ringdown_params,
+    Ω_threshold=0.1,
+    ladder_threshold=0.1,
+    start_peaks=3,
+    bifurcations=2,
+)
+
+
+for index, type in ladder:
+    freq_index = peak_info.peaks[index]
+    print(
+        type,
+        index,
+        (peak_info.peak_freqs[index] - peak_info.peak_freqs[ladder[0][0]]) / Ω,
+    )
+    ax_spectrum.plot(
+        peak_info.freq[freq_index],
+        peak_info.normalized_power[freq_index],
+        "o" if type.value == StepType.BATH.value else "*",
+        color="C4" if type.value == StepType.BATH.value else "C5",
+        label=type,
+    )

scripts/experiments/009_analyze_12_07_transients.py

@@ -0,0 +1,88 @@
+from rabifun.system import *
+from rabifun.plots import *
+from rabifun.utilities import *
+from rabifun.analysis import *
+from ringfit.data import ScanData
+from ringfit.plotting import *
+import functools
+import gc
+
+# %% load data
+path = "../../data/15_07_24/finite_life_only_bath"
+scan = ScanData.from_dir(path, extension="npz")
+
+
+# %% plot scan
+gc.collect()
+fig = plt.figure("interactive", constrained_layout=True, figsize=(20, 3 * 5))
+fig.clf()
+(ax, ax2, ax_signal, ax_window, ax_spectrum) = fig.subplot_mosaic("AB;CC;DE").values()
+plot_scan(scan, ax=ax_signal, smoothe_output=1e-8, linewidth=0.5)
+
+
+# %% window
+T_step = 0.00010416666666666667
+t_peak = 0.055516 + 2 * T_step - 0.2 * 0.28e-6 - 200 * T_step
+t_scan_peak = t_peak
+win_length = 2.5e-05 * 0.3
+window = t_scan_peak, t_scan_peak + win_length
+
+ax_signal.axvline(t_peak, color="r", linestyle="--")
+ax_signal.axvline(t_scan_peak, color="r", linestyle="--")
+ax_signal.axvspan(
+    *window,
+    color="r",
+    linestyle="--",
+)
+mask = (scan.time > window[0]) & (scan.time < window[1])
+
+ax_window.clear()
+ax_window.plot(scan.time[mask], scan.output[mask], linewidth=0.1)
+
+freq, fft = fourier_transform(
+    scan.time, scan.output, window=window, low_cutoff=3e6, high_cutoff=90e6
+)
+ax.clear()
+ax.plot(freq, np.abs(fft) ** 2)
+
+ringdown_params = RingdownParams(
+    fω_shift=0,
+    mode_window=(0, 4),
+    fΩ_guess=13e6,
+    fδ_guess=0.2 * 13e6,
+    η_guess=4e6,
+    absolute_low_cutoff=4e6,
+)
+
+
+peak_info = find_peaks(scan, ringdown_params, window, prominence=0.1)
+peak_info = refine_peaks(peak_info, ringdown_params)
+
+
+plot_spectrum_and_peak_info(ax_spectrum, peak_info, ringdown_params)
+# %%
+
+Ω, ΔΩ, δ, Δδ, ladder = extract_Ω_δ(
+    peak_info,
+    ringdown_params,
+    Ω_threshold=0.1,
+    ladder_threshold=0.1,
+    start_peaks=3,
+    bifurcations=2,
+)
+
+
+for index, type in ladder:
+    freq_index = peak_info.peaks[index]
+    print(
+        type,
+        index,
+        (peak_info.peak_freqs[index] - peak_info.peak_freqs[ladder[0][0]]) / Ω,
+    )
+    ax_spectrum.plot(
+        peak_info.freq[freq_index],
+        peak_info.normalized_power[freq_index],
+        "o" if type.value == StepType.BATH.value else "*",
+        color="C4" if type.value == StepType.BATH.value else "C5",
+        label=type,
+    )

scripts/ringdown_spectrum_analysis/001_first_exploration.py

@@ -13,7 +13,7 @@ import pickle
 
 # %% load data
 path = "../../data/22_05_24/ringdown_try_2"
-scan = ScanData.from_dir(path)
+scan = ScanData.from_dir(path, extension="npz")
 
 
 # %% Set Window

src/rabifun/analysis.py

@@ -90,10 +90,19 @@ class RingdownParams:
     mode_window: tuple[int, int] = (10, 10)
     """How many FSRs of frequency to consider around :any`fω_shift`."""
 
+    absolute_low_cutoff: float = 0e6
+    """
+    The absolute lowest frequency to consider in Hz.  This has
+    precedence over the ``mode_window``.
+    """
+
     @property
     def low_cutoff(self) -> float:
         """The low cutoff frequency of the ringdown spectrum fft."""
-        return max(self.fω_shift - self.mode_window[0] * self.fΩ_guess, 0)
+        return max(
+            self.fω_shift - self.mode_window[0] * self.fΩ_guess,
+            self.absolute_low_cutoff,
+        )
 
     @property
     def high_cutoff(self) -> float:
@@ -143,7 +152,7 @@ def find_peaks(
     params: RingdownParams,
     window: tuple[float, float],
     prominence: float = 0.005,
-):
+) -> RingdownPeakData:
     """Determine the peaks of the normalized power spectrum of the
     ringdown data.
 
@@ -321,7 +330,7 @@ def extract_Ω_δ(
     # first step: we extract the most common frequency spacing
     all_diff = np.abs(peak_freqs[:, None] - peak_freqs[None, :])
     all_diff = np.triu(all_diff)
-    all_ΔΩ = np.abs(Δpeak_freqs[:, None] ** 2 + Δpeak_freqs[None, :] ** 2)
+    all_ΔΩ = np.sqrt(Δpeak_freqs[:, None] ** 2 + Δpeak_freqs[None, :] ** 2)
 
     bath_mask = (np.abs((all_diff - Ω_guess)) / Ω_guess < Ω_threshold) & (all_diff > 0)
     candidates = all_diff[bath_mask]

src/ringfit/data.py

@@ -56,23 +56,30 @@ class ScanData:
         gc.collect()
 
     @classmethod
-    def from_dir(cls, directory: str | Path, **kwargs):
+    def from_dir(cls, directory: str | Path, extension: str = "npy", **kwargs):
         """Load and parse the oscilloscope data from the
-        ``directory``.  The ``**kwargs`` are passed to the
-        constructor.
+        ``directory``.  The ``extension`` determines the file
+        extension of the data dumps.  The ``**kwargs`` are passed to
+        the constructor.
 
         The directory should contain ``signal_laser.npy``,
         ``signal_outp.npy``, and ``time.npy``.
         """
 
         directory = Path(directory)
-        laserpath = directory / "signal_laser.npy"
+        laserpath = directory / ("signal_laser." + extension)
 
-        output = np.load(directory / "signal_outp.npy")
-        time = np.load(directory / "time.npy")
-        laser = np.load(laserpath) if laserpath.exists() else np.zeros_like(time)
+        output = np.load(directory / ("signal_outp." + extension))
+        time = np.load(directory / ("time." + extension))
+        laser = np.load(laserpath) if laserpath.exists() else None
 
-        return cls(laser, output, time, **kwargs)
+        files = [laser, output, time]
+        for i, file in enumerate(files):
+            if file is not None and isinstance(file, np.lib.npyio.NpzFile):
+                files[i] = file["arr_0"]
+
+        files[0] = laser or np.zeros_like(files[-1])
+        return cls(*files, **kwargs)
 
     @property
     def laser(self):