STUCK: look at first ringdown data

poetry.lock

@@ -701,6 +701,24 @@ files = [
     {file = "nest_asyncio-1.6.0.tar.gz", hash = "sha256:6f172d5449aca15afd6c646851f4e31e02c598d553a667e38cafa997cfec55fe"},
 ]
 
+[[package]]
+name = "networkx"
+version = "3.3"
+description = "Python package for creating and manipulating graphs and networks"
+optional = false
+python-versions = ">=3.10"
+files = [
+    {file = "networkx-3.3-py3-none-any.whl", hash = "sha256:28575580c6ebdaf4505b22c6256a2b9de86b316dc63ba9e93abde3d78dfdbcf2"},
+    {file = "networkx-3.3.tar.gz", hash = "sha256:0c127d8b2f4865f59ae9cb8aafcd60b5c70f3241ebd66f7defad7c4ab90126c9"},
+]
+
+[package.extras]
+default = ["matplotlib (>=3.6)", "numpy (>=1.23)", "pandas (>=1.4)", "scipy (>=1.9,!=1.11.0,!=1.11.1)"]
+developer = ["changelist (==0.5)", "mypy (>=1.1)", "pre-commit (>=3.2)", "rtoml"]
+doc = ["myst-nb (>=1.0)", "numpydoc (>=1.7)", "pillow (>=9.4)", "pydata-sphinx-theme (>=0.14)", "sphinx (>=7)", "sphinx-gallery (>=0.14)", "texext (>=0.6.7)"]
+extra = ["lxml (>=4.6)", "pydot (>=2.0)", "pygraphviz (>=1.12)", "sympy (>=1.10)"]
+test = ["pytest (>=7.2)", "pytest-cov (>=4.0)"]
+
 [[package]]
 name = "numpy"
 version = "1.26.4"
