From bd9f8e17048a6d1ddc8d9193aca253293180af7e Mon Sep 17 00:00:00 2001 From: Valentin Boettcher Date: Tue, 16 Jul 2024 09:51:18 -0400 Subject: [PATCH] some experimenatal data analysis --- .../007_generate_fake_analysis_data.py | 2 +- .../008_analyze_11_07_transients.py | 90 +++++++++++++++++++ .../009_analyze_12_07_transients.py | 88 ++++++++++++++++++ .../001_first_exploration.py | 2 +- src/rabifun/analysis.py | 15 +++- src/ringfit/data.py | 23 +++-- 6 files changed, 207 insertions(+), 13 deletions(-) create mode 100644 scripts/experiments/008_analyze_11_07_transients.py create mode 100644 scripts/experiments/009_analyze_12_07_transients.py diff --git a/scripts/experiments/007_generate_fake_analysis_data.py b/scripts/experiments/007_generate_fake_analysis_data.py index ffdf91c..0f189be 100644 --- a/scripts/experiments/007_generate_fake_analysis_data.py +++ b/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, diff --git a/scripts/experiments/008_analyze_11_07_transients.py b/scripts/experiments/008_analyze_11_07_transients.py new file mode 100644 index 0000000..f4f4ce7 --- /dev/null +++ b/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, + ) diff --git a/scripts/experiments/009_analyze_12_07_transients.py b/scripts/experiments/009_analyze_12_07_transients.py new file mode 100644 index 0000000..d14afc2 --- /dev/null +++ b/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, + ) diff --git a/scripts/ringdown_spectrum_analysis/001_first_exploration.py b/scripts/ringdown_spectrum_analysis/001_first_exploration.py index a023e00..460275b 100644 --- a/scripts/ringdown_spectrum_analysis/001_first_exploration.py +++ b/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 diff --git a/src/rabifun/analysis.py b/src/rabifun/analysis.py index ccc8e12..d407e0f 100644 --- a/src/rabifun/analysis.py +++ b/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] diff --git a/src/ringfit/data.py b/src/ringfit/data.py index 06ad371..bacddce 100644 --- a/src/ringfit/data.py +++ b/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):