Add sampling rate to the waveform tools (#185)

* add sampling_rate arg for OPX1k

* Fix DRAG + update tests

* changelog

* sampling with 1Gs/s resolution for flattop pulses

CHANGELOG.md

@@ -4,6 +4,8 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
 ## [Unreleased]
+### Changed
+- config/waveform_tools - Added sampling rate argument with default value set to 1GS/s to the waveforms.
 
 ## [0.16.0] - 2024-01-25
 ### Fixed

qualang_tools/config/waveform_tools.py

@@ -3,7 +3,7 @@
 
 
 def drag_gaussian_pulse_waveforms(
-    amplitude, length, sigma, alpha, anharmonicity, detuning=0.0, subtracted=True, **kwargs
+    amplitude, length, sigma, alpha, anharmonicity, detuning=0.0, subtracted=True, sampling_rate=1e9, **kwargs
 ):
     """
     Creates Gaussian based DRAG waveforms that compensate for the leakage and for the AC stark shift.
@@ -22,16 +22,17 @@ def drag_gaussian_pulse_waveforms(
     :param float detuning: The frequency shift to correct for AC stark shift, in Hz.
     :param bool subtracted: If true, returns a subtracted Gaussian, such that the first and last points will be at 0
         volts. This reduces high-frequency components due to the initial and final points offset. Default is true.
-    :return: Returns a tuple of two lists. The first list is the I waveform (real part) and the second is the
-        Q waveform (imaginary part)
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Default is 1G samples/s.
+    :return: Returns a tuple of two lists. The first list is the 'I' waveform (real part) and the second is the
+        'Q' waveform (imaginary part)
     """
     delta = kwargs.get("delta", None)
     if delta is not None:
         print("'delta' has been replaced by 'anharmonicity' and will be deprecated in the future. ")
         if alpha != 0 and delta == 0:
             raise Exception("Cannot create a DRAG pulse with `anharmonicity=0`")
-        t = np.arange(length, dtype=int)  # An array of size pulse length in ns
-        center = (length - 1) / 2
+        t = np.arange(length, step=1e9 / sampling_rate)  # An array of size pulse length in ns
+        center = (length - 1e9 / sampling_rate) / 2
         gauss_wave = amplitude * np.exp(-((t - center) ** 2) / (2 * sigma**2))  # The gaussian function
         gauss_der_wave = (
             amplitude * (-2 * 1e9 * (t - center) / (2 * sigma**2)) * np.exp(-((t - center) ** 2) / (2 * sigma**2))
@@ -49,8 +50,8 @@ def drag_gaussian_pulse_waveforms(
     else:
         if alpha != 0 and anharmonicity == 0:
             raise Exception("Cannot create a DRAG pulse with `anharmonicity=0`")
-        t = np.arange(length, dtype=int)  # An array of size pulse length in ns
-        center = (length - 1) / 2
+        t = np.arange(length, step=1e9 / sampling_rate)  # An array of size pulse length in ns
+        center = (length - 1e9 / sampling_rate) / 2
         gauss_wave = amplitude * np.exp(-((t - center) ** 2) / (2 * sigma**2))  # The gaussian function
         gauss_der_wave = (
             amplitude * (-2 * 1e9 * (t - center) / (2 * sigma**2)) * np.exp(-((t - center) ** 2) / (2 * sigma**2))
@@ -68,7 +69,7 @@ def drag_gaussian_pulse_waveforms(
     return I_wf, Q_wf
 
 
-def drag_cosine_pulse_waveforms(amplitude, length, alpha, anharmonicity, detuning=0.0, **kwargs):
+def drag_cosine_pulse_waveforms(amplitude, length, alpha, anharmonicity, detuning=0.0, sampling_rate=1e9, **kwargs):
     """
     Creates Cosine based DRAG waveforms that compensate for the leakage and for the AC stark shift.
 
@@ -83,16 +84,17 @@ def drag_cosine_pulse_waveforms(amplitude, length, alpha, anharmonicity, detunin
     :param float alpha: The DRAG coefficient.
     :param float anharmonicity: f_21 - f_10 - The differences in energy between the 2-1 and the 1-0 energy levels, in Hz.
     :param float detuning: The frequency shift to correct for AC stark shift, in Hz.
-    :return: Returns a tuple of two lists. The first list is the I waveform (real part) and the second is the
-        Q waveform (imaginary part)
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Default is 1G samples/s.
+    :return: Returns a tuple of two lists. The first list is the 'I' waveform (real part) and the second is the
+        'Q' waveform (imaginary part)
     """
     delta = kwargs.get("delta", None)
     if delta is not None:
         print("'delta' has been replaced by 'anharmonicity' and will be deprecated in the future.")
         if alpha != 0 and anharmonicity == 0:
             raise Exception("Cannot create a DRAG pulse with `anharmonicity=0`")
-        t = np.arange(length, dtype=int)  # An array of size pulse length in ns
-        end_point = length - 1
+        t = np.arange(length, step=1e9 / sampling_rate)  # An array of size pulse length in ns
+        end_point = length - 1e9 / sampling_rate
         cos_wave = 0.5 * amplitude * (1 - np.cos(t * 2 * np.pi / end_point))  # The cosine function
         sin_wave = (
             0.5 * amplitude * (2 * np.pi / end_point * 1e9) * np.sin(t * 2 * np.pi / end_point)
@@ -108,8 +110,8 @@ def drag_cosine_pulse_waveforms(amplitude, length, alpha, anharmonicity, detunin
     else:
         if alpha != 0 and anharmonicity == 0:
             raise Exception("Cannot create a DRAG pulse with `anharmonicity=0`")
-        t = np.arange(length, dtype=int)  # An array of size pulse length in ns
-        end_point = length - 1
+        t = np.arange(length, step=1e9 / sampling_rate)  # An array of size pulse length in ns
+        end_point = length - 1e9 / sampling_rate
         cos_wave = 0.5 * amplitude * (1 - np.cos(t * 2 * np.pi / end_point))  # The cosine function
         sin_wave = (
             0.5 * amplitude * (2 * np.pi / end_point * 1e9) * np.sin(t * 2 * np.pi / end_point)
@@ -125,7 +127,7 @@ def drag_cosine_pulse_waveforms(amplitude, length, alpha, anharmonicity, detunin
     return I_wf, Q_wf
 
 
-def flattop_gaussian_waveform(amplitude, flat_length, rise_fall_length, return_part="all"):
+def flattop_gaussian_waveform(amplitude, flat_length, rise_fall_length, return_part="all", sampling_rate=1e9):
     """
     Returns a flat top Gaussian waveform. This is a square pulse with a rise and fall of a Gaussian with the given
     sigma. It is possible to only get the rising or falling parts, which allows scanning the flat part length from QUA.
@@ -138,14 +140,19 @@ def flattop_gaussian_waveform(amplitude, flat_length, rise_fall_length, return_p
     :param str return_part: When set to 'all', returns the complete waveform. Default is 'all'. When set to 'rise',
     returns only the rising part. When set to 'fall', returns only the falling part. This is useful for separating
     the three parts which allows scanning the duration of the flat part is to scanned from QUA
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Must be an integer multiple of 1e9 samples per seconds. Default is 1G samples/s.
     :return: Returns the waveform as a list of values with 1ns spacing
     """
-    gauss_wave = amplitude * gaussian(2 * rise_fall_length, rise_fall_length / 5)
 
-    rise_part = gauss_wave[:rise_fall_length]
+    assert sampling_rate % 1e9 == 0, "The sampling rate must be an integer multiple of 1e9 samples per second."
+
+    gauss_wave = amplitude * gaussian(
+        int(np.round(2 * rise_fall_length * sampling_rate / 1e9)), rise_fall_length / 5 * sampling_rate / 1e9
+    )
+    rise_part = gauss_wave[: int(rise_fall_length * sampling_rate / 1e9)]
     rise_part = rise_part.tolist()
     if return_part == "all":
-        return rise_part + [amplitude] * flat_length + rise_part[::-1]
+        return rise_part + [amplitude] * int(flat_length * sampling_rate / 1e9) + rise_part[::-1]
     elif return_part == "rise":
         return rise_part
     elif return_part == "fall":
@@ -154,7 +161,7 @@ def flattop_gaussian_waveform(amplitude, flat_length, rise_fall_length, return_p
         raise Exception("'return_part' must be either 'all', 'rise' or 'fall'")
 
 
-def flattop_cosine_waveform(amplitude, flat_length, rise_fall_length, return_part="all"):
+def flattop_cosine_waveform(amplitude, flat_length, rise_fall_length, return_part="all", sampling_rate=1e9):
     """
     Returns a flat top cosine waveform. This is a square pulse with a rise and fall with cosine shape with the given
     sigma. It is possible to only get the rising or falling parts, which allows scanning the flat part length from QUA.
@@ -167,12 +174,14 @@ def flattop_cosine_waveform(amplitude, flat_length, rise_fall_length, return_par
     :param str return_part: When set to 'all', returns the complete waveform. Default is 'all'. When set to 'rise',
     returns only the rising part. When set to 'fall', returns only the falling part. This is useful for separating
     the three parts which allows scanning the duration of the flat part is to scanned from QUA
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Must be an integer multiple of 1e9 samples per seconds. Default is 1G samples/s.
     :return: Returns the waveform as a list of values with 1ns spacing
     """
-    rise_part = amplitude * 0.5 * (1 - np.cos(np.linspace(0, np.pi, rise_fall_length)))
+    assert sampling_rate % 1e9 == 0, "The sampling rate must be an integer multiple of 1e9 samples per second."
+    rise_part = amplitude * 0.5 * (1 - np.cos(np.linspace(0, np.pi, int(rise_fall_length * sampling_rate / 1e9))))
     rise_part = rise_part.tolist()
     if return_part == "all":
-        return rise_part + [amplitude] * flat_length + rise_part[::-1]
+        return rise_part + [amplitude] * int(flat_length * sampling_rate / 1e9) + rise_part[::-1]
     elif return_part == "rise":
         return rise_part
     elif return_part == "fall":
@@ -181,7 +190,7 @@ def flattop_cosine_waveform(amplitude, flat_length, rise_fall_length, return_par
         raise Exception("'return_part' must be either 'all', 'rise' or 'fall'")
 
 
-def flattop_tanh_waveform(amplitude, flat_length, rise_fall_length, return_part="all"):
+def flattop_tanh_waveform(amplitude, flat_length, rise_fall_length, return_part="all", sampling_rate=1e9):
     """
     Returns a flat top tanh waveform. This is a square pulse with a rise and fall with tanh shape with the given
     sigma. It is possible to only get the rising or falling parts, which allows scanning the flat part length from QUA.
@@ -194,12 +203,14 @@ def flattop_tanh_waveform(amplitude, flat_length, rise_fall_length, return_part=
     :param str return_part: When set to 'all', returns the complete waveform. Default is 'all'. When set to 'rise',
     returns only the rising part. When set to 'fall', returns only the falling part. This is useful for separating
     the three parts which allows scanning the duration of the flat part is to scanned from QUA
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Must be an integer multiple of 1e9 samples per seconds. Default is 1G samples/s.
     :return: Returns the waveform as a list of values with 1ns spacing
     """
-    rise_part = amplitude * 0.5 * (1 + np.tanh(np.linspace(-4, 4, rise_fall_length)))
+    assert sampling_rate % 1e9 == 0, "The sampling rate must be an integer multiple of 1e9 samples per second."
+    rise_part = amplitude * 0.5 * (1 + np.tanh(np.linspace(-4, 4, int(rise_fall_length * sampling_rate / 1e9))))
     rise_part = rise_part.tolist()
     if return_part == "all":
-        return rise_part + [amplitude] * flat_length + rise_part[::-1]
+        return rise_part + [amplitude] * int(flat_length * sampling_rate / 1e9) + rise_part[::-1]
     elif return_part == "rise":
         return rise_part
     elif return_part == "fall":
@@ -208,7 +219,7 @@ def flattop_tanh_waveform(amplitude, flat_length, rise_fall_length, return_part=
         raise Exception("'return_part' must be either 'all', 'rise' or 'fall'")
 
 
-def flattop_blackman_waveform(amplitude, flat_length, rise_fall_length, return_part="all"):
+def flattop_blackman_waveform(amplitude, flat_length, rise_fall_length, return_part="all", sampling_rate=1e9):
     """
     Returns a flat top Blackman waveform. This is a square pulse with a rise and fall with Blackman shape with the given
     length. It is possible to only get the rising or falling parts, which allows scanning the flat part length from QUA.
@@ -220,13 +231,15 @@ def flattop_blackman_waveform(amplitude, flat_length, rise_fall_length, return_p
     :param str return_part: When set to 'all', returns the complete waveform. Default is 'all'. When set to 'rise',
     returns only the rising part. When set to 'fall', returns only the falling part. This is useful for separating
     the three parts which allows scanning the duration of the  flat part is to scanned from QUA
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Must be an integer multiple of 1e9 samples per seconds. Default is 1G samples/s.
     :return: Returns the waveform as a list
     """
-    backman_wave = amplitude * blackman(2 * rise_fall_length)
-    rise_part = backman_wave[:rise_fall_length]
+    assert sampling_rate % 1e9 == 0, "The sampling rate must be an integer multiple of 1e9 samples per second."
+    backman_wave = amplitude * blackman(2 * int(rise_fall_length * sampling_rate / 1e9))
+    rise_part = backman_wave[: int(rise_fall_length * sampling_rate / 1e9)]
     rise_part = rise_part.tolist()
     if return_part == "all":
-        return rise_part + [amplitude] * flat_length + rise_part[::-1]
+        return rise_part + [amplitude] * int(flat_length * sampling_rate / 1e9) + rise_part[::-1]
     elif return_part == "rise":
         return rise_part
     elif return_part == "fall":
@@ -235,17 +248,19 @@ def flattop_blackman_waveform(amplitude, flat_length, rise_fall_length, return_p
         raise Exception("'return_part' must be either 'all', 'rise' or 'fall'")
 
 
-def blackman_integral_waveform(pulse_length, v_start, v_end):
+def blackman_integral_waveform(pulse_length, v_start, v_end, sampling_rate=1e9):
     """
     Returns a Blackman integral waveform. This is the integral of a Blackman waveform, adiabatically going from
     'v_start' to 'v_end' in 'pulse_length' ns.
 
     :param int pulse_length: The pulse length in ns.
     :param float v_start: The starting amplitude in volts.
     :param float v_end: The ending amplitude in volts.
+    :param float sampling_rate: The sampling rate used to describe the waveform, in samples/s. Must be an integer multiple of 1e9 samples per seconds. Default is 1G samples/s.
     :return: Returns the waveform as a list
     """
-    time = np.asarray([x * 1.0 for x in range(int(pulse_length))])
+    assert sampling_rate % 1e9 == 0, "The sampling rate must be an integer multiple of 1e9 samples per second."
+    time = np.linspace(0, pulse_length - 1, int(pulse_length * sampling_rate / 1e9))
     black_wave = v_start + (
         time / (pulse_length - 1)
         - (25 / (42 * np.pi)) * np.sin(2 * np.pi * time / (pulse_length - 1))