diff --git a/src/qibocal/protocols/characterization/rabi/amplitude.py b/src/qibocal/protocols/characterization/rabi/amplitude.py
index 59db25485..59f673b7c 100644
--- a/src/qibocal/protocols/characterization/rabi/amplitude.py
+++ b/src/qibocal/protocols/characterization/rabi/amplitude.py
@@ -9,6 +9,7 @@
 from qibolab.qubits import QubitId
 from qibolab.sweeper import Parameter, Sweeper, SweeperType
 from scipy.optimize import curve_fit
+from scipy.signal import find_peaks
 
 from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
 from qibocal.config import log
@@ -161,9 +162,10 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:
         # Guessing period using fourier transform
         ft = np.fft.rfft(y)
         mags = abs(ft)
-        index = np.argmax(mags) if np.argmax(mags) != 0 else np.argmax(mags[1:]) + 1
-        f = x[index] / (x[1] - x[0])
-
+        local_maxima = find_peaks(mags, threshold=10)[0]
+        index = local_maxima[0] if len(local_maxima) > 0 else None
+        # 0.5 hardcoded guess for less than one oscillation
+        f = x[index] / (x[1] - x[0]) if index is not None else 0.5
         pguess = [0.5, 1, f, np.pi / 2]
         try:
             popt, _ = curve_fit(
diff --git a/src/qibocal/protocols/characterization/rabi/length.py b/src/qibocal/protocols/characterization/rabi/length.py
index 333e9ed91..895f86f7d 100644
--- a/src/qibocal/protocols/characterization/rabi/length.py
+++ b/src/qibocal/protocols/characterization/rabi/length.py
@@ -9,6 +9,7 @@
 from qibolab.qubits import QubitId
 from qibolab.sweeper import Parameter, Sweeper, SweeperType
 from scipy.optimize import curve_fit
+from scipy.signal import find_peaks
 
 from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
 from qibocal.config import log
@@ -168,8 +169,10 @@ def _fit(data: RabiLengthData) -> RabiLengthResults:
         # Guessing period using fourier transform
         ft = np.fft.rfft(y)
         mags = abs(ft)
-        index = np.argmax(mags) if np.argmax(mags) != 0 else np.argmax(mags[1:]) + 1
-        f = x[index] / (x[1] - x[0])
+        local_maxima = find_peaks(mags, threshold=1)[0]
+        index = local_maxima[0] if len(local_maxima) > 0 else None
+        # 0.5 hardcoded guess for less than one oscillation
+        f = x[index] / (x[1] - x[0]) if index is not None else 0.5
 
         pguess = [1, 1, f, np.pi / 2, 0]
         try:
diff --git a/src/qibocal/protocols/characterization/ramsey.py b/src/qibocal/protocols/characterization/ramsey.py
index b996905f3..ab01e20f3 100644
--- a/src/qibocal/protocols/characterization/ramsey.py
+++ b/src/qibocal/protocols/characterization/ramsey.py
@@ -10,6 +10,7 @@
 from qibolab.qubits import QubitId
 from qibolab.sweeper import Parameter, Sweeper, SweeperType
 from scipy.optimize import curve_fit
+from scipy.signal import find_peaks
 
 from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
 from qibocal.config import log
@@ -227,8 +228,11 @@ def _fit(data: RamseyData) -> RamseyResults:
             ft = np.fft.rfft(y)
             freqs = np.fft.rfftfreq(len(y), x[1] - x[0])
             mags = abs(ft)
-            index = np.argmax(mags) if np.argmax(mags) != 0 else np.argmax(mags[1:]) + 1
-            f = freqs[index] * 2 * np.pi
+
+            local_maxima = find_peaks(mags, threshold=10)[0]
+            index = local_maxima[0] if len(local_maxima) > 0 else None
+            # 0.5 hardcoded guess for less than one oscillation
+            f = freqs[index] * 2 * np.pi if index is not None else 0.5
             p0 = [
                 0.5,
                 0.5,