@@ -1284,4 +1302,4 @@ files = [
 [metadata]
 lock-version = "2.0"
 python-versions = "^3.11"
-content-hash = "145e0c6ca0ebf2f1f6493a208a0a08871fe361b4458d9a7d4d32d83930aa0f4d"
+content-hash = "7cbb6668719eef0958556b0c35997d0c18a2611c0634994da0a3ad2bc00e8fd6"

pyproject.toml

@@ -10,6 +10,7 @@ python = "^3.11"
 scipy = "^1.13.0"
 click = "^8.1.7"
 matplotlib = "^3.8.4"
+networkx = "^3.3"
 
 
 [tool.poetry.group.dev.dependencies]

scripts/experiments/*python:main*

@@ -0,0 +1,151 @@
+Python 3.11.9 (main, Apr  2 2024, 08:25:04) [GCC 13.2.0]
+Type 'copyright', 'credits' or 'license' for more information
+IPython 8.24.0 -- An enhanced Interactive Python. Type '?' for help.
+
+Python 3.11.9 (main, Apr  2 2024, 08:25:04) [GCC 13.2.0]
+Type 'copyright', 'credits' or 'license' for more information
+IPython 8.24.0 -- An enhanced Interactive Python. Type '?' for help.
+
+In [1]: 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [2]: %run -i /home/hiro/Documents/org/roam/code/fitting_ringdown/scripts/experiments/003_experimental_ringdown_fft.py # load script buffer
+
+In [3]: 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [3]: %run -i /tmp/python-vterm1EdyiK.py # Fourier
+
+In [4]: %run -i /tmp/python-vtermHclYEt.py # Fourier
+
+In [5]: %run -i /tmp/python-vtermUXKyDJ.py # CSV
+
+In [6]: %run -i /tmp/python-vtermIjaHIZ.py # Fourier
+
+In [7]: %run -i /tmp/python-vtermHYpcEb.py # Fourier
+
+In [8]: ax3 = ax.twinx()
+
+In [9]: ax3.plot(freq, np.gradient(np.angle(fft)), linewidth=0.1, color="red", zorder=-10)
+Out[9]: [<matplotlib.lines.Line2D at 0x7f83e3008dd0>]
+
+In [10]: %run -i /tmp/python-vtermIhW5K3.py # peaks
+
+In [11]: %run -i /tmp/python-vterm3zwD7E.py # peaks
+
+In [12]: %run -i /tmp/python-vtermJ1uevg.py # peaks
+---------------------------------------------------------------------------
+ValueError                                Traceback (most recent call last)
+File /tmp/python-vtermJ1uevg.py:4
+      1 # %% peaks
+      2 freq_step = freq[1] - freq[0]
+----> 4 peaks, peak_info = scipy.signal.find_peaks(np.abs(fft) ** 2, distance=1e4 / freq_step)
+      5 ax.plot(freq[peaks], np.abs(fft[peaks]) ** 2, "x")
+
+File /nix/store/dbxrivq1l53hnfjsp491yrjnzy651sqg-python3.11-scipy-1.13.0/lib/python3.11/site-packages/scipy/signal/_peak_finding.py:941, in find_peaks(x, height, threshol
+d, distance, prominence, width, wlen, rel_height, plateau_size)
+    939 x = _arg_x_as_expected(x)
+    940 if distance is not None and distance < 1:
+--> 941     raise ValueError('`distance` must be greater or equal to 1')
+    943 peaks, left_edges, right_edges = _local_maxima_1d(x)
+    944 properties = {}
+
+ValueError: `distance` must be greater or equal to 1
+
+In [13]: %run -i /tmp/python-vtermZ1ZHBs.py # peaks
+---------------------------------------------------------------------------
+ValueError                                Traceback (most recent call last)
+File /tmp/python-vtermZ1ZHBs.py:4
+      1 # %% peaks
+      2 freq_step = freq[1] - freq[0]
+----> 4 peaks, peak_info = scipy.signal.find_peaks(np.abs(fft) ** 2, distance=1e4 / freq_step)
+      5 ax.plot(freq[peaks], np.abs(fft[peaks]) ** 2, "x")
+
+File /nix/store/dbxrivq1l53hnfjsp491yrjnzy651sqg-python3.11-scipy-1.13.0/lib/python3.11/site-packages/scipy/signal/_peak_finding.py:941, in find_peaks(x, height, threshol
+d, distance, prominence, width, wlen, rel_height, plateau_size)
+    939 x = _arg_x_as_expected(x)
+    940 if distance is not None and distance < 1:
+--> 941     raise ValueError('`distance` must be greater or equal to 1')
+    943 peaks, left_edges, right_edges = _local_maxima_1d(x)
+    944 properties = {}
+
+ValueError: `distance` must be greater or equal to 1
+
+In [14]: %run -i /tmp/python-vterm220t3Y.py # peaks
+
+In [15]: %run -i /tmp/python-vtermoezkoM.py # Fourier
+
+In [16]: %run -i /tmp/python-vterm3uhE2k.py # plot
+

scripts/experiments/003_experimental_ringdown_fft.py

@@ -5,22 +5,147 @@ from ringfit.plotting import *
 from ringfit.fit import *
 from rabifun.analysis import *
 import numpy as np
-
+import scipy
+from collections import OrderedDict
+import networkx as nx
+from functools import reduce
 
 # %% load data
-path = "../../data/08_05_24/characterization_first"
+path = "../../data/22_05_24/ringdown_with_hybridized_modes"
 scan = ScanData.from_dir(path)
 
 
 # %% Fourier
+
+# window = (0.027751370026589985, 0.027751370026589985 + 0.00001 / 2)
+window = tuple(
+    np.array([0.03075207891902308, 0.03075207891902308 + 0.00001]) - 1e-3 - 0.82e-6
+)
 freq, fft = fourier_transform(
-    scan.time, scan.output, low_cutoff=1000, high_cutoff=10**8
+    scan.time, scan.output, window=window, low_cutoff=20e5, high_cutoff=100e6
 )
 
+
+# %% plot
+gc.collect()
 fig = plt.figure("interactive")
 fig.clf()
 ax, ax2 = fig.subplots(1, 2)
-ax.set_yscale("log")
+# ax.set_yscale("log")
-ax.plot(freq * 10 ** (-6), np.abs(fft) ** 2)
+ax.plot(freq, np.abs(fft))
-plot_scan(scan, ax=ax2, linewidth=0.1)
+ax3 = ax.twinx()
-ax2.set_xlim(0.02, 0.020001)
+ax.plot(freq, np.abs(fft.real), linewidth=1, color="red")
+
+# ax.plot(freq, fft.imag, linewidth=1, color="green")
+# ax3.plot(
+#     freq,
+#     np.gradient(np.unwrap(np.angle(fft) + np.pi, period=2 * np.pi)),
+#     linestyle="--",
+# )
+
+ax2.set_xlim(*window)
+
+
+# plot_scan(scan, ax=ax2, linewidth=0.5, smoothe_output=1e-8)
+
+freq_step = freq[1] - freq[0]
+
+peaks, peak_info = scipy.signal.find_peaks(
+    np.abs(fft) ** 2, distance=2e6 / freq_step, prominence=3000
+)
+
+anglegrad = np.gradient(np.unwrap(np.angle(fft) + np.pi, period=2 * np.pi))
+neg_peaks = peaks[anglegrad[peaks] < 0]
+pos_peaks = peaks[anglegrad[peaks] > 0]
+ax.plot(freq[neg_peaks], np.abs(fft[neg_peaks]), "x")
+ax.plot(freq[pos_peaks], np.abs(fft[pos_peaks]), "o")
+
+# phase_detuning = np.angle(fft[peaks])
+
+
+def extract_peak(index, width, sign=1):
+    return sign * freq[index - width : index + width], np.abs(
+        fft[index - width : index + width]
+    )
+
+
+# for peak in neg_peaks:
+#     ax2.plot(*extract_peak(peak, 10, -1))
+# for peak in pos_peaks:
+#     ax2.plot(*extract_peak(peak, 10, 1))
+
+
+Ω_guess = 13e6
+δ_guess = 2.69e6
+
+N = np.arange(1, 3)
+possibilieties = np.concatenate([[2 * δ_guess], N * Ω_guess, N * Ω_guess - δ_guess])
+
+
+abs(
+    np.array(
+        [Ω_guess, 2 * δ_guess, Ω_guess - δ_guess, 2 * Ω_guess, 2 * Ω_guess - δ_guess]
+    )
+)
+
+mode_freqs = freq[peaks]
+
+final_freqs = [mode_freqs[0]]
+
+
+all_diffs = np.abs(
+    np.abs(mode_freqs[:, None] - mode_freqs[None, :])[:, :, None]
+    - possibilieties[None, None, :]
+)
+
+all_diffs[all_diffs == 0] = np.inf
+all_diffs[all_diffs > δ_guess / 2] = np.inf
+matches = np.asarray(all_diffs < np.inf).nonzero()
+
+pairs = np.array(list(zip(*matches, all_diffs[matches])), dtype=int)
+
+relationships = nx.DiGraph()
+for node, peak in enumerate(peaks):
+    relationships.add_node(node, freq=freq[peak])
+
+for left, right, relationship, diff in pairs:
+    if freq[left] > freq[right]:
+        left, right = right, left
+
+    if (
+        not relationships.has_edge(left, right)
+        or relationships[left][right]["weight"] > diff * 1e-6
+    ):
+        relationships.add_edge(
+            left, right, weight=diff * 1e-6, type=relationship, freqdis=right - left
+        )
+
+
+pos = {}
+for node, peak in enumerate(peaks):
+    pos[node] = [freq[peak], abs(fft[peak])]
+# nx.draw(relationships, pos, with_labels=True)
+# nx.draw_networkx_edge_labels(relationships, pos)
+
+# %%
+cycle = nx.find_cycle(relationships, orientation="ignore")
+
+total = 0
+for s, t, direction in cycle:
+    difference = possibilieties[relationships[s][t]["type"]]
+    if direction == "reverse":
+        difference *= -1
+    total += difference
+
+# %%
+relationships.remove_nodes_from(list(nx.isolates(relationships)))
+
+spectrum = list(
+    sorted(nx.all_simple_paths(relationships, 1, 9), key=lambda x: -len(x))
+)[0]
+
+for s, t in zip(spectrum[:-1], spectrum[1:]):
+    print(s, relationships[s][t])
+
+for node in spectrum:
+    plt.plot(freq[node], abs(fft[node]), "*")

src/rabifun/analysis.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy
 
 
 def fourier_transform(
@@ -20,9 +21,9 @@ def fourier_transform(
     if window:
         mask = (window[1] > t) & (t > window[0])
         t = t[mask]
-        signal = signal[mask]
+        signal = signal[mask]  # * scipy.signal.windows.hamming(len(t))
 
-    freq = np.fft.rfftfreq(len(t), t[1] - t[0])
+    freq = np.fft.rfftfreq(len(t), t[2] - t[1])
     fft = np.fft.rfft(signal)
 
     mask = (freq > low_cutoff) & (freq < high_cutoff)

src/ringfit/data.py

@@ -66,9 +66,11 @@ class ScanData:
         """
 
         directory = Path(directory)
-        laser = np.load(directory / "signal_laser.npy")
+        laserpath = directory / "signal_laser.npy"
+
         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)
 
         return cls(laser, output, time, **kwargs)
 
@@ -119,7 +121,7 @@ class ScanData:
         """
 
         return ScanData(
-            self._laser, self._output, self._time, max_frequency=max_frequency
+            self._laser, self._output.copy(), self._time, max_frequency=max_frequency
         )
 
     @functools.cache

src/ringfit/plotting.py

@@ -32,12 +32,16 @@ def plot_scan(
     steps: bool | int = False,
     normalize=False,
     smoothe_output: bool | float = False,
+    max_frequency: float | bool = 0,
     ax=None,
     **kwargs,
 ):
     if not (laser or output):
         raise ValueError("At least one of 'laser' or 'output' must be True.")
 
+    if max_frequency:
+        data = data.sparsified(max_frequency)
+
     if smoothe_output:
         if not isinstance(smoothe_output, float):
             smoothe_output = 10 ** (-7)