Add alternative Ne/Si calibration method based on argmin2d

src/ramanchada2/protocols/calib_ne_si_argmin2d_iter_gg.py

@@ -0,0 +1,44 @@
+from typing import Literal
+
+import numpy as np
+from scipy import interpolate
+
+from ..misc import constants as rc2const
+from ..spectrum.spectrum import Spectrum
+
+
+def neon_calibration(ne_cm_1: Spectrum,
+                     wl: Literal[514, 532, 633, 785]):
+    ref = rc2const.neon_wl_dict[wl]
+    ne_nm = ne_cm_1.subtract_moving_minimum(200).shift_cm_1_to_abs_nm_filter(wl).normalize()
+
+    ne_cal = ne_nm.xcal_argmin2d_iter_lowpass(ref=ref)
+    spline = interpolate.Akima1DInterpolator(ne_cm_1.x, ne_cal.x, method='makima')
+    return spline
+
+
+def silicon_calibration(si_nm: Spectrum,
+                        wl: Literal[514, 532, 633, 785]):
+    si_nm = si_nm.dropna().normalize()
+    peaks = si_nm.find_peak_multipeak(prominence=.05, width=2, wlen=50)
+    fitres = si_nm.fit_peak_multimodel(candidates=peaks, profile='Pearson4')
+    si_wl = fitres.centers
+    laser_wl = 1/(520.45/1e7 + 1/si_wl)
+
+    laser_wl = laser_wl[np.argmin(np.abs(laser_wl-wl))]
+    x_cm_1 = 1e7*(1/laser_wl-1/si_nm.x)
+
+    spline = interpolate.Akima1DInterpolator(si_nm.x, x_cm_1, method='makima')
+    return spline, laser_wl
+
+
+def neon_silicon_calibration(ne_cm_1: Spectrum,
+                             si_cm_1: Spectrum,
+                             wl: Literal[514, 532, 633, 785]):
+    ne_spline = neon_calibration(ne_cm_1, wl)
+    si_nm = si_cm_1.scale_xaxis_fun(ne_spline)
+    si_spline, wl = silicon_calibration(si_nm, wl)
+    ne_nm = ne_cm_1.scale_xaxis_fun(ne_spline)
+    ne_cal_cm_1 = ne_nm.abs_nm_to_shift_cm_1_filter(wl)
+    spline = interpolate.Akima1DInterpolator(ne_cm_1.x, ne_cal_cm_1.x, method='makima')
+    return spline

src/ramanchada2/spectrum/calibration/by_deltas.py

@@ -5,6 +5,7 @@
 import numpy.typing as npt
 from pydantic import BaseModel, NonNegativeInt, validate_call
 from scipy import interpolate
+from scipy import fft, signal
 
 from ramanchada2.misc.spectrum_deco import (add_spectrum_filter,
                                             add_spectrum_method)
@@ -251,3 +252,43 @@ def xcal_fine_RBF(old_spe: Spectrum,
         kwargs["kernel"] = kernel
         interp = interpolate.RBFInterpolator(spe_cent[spe_idx].reshape(-1, 1), ref_pos[ref_idx], **kwargs)
         new_spe.x = interp(old_spe.x.reshape(-1, 1))
+
+
+@add_spectrum_filter
+@validate_call(config=dict(arbitrary_types_allowed=True))
+def xcal_argmin2d_iter_lowpass(old_spe: Spectrum,
+                               new_spe: Spectrum, /, *,
+                               ref: Dict[float, float],
+                               low_pass_nfreqs: List[int] = [100, 500]):
+    def semi_spe_from_dict(deltas: dict, xaxis):
+        y = np.zeros_like(xaxis)
+        for pos, ampl in deltas.items():
+            idx = np.argmin(np.abs(xaxis - pos))
+            y[idx] += ampl
+        # remove overflows and underflows
+        y[0] = 0
+        y[-1] = 0
+        return y
+
+    def low_pass(x, nbin, window=signal.windows.blackmanharris):
+        h = window(nbin*2-1)[nbin-1:]
+        X = fft.rfft(x)
+        X[:nbin] *= h  # apply the window
+        X[nbin:] = 0  # clear upper frequencies
+        return fft.irfft(X, n=len(x))
+
+    spe = old_spe.__copy__()
+    for low_pass_i in low_pass_nfreqs:
+        xaxis = spe.x
+        y_ref_semi_spe = semi_spe_from_dict(ref, spe.x)
+        y_ref_semi_spe = low_pass(y_ref_semi_spe, low_pass_i)
+
+        r = xaxis[signal.find_peaks(y_ref_semi_spe)[0]]
+
+        spe_low = spe.__copy__()
+        spe_low.y = low_pass(spe.y, low_pass_i)
+
+        spe_cal = spe_low.xcal_fine(ref=r, should_fit=False, poly_order=2)
+        spe.x = spe_cal.x
+    spe_cal_fin = spe.xcal_fine(ref=ref, should_fit=False, poly_order=2)
+    new_spe.x = spe_cal_fin.x