Repository: pasqal-io/Pulser

Issue:
Adapt pulses to channel's response time
Physical HW channel's have finite response times and are thus incapable of arbitrarily steep changes. It would be good to have this effect represented, such that the user knows the difference between what they want to execute and an approximation of what actually comes out of the channel.

This feature encompasses two main aspects:
- The inclusion of the response time of a channel as another of its attributes
- A method to adapt a waveform/pulse to a given channel, such that the resulting waveform/pulse is, in good approximation, what will happen on the device. This would bring the representation of a sequence closer to a "what you see is what you get" scenario.

Discussion:
review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/pasqal-io/Pulser/pull/312"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>

Patch:
diff --git a/pulser/_seq_drawer.py b/pulser/_seq_drawer.py
index c6f6a2f..d65d526 100644
--- a/pulser/_seq_drawer.py
+++ b/pulser/_seq_drawer.py
@@ -111,6 +111,8 @@ def draw_sequence(
     draw_interp_pts: bool = True,
     draw_phase_shifts: bool = False,
     draw_register: bool = False,
+    draw_input: bool = True,
+    draw_modulation: bool = False,
 ) -> tuple[Figure, Figure]:
     """Draws the entire sequence.
 
@@ -129,6 +131,11 @@ def draw_sequence(
         draw_register (bool): Whether to draw the register before the pulse
             sequence, with a visual indication (square halo) around the qubits
             masked by the SLM, defaults to False.
+        draw_input(bool): Draws the programmed pulses on the channels, defaults
+            to True.
+        draw_modulation(bool): Draws the expected channel output, defaults to
+            False. If the channel does not have a defined 'mod_bandwidth', this
+            is skipped unless 'draw_input=False'.
     """
 
     def phase_str(phi: float) -> str:
@@ -256,9 +263,22 @@ def draw_sequence(
         # Compare pulse with an interpolated pulse with 100 times more samples
         teff = np.arange(0, max(solver_time), delta_t / 100)
 
+    # Make sure the time axis of all channels are aligned
+    final_t = total_duration / time_scale
+    if draw_modulation:
+        for ch, ch_obj in seq._channels.items():
+            final_t = max(
+                final_t,
+                (seq.get_duration(ch) + 2 * ch_obj.rise_time) / time_scale,
+            )
+    t_min = -final_t * 0.03
+    t_max = final_t * 1.05
+
     for ch, (a, b) in ch_axes.items():
-        basis = seq._channels[ch].basis
-        t = np.array(data[ch]["time"]) / time_scale
+        ch_obj = seq._channels[ch]
+        basis = ch_obj.basis
+        times = np.array(data[ch]["time"])
+        t = times / time_scale
         ya = data[ch]["amp"]
         yb = data[ch]["detuning"]
         if sampling_rate:
@@ -276,8 +296,17 @@ def draw_sequence(
             yaeff = cs_amp(teff)
             ybeff = cs_detuning(teff)
 
-        t_min = -t[-1] * 0.03
-        t_max = t[-1] * 1.05
+        draw_output = draw_modulation and (
+            ch_obj.mod_bandwidth or not draw_input
+        )
+        if draw_output:
+            t_diffs = np.diff(times)
+            input_a = np.repeat(ya[1:], t_diffs)
+            input_b = np.repeat(yb[1:], t_diffs)
+            end_index = int(final_t * time_scale)
+            ya_mod = ch_obj.modulate(input_a)[:end_index]
+            yb_mod = ch_obj.modulate(input_b, keep_ends=True)[:end_index]
+
         a.set_xlim(t_min, t_max)
         b.set_xlim(t_min, t_max)
 
@@ -294,16 +323,26 @@ def draw_sequence(
         det_bottom = det_min - det_range * 0.05
         b.set_ylim(det_bottom, det_top)
 
-        a.plot(t, ya, color="darkgreen", linewidth=0.8)
-        b.plot(t, yb, color="indigo", linewidth=0.8)
+        if draw_input:
+            a.plot(t, ya, color="darkgreen", linewidth=0.8)
+            b.plot(t, yb, color="indigo", linewidth=0.8)
         if sampling_rate:
             a.plot(teff, yaeff, color="darkgreen", linewidth=0.8)
             b.plot(teff, ybeff, color="indigo", linewidth=0.8, ls="-")
             a.fill_between(teff, 0, yaeff, color="darkgreen", alpha=0.3)
             b.fill_between(teff, 0, ybeff, color="indigo", alpha=0.3)
-        else:
+        elif draw_input:
             a.fill_between(t, 0, ya, color="darkgreen", alpha=0.3)
             b.fill_between(t, 0, yb, color="indigo", alpha=0.3)
+        if draw_output:
+            a.plot(ya_mod, color="darkred", linewidth=0.8)
+            b.plot(yb_mod, color="gold", linewidth=0.8)
+            a.fill_between(
+                np.arange(ya_mod.size), 0, ya_mod, color="darkred", alpha=0.3
+            )
+            b.fill_between(
+                np.arange(yb_mod.size), 0, yb_mod, color="gold", alpha=0.3
+            )
         a.set_ylabel(r"$\Omega$ (rad/µs)", fontsize=14, labelpad=10)
         b.set_ylabel(r"$\delta$ (rad/µs)", fontsize=14)
 
@@ -367,7 +406,7 @@ def draw_sequence(
             tgt_txt_y = max_amp * 1.1 - 0.25 * (len(targets) - 1)
             tgt_str = "\n".join(tgt_strs)
             if coords == "initial":
-                x = t_min + t[-1] * 0.005
+                x = t_min + final_t * 0.005
                 target_regions.append([0, targets])
                 if seq._channels[ch].addressing == "Global":
                     a.text(
@@ -410,7 +449,7 @@ def draw_sequence(
                 a.axvspan(ti, tf, alpha=0.4, color="grey", hatch="//")
                 b.axvspan(ti, tf, alpha=0.4, color="grey", hatch="//")
                 a.text(
-                    tf + t[-1] * 5e-3,
+                    tf + final_t * 5e-3,
                     tgt_txt_y,
                     tgt_str,
                     ha="left",
@@ -420,7 +459,7 @@ def draw_sequence(
                 if phase and draw_phase_shifts:
                     msg = r"$\phi=$" + phase_str(phase)
                     wrd_len = len(max(tgt_strs, key=len))
-                    x = tf + t[-1] * 0.01 * (wrd_len + 1)
+                    x = tf + final_t * 0.01 * (wrd_len + 1)
                     a.text(
                         x,
                         max_amp * 1.1,
@@ -431,7 +470,7 @@ def draw_sequence(
                     )
         # Terminate the last open regions
         if target_regions:
-            target_regions[-1].append(t[-1])
+            target_regions[-1].append(final_t)
         for start, targets_, end in (
             target_regions if draw_phase_shifts else []
         ):
@@ -449,7 +488,7 @@ def draw_sequence(
                 b.axvline(t_, **conf)
                 msg = "\u27F2 " + phase_str(delta)
                 a.text(
-                    t_ - t[-1] * 8e-3,
+                    t_ - final_t * 8e-3,
                     max_amp * 1.1,
                     msg,
                     ha="right",
@@ -463,7 +502,7 @@ def draw_sequence(
             a.axvspan(0, tf_m, color="black", alpha=0.1, zorder=-100)
             b.axvspan(0, tf_m, color="black", alpha=0.1, zorder=-100)
             tgt_strs = [str(q) for q in seq._slm_mask_targets]
-            tgt_txt_x = t[-1] * 0.005
+            tgt_txt_x = final_t * 0.005
             tgt_txt_y = b.get_ylim()[0]
             tgt_str = "\n".join(tgt_strs)
             b.text(
@@ -478,7 +517,7 @@ def draw_sequence(
         if "measurement" in data[ch]:
             msg = f"Basis: {data[ch]['measurement']}"
             b.text(
-                t[-1] * 1.025,
+                final_t * 1.025,
                 det_top,
                 msg,
                 ha="center",
@@ -487,8 +526,8 @@ def draw_sequence(
                 color="white",
                 rotation=90,
             )
-            a.axvspan(t[-1], t_max, color="midnightblue", alpha=1)
-            b.axvspan(t[-1], t_max, color="midnightblue", alpha=1)
+            a.axvspan(final_t, t_max, color="midnightblue", alpha=1)
+            b.axvspan(final_t, t_max, color="midnightblue", alpha=1)
             a.axhline(0, xmax=0.95, linestyle="-", linewidth=0.5, color="grey")
             b.axhline(0, xmax=0.95, linestyle=":", linewidth=0.5, color="grey")
         else:
diff --git a/pulser/channels.py b/pulser/channels.py
index 694779e..89015e9 100644
--- a/pulser/channels.py
+++ b/pulser/channels.py
@@ -19,9 +19,17 @@ from dataclasses import dataclass
 from typing import cast, ClassVar, Optional
 import warnings
 
+import numpy as np
+from numpy.typing import ArrayLike
+from scipy.fft import fft, ifft, fftfreq
+
 # Warnings of adjusted waveform duration appear just once
 warnings.filterwarnings("once", "A duration of")
 
+# Conversion factor from modulation bandwith to rise time
+# For more info, see https://tinyurl.com/bdeumc8k
+MODBW_TO_TR = 0.48
+
 
 @dataclass(init=True, repr=False, frozen=True)
 class Channel:
@@ -48,6 +56,7 @@ class Channel:
             clock cycle.
         min_duration: The shortest duration an instruction can take.
         max_duration: The longest duration an instruction can take.
+        mod_bandwidth: The modulation bandwidth at -3dB (50% redution), in MHz.
 
     Example:
         To create a channel targeting the 'ground-rydberg' transition globally,
@@ -66,6 +75,19 @@ class Channel:
     clock_period: int = 4  # ns
     min_duration: int = 16  # ns
     max_duration: int = 67108864  # ns
+    mod_bandwidth: Optional[float] = None  # MHz
+
+    @property
+    def rise_time(self) -> int:
+        """The rise time (in ns).
+
+        Defined as the time taken to go from 10% to 90% output in response to
+        a step change in the input.
+        """
+        if self.mod_bandwidth:
+            return int(MODBW_TO_TR / self.mod_bandwidth * 1e3)
+        else:
+            return 0
 
     @classmethod
     def Local(
@@ -161,6 +183,88 @@ class Channel:
             )
         return _duration
 
+    def modulate(
+        self, input_samples: np.ndarray, keep_ends: bool = False
+    ) -> np.ndarray:
+        """Modulates the input according to the channel's modulation bandwidth.
+
+        Args:
+            input_samples (np.ndarray): The samples to modulate.
+            keep_ends (bool): Assume the end values of the samples were kept
+                constant (i.e. there is no ramp from zero on the ends).
+
+        Returns:
+            np.ndarray: The modulated output signal.
+        """
+        if not self.mod_bandwidth:
+            warnings.warn(
+                f"No modulation bandwidth defined for channel '{self}',"
+                " 'Channel.modulate()' returns the 'input_samples' unchanged.",
+                stacklevel=2,
+            )
+            return input_samples
+
+        # The cutoff frequency (fc) and the modulation transfer function
+        # are defined in https://tinyurl.com/bdeumc8k
+        fc = self.mod_bandwidth * 1e-3 / np.sqrt(np.log(2))
+        if keep_ends:
+            samples = np.pad(input_samples, 2 * self.rise_time, mode="edge")
+        else:
+            samples = np.pad(input_samples, self.rise_time)
+        freqs = fftfreq(samples.size)
+        modulation = np.exp(-(freqs**2) / fc**2)
+        mod_samples = ifft(fft(samples) * modulation).real
+        if keep_ends:
+            # Cut off the extra ends
+            return cast(
+                np.ndarray, mod_samples[self.rise_time : -self.rise_time]
+            )
+        return cast(np.ndarray, mod_samples)
+
+    def calc_modulation_buffer(
+        self,
+        input_samples: ArrayLike,
+        mod_samples: ArrayLike,
+        max_allowed_diff: float = 1e-2,
+    ) -> tuple[int, int]:
+        """Calculates the minimal buffers needed around a modulated waveform.
+
+        Args:
+            input_samples (ArrayLike): The input samples.
+            mod_samples (ArrayLike): The modulated samples. Must be of size
+                ``len(input_samples) + 2 * self.rise_time``.
+            max_allowed_diff (float): The maximum allowed difference between
+                the input and modulated samples at the end points.
+
+        Returns:
+            tuple[int, int]: The minimum buffer times at the start and end of
+            the samples, in ns.
+        """
+        if not self.mod_bandwidth:
+            raise TypeError(
+                f"The channel {self} doesn't have a modulation bandwidth."
+            )
+
+        tr = self.rise_time
+        samples = np.pad(input_samples, tr)
+        diffs = np.abs(samples - mod_samples) <= max_allowed_diff
+        try:
+            # Finds the last index in the start buffer that's below the max
+            # allowed diff. Considers that the waveform could start at the next
+            # indice (hence the -1, since we are subtracting from tr)
+            start = tr - np.argwhere(diffs[:tr])[-1][0] - 1
+        except IndexError:
+            start = tr
+        try:
+            # Finds the first index in the end buffer that's below the max
+            # allowed diff. The index value found matches the minimum length
+            # for this end buffer.
+            end = np.argwhere(diffs[-tr:])[0][0]
+        except IndexError:
+            end = tr
+
+        return start, end
+
     def __repr__(self) -> str:
         config = (
             f".{self.addressing}(Max Absolute Detuning: "
@@ -175,6 +279,8 @@ class Channel:
             if cast(int, self.max_targets) > 1:
                 config += f", Max targets: {self.max_targets}"
         config += f", Basis: '{self.basis}'"
+        if self.mod_bandwidth:
+            config += f", Modulation Bandwidth: {self.mod_bandwidth} MHz"
         return self.name + config + ")"
 
 
diff --git a/pulser/devices/__init__.py b/pulser/devices/__init__.py
index 40c88be..0660fdc 100644
--- a/pulser/devices/__init__.py
+++ b/pulser/devices/__init__.py
@@ -13,12 +13,12 @@
 # limitations under the License.
 """Valid devices for Pulser Sequence execution."""
 
-from pulser.devices._devices import (
-    Chadoq2,
-)
+from pulser.devices._devices import Chadoq2
 
 from pulser.devices._mock_device import MockDevice
 
+from pulser.devices._device_datacls import Device
+
 # Registers which devices can be used to avoid definition of custom devices
 _mock_devices = (MockDevice,)
 _valid_devices = (Chadoq2,)
diff --git a/pulser/pulse.py b/pulser/pulse.py
index db28249..c3a2a75 100644
--- a/pulser/pulse.py
+++ b/pulser/pulse.py
@@ -23,6 +23,7 @@ from typing import Any, cast, Union
 import matplotlib.pyplot as plt
 import numpy as np
 
+from pulser.channels import Channel
 from pulser.parametrized import Parametrized, ParamObj
 from pulser.parametrized.decorators import parametrize
 from pulser.waveforms import Waveform, ConstantWaveform
@@ -187,6 +188,15 @@ class Pulse:
         fig.tight_layout()
         plt.show()
 
+    def fall_time(self, channel: Channel) -> int:
+        """Calculates the extra time needed to ramp down to zero."""
+        aligned_start_extra_time = channel.rise_time
+        end_extra_time = max(
+            self.amplitude.modulation_buffers(channel)[1],
+            self.detuning.modulation_buffers(channel)[1],
+        )
+        return aligned_start_extra_time + end_extra_time
+
     def _to_dict(self) -> dict[str, Any]:
         return obj_to_dict(
             self,
diff --git a/pulser/sequence.py b/pulser/sequence.py
index 78d0944..cdd85f1 100644
--- a/pulser/sequence.py
+++ b/pulser/sequence.py
@@ -282,27 +282,50 @@ class Sequence:
         return not self._building
 
     @_screen
-    def get_duration(self, channel: Optional[str] = None) -> int:
+    def get_duration(
+        self, channel: Optional[str] = None, include_fall_time: bool = False
+    ) -> int:
         """Returns the current duration of a channel or the whole sequence.
 
         Keyword Args:
             channel (Optional[str]): A specific channel to return the duration
                 of. If left as None, it will return the duration of the whole
                 sequence.
+            include_fall_time (bool): Whether to include in the duration the
+                extra time needed by the last pulse to finish, if there is
+                modulation.
 
         Returns:
             int: The duration of the channel or sequence, in ns.
         """
         if channel is None:
-            durations = [
-                self._last(ch).tf
-                for ch in self._schedule
-                if self._schedule[ch]
-            ]
-            return 0 if not durations else max(durations)
+            channels = tuple(self._channels.keys())
+            if not channels:
+                return 0
+        else:
+            self._validate_channel(channel)
+            channels = (channel,)
+        last_ts = {}
+        for id in channels:
+            this_chobj = self._channels[id]
+            temp_tf = 0
+            for i, op in enumerate(self._schedule[id][::-1]):
+                if i == 0:
+                    # Start with the last slot found
+                    temp_tf = op.tf
+                    if not include_fall_time:
+                        break
+                if isinstance(op.type, Pulse):
+                    temp_tf = max(
+                        temp_tf, op.tf + op.type.fall_time(this_chobj)
+                    )
+                    break
+                elif temp_tf - op.tf >= 2 * this_chobj.rise_time:
+                    # No pulse behind 'op' with a long enough fall time
+                    break
+            last_ts[id] = temp_tf
 
-        self._validate_channel(channel)
-        return self._last(channel).tf if self._schedule[channel] else 0
+        return max(last_ts.values())
 
     @_screen
     def current_phase_ref(
@@ -652,13 +675,20 @@ class Sequence:
             for ch, seq in self._schedule.items():
                 if ch == channel:
                     continue
+                this_chobj = self._channels[ch]
                 for op in self._schedule[ch][::-1]:
-                    if op.tf <= current_max_t:
-                        break
                     if not isinstance(op.type, Pulse):
-                        continue
-                    if op.targets & last.targets or protocol == "wait-for-all":
-                        current_max_t = op.tf
+                        if op.tf + 2 * this_chobj.rise_time <= current_max_t:
+                            # No pulse behind 'op' needing a delay
+                            break
+                    elif (
+                        op.tf + op.type.fall_time(this_chobj) <= current_max_t
+                    ):
+                        break
+                    elif (
+                        op.targets & last.targets or protocol == "wait-for-all"
+                    ):
+                        current_max_t = op.tf + op.type.fall_time(this_chobj)
                         break
 
         delay_duration = current_max_t - t0
@@ -674,8 +704,13 @@ class Sequence:
                 break
 
         if delay_duration > 0:
-            # Delay must not be shorter than the min duration for this channel
-            delay_duration = max(delay_duration, channel_obj.min_duration)
+            # Delay must not be shorter than the min duration of this channel
+            # and a multiple of the clock period (forced by validate_duration)
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                delay_duration = channel_obj.validate_duration(
+                    max(delay_duration, channel_obj.min_duration)
+                )
             self._delay(delay_duration, channel)
 
         ti = t0 + delay_duration
@@ -683,9 +718,10 @@ class Sequence:
 
         self._add_to_schedule(channel, _TimeSlot(pulse, ti, tf, last.targets))
 
+        true_finish = tf + pulse.fall_time(channel_obj)
         for qubit in last.targets:
-            if self._last_used[basis][qubit] < tf:
-                self._last_used[basis][qubit] = tf
+            if self._last_used[basis][qubit] < true_finish:
+                self._last_used[basis][qubit] = true_finish
 
         if pulse.post_phase_shift:
             self._phase_shift(
@@ -834,13 +870,23 @@ class Sequence:
         if self.is_parametrized():
             return
 
-        last_ts = {id: self._last(cast(str, id)).tf for id in channels}
+        channels = cast(Tuple[str], channels)
+        last_ts = {
+            id: self.get_duration(id, include_fall_time=True)
+            for id in channels
+        }
         tf = max(last_ts.values())
 
         for id in channels:
             delta = tf - last_ts[id]
             if delta > 0:
-                self._delay(delta, cast(str, id))
+                channel_obj = self._channels[id]
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore")
+                    delta = channel_obj.validate_duration(
+                        max(delta, channel_obj.min_duration)
+                    )
+                self._delay(delta, id)
 
     def build(self, **vars: Union[ArrayLike, float, int, str]) -> Sequence:
         """Builds a sequence from the programmed instructions.
@@ -957,6 +1003,7 @@ class Sequence:
     @_screen
     def draw(
         self,
+        mode: str = "input+output",
         draw_phase_area: bool = False,
         draw_interp_pts: bool = True,
         draw_phase_shifts: bool = False,
@@ -967,11 +1014,18 @@ class Sequence:
         """Draws the sequence in its current state.
 
         Keyword Args:
+            mode (str, default="input+output"): The curves to draw. 'input'
+                draws only the programmed curves, 'output' the excepted curves
+                after modulation. 'input+output' will draw both curves except
+                for channels without a defined modulation bandwidth, in which
+                case only the input is drawn.
             draw_phase_area (bool): Whether phase and area values need to be
-                shown as text on the plot, defaults to False.
+                shown as text on the plot, defaults to False. Doesn't work in
+                'output' mode.
             draw_interp_pts (bool): When the sequence has pulses with waveforms
                 of type InterpolatedWaveform, draws the points of interpolation
-                on top of the respective waveforms (defaults to True).
+                on top of the respective input waveforms (defaults to True).
+                Doesn't work in 'output' mode.
             draw_phase_shifts (bool): Whether phase shift and reference
                 information should be added to the plot, defaults to False.
             draw_register (bool): Whether to draw the register before the pulse
@@ -991,12 +1045,34 @@ class Sequence:
             Simulation.draw(): Draws the provided sequence and the one used by
             the solver.
         """
+        valid_modes = ("input", "output", "input+output")
+        if mode not in valid_modes:
+            raise ValueError(
+                f"'mode' must be one of {valid_modes}, not '{mode}'."
+            )
+        if mode == "output":
+            if draw_phase_area:
+                warnings.warn(
+                    "'draw_phase_area' doesn't work in 'output' mode, so it "
+                    "will default to 'False'.",
+                    stacklevel=2,
+                )
+                draw_phase_area = False
+            if draw_interp_pts:
+                warnings.warn(
+                    "'draw_interp_pts' doesn't work in 'output' mode, so it "
+                    "will default to 'False'.",
+                    stacklevel=2,
+                )
+                draw_interp_pts = False
         fig_reg, fig = draw_sequence(
             self,
             draw_phase_area=draw_phase_area,
             draw_interp_pts=draw_interp_pts,
             draw_phase_shifts=draw_phase_shifts,
             draw_register=draw_register,
+            draw_input="input" in mode,
+            draw_modulation="output" in mode,
         )
         if fig_name is not None and draw_register:
             name, ext = os.path.splitext(fig_name)
@@ -1010,7 +1086,7 @@ class Sequence:
         self, qubits: Union[Iterable[QubitId], QubitId], channel: str
     ) -> None:
         self._validate_channel(channel)
-
+        channel_obj = self._channels[channel]
         try:
             qubits_set = (
                 set(cast(Iterable, qubits))
@@ -1020,12 +1096,12 @@ class Sequence:
         except TypeError:
             qubits_set = {qubits}
 
-        if self._channels[channel].addressing != "Local":
+        if channel_obj.addressing != "Local":
             raise ValueError("Can only choose target of 'Local' channels.")
-        elif len(qubits_set) > cast(int, self._channels[channel].max_targets):
+        elif len(qubits_set) > cast(int, channel_obj.max_targets):
             raise ValueError(
-                "This channel can target at most "
-                f"{self._channels[channel].max_targets} qubits at a time"
+                f"This channel can target at most {channel_obj.max_targets} "
+                "qubits at a time."
             )
 
         if self.is_parametrized():
@@ -1039,7 +1115,7 @@ class Sequence:
         elif not qubits_set.issubset(self._qids):
             raise ValueError("All given qubits must belong to the register.")
 
-        basis = self._channels[channel].basis
+        basis = channel_obj.basis
         phase_refs = {self._phase_ref[basis][q].last_phase for q in qubits_set}
         if len(phase_refs) != 1:
             raise ValueError(
@@ -1048,22 +1124,33 @@ class Sequence:
             )
 
         try:
+            fall_time = self.get_duration(
+                channel, include_fall_time=True
+            ) - self.get_duration(channel)
+            if fall_time > 0:
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore")
+                    self.delay(
+                        max(fall_time, channel_obj.min_duration),
+                        channel,
+                    )
+
             last = self._last(channel)
             if last.targets == qubits_set:
                 return
             ti = last.tf
-            retarget = cast(int, self._channels[channel].min_retarget_interval)
+            retarget = cast(int, channel_obj.min_retarget_interval)
             elapsed = ti - self._last_target[channel]
             delta = cast(int, np.clip(retarget - elapsed, 0, retarget))
+            if channel_obj.fixed_retarget_t:
+                delta = max(delta, channel_obj.fixed_retarget_t)
             if delta != 0:
                 with warnings.catch_warnings():
                     warnings.simplefilter("ignore")
-                    delta = self._channels[channel].validate_duration(
-                        16 if delta < 16 else delta
+                    delta = channel_obj.validate_duration(
+                        max(delta, channel_obj.min_duration)
                     )
-            tf = ti + max(
-                delta, cast(int, self._channels[channel].fixed_retarget_t)
-            )
+            tf = ti + delta
 
         except ValueError:
             ti = -1
diff --git a/pulser/waveforms.py b/pulser/waveforms.py
index b21ce77..7daddbf 100644
--- a/pulser/waveforms.py
+++ b/pulser/waveforms.py
@@ -31,6 +31,7 @@ import numpy as np
 from numpy.typing import ArrayLike
 import scipy.interpolate as interpolate
 
+from pulser.channels import Channel
 from pulser.parametrized import Parametrized, ParamObj
 from pulser.parametrized.decorators import parametrize
 from pulser.json.utils import obj_to_dict
@@ -123,11 +124,28 @@ class Waveform(ABC):
         """Integral of the waveform (time in ns, value in rad/µs)."""
         return float(np.sum(self.samples)) * 1e-3  # ns * rad/µs = 1e-3
 
-    def draw(self) -> None:
-        """Draws the waveform."""
-        fig, ax = plt.subplots()
-        self._plot(ax, "rad/µs")
+    def draw(self, output_channel: Optional[Channel] = None) -> None:
+        """Draws the waveform.
 
+        Args:
+            output_channel: The output channel. If given, will draw the
+                modulated waveform on top of the input one.
+        """
+        fig, ax = plt.subplots()
+        if not output_channel:
+            self._plot(ax, "rad/µs")
+        else:
+            self._plot(
+                ax,
+                "rad/µs",
+                label="Input",
+                start_t=self.modulation_buffers(output_channel)[0],
+            )
+            self._plot(
+                ax,
+                channel=output_channel,
+                label="Output",
+            )
         plt.show()
 
     def change_duration(self, new_duration: int) -> Waveform:
@@ -141,6 +159,56 @@ class Waveform(ABC):
             " modifications to its duration."
         )
 
+    def modulated_samples(self, channel: Channel) -> np.ndarray:
+        """The waveform samples as output of a given channel.
+
+        This duration is adjusted according to the minimal buffer times.
+
+        Args:
+            channel (Channel): The channel modulating the waveform.
+
+        Returns:
+            numpy.ndarray: The array of samples after modulation.
+        """
+        start, end = self.modulation_buffers(channel)
+        mod_samples = self._modulated_samples(channel)
+        tr = channel.rise_time
+        trim = slice(tr - start, len(mod_samples) - tr + end)
+        return cast(np.ndarray, mod_samples[trim])
+
+    @functools.lru_cache()
+    def modulation_buffers(self, channel: Channel) -> tuple[int, int]:
+        """The minimal buffers needed around a modulated waveform.
+
+        Args:
+            channel (Channel): The channel modulating the waveform.
+
+        Returns:
+            tuple[int, int]: The minimum buffer times at the start and end of
+            the samples, in ns.
+        """
+        if not channel.mod_bandwidth:
+            return 0, 0
+
+        return channel.calc_modulation_buffer(
+            self._samples, self._modulated_samples(channel)
+        )
+
+    @functools.lru_cache()
+    def _modulated_samples(self, channel: Channel) -> np.ndarray:
+        """The waveform samples as output of a given channel.
+
+        This is not adjusted to the minimal buffer times. Use
+        ``Waveform.modulated_samples()`` to get the output already truncated.
+
+        Args:
+            channel (Channel): The channel modulating the waveform.
+
+        Returns:
+            numpy.ndarray: The array of samples after modulation.
+        """
+        return channel.modulate(self._samples)
+
     @abstractmethod
     def _to_dict(self) -> dict[str, Any]:
         pass
@@ -229,19 +297,42 @@ class Waveform(ABC):
         return hash(tuple(self.samples))
 
     def _plot(
-        self, ax: Axes, ylabel: str, color: Optional[str] = None
+        self,
+        ax: Axes,
+        ylabel: Optional[str] = None,
+        color: Optional[str] = None,
+        channel: Optional[Channel] = None,
+        label: str = "",
+        start_t: int = 0,
     ) -> None:
         ax.set_xlabel("t (ns)")
-        ts = np.arange(self.duration)
+        samples = (
+            self.samples
+            if channel is None
+            else self.modulated_samples(channel)
+        )
+        ts = np.arange(len(samples)) + start_t
+        if not channel and start_t:
+            # Adds zero on both ends to show rise and fall
+            samples = np.pad(samples, 1)
+            # Repeats the times on the edges once
+            ts = np.pad(ts, 1, mode="edge")
+
         if color:
-            ax.set_ylabel(ylabel, color=color, fontsize=14)
-            ax.plot(ts, self.samples, color=color)
+            color_dict = {"color": color}
+            hline_color = color
             ax.tick_params(axis="y", labelcolor=color)
-            ax.axhline(0, color=color, linestyle=":", linewidth=0.5)
         else:
-            ax.set_ylabel(ylabel, fontsize=14)
-            ax.plot(ts, self.samples)
-            ax.axhline(0, color="black", linestyle=":", linewidth=0.5)
+            color_dict = {}
+            hline_color = "black"
+
+        if ylabel:
+            ax.set_ylabel(ylabel, fontsize=14, **color_dict)
+        ax.plot(ts, samples, label=label, **color_dict)
+        ax.axhline(0, color=hline_color, linestyle=":", linewidth=0.5)
+
+        if label:
+            plt.legend()
 
 
 class CompositeWaveform(Waveform):
@@ -723,10 +814,26 @@ class InterpolatedWaveform(Waveform):
         return InterpolatedWaveform(new_duration, self._values, **self._kwargs)
 
     def _plot(
-        self, ax: Axes, ylabel: str, color: Optional[str] = None
+        self,
+        ax: Axes,
+        ylabel: Optional[str] = None,
+        color: Optional[str] = None,
+        channel: Optional[Channel] = None,
+        label: str = "",
+        start_t: int = 0,
     ) -> None:
-        super()._plot(ax, ylabel, color=color)
-        ax.scatter(self._data_pts[:, 0], self._data_pts[:, 1], c=color)
+        super()._plot(
+            ax,
+            ylabel,
+            color=color,
+            channel=channel,
+            label=label,
+            start_t=start_t,
+        )
+        if not channel:
+            ax.scatter(
+                self._data_pts[:, 0] + start_t, self._data_pts[:, 1], c=color
+            )
 
     def _to_dict(self) -> dict[str, Any]:
         return obj_to_dict(self, self._duration, self._values, **self._kwargs)
diff --git a/tutorials/applications/Using QAOA to solve a MIS problem.ipynb b/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
index cbc101d..d2f72c2 100644
--- a/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
+++ b/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
@@ -61,7 +61,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "For example, assume an ensemble of identical radio transmitters over French cities that each have the same radius of transmission. It was quickly realized that two transmitters with close or equal frequencies could interfere with one another, hence the necessity to assign non-interfering frequencies to overlapping transmiting towers. Because of the limited amount of bandwith space, some towers have to be assigned the same or close frequencies. The MIS of a graph of towers indicate the maximum number of towers that can have close or equal given frequency (red points). \n",
+    "For example, assume an ensemble of identical radio transmitters over French cities that each have the same radius of transmission. It was quickly realized that two transmitters with close or equal frequencies could interfere with one another, hence the necessity to assign non-interfering frequencies to overlapping transmiting towers. Because of the limited amount of bandwidth space, some towers have to be assigned the same or close frequencies. The MIS of a graph of towers indicate the maximum number of towers that can have close or equal given frequency (red points). \n",
     "\n",
     "<center>\n",
     "<img src=\"attachment:France_MIS.png\" alt=\"MIS map France\" width=\"650\">\n",


Test patch:
diff --git a/pulser/tests/test_channels.py b/pulser/tests/test_channels.py
index ce56cf1..4c25202 100644
--- a/pulser/tests/test_channels.py
+++ b/pulser/tests/test_channels.py
@@ -12,10 +12,12 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+import numpy as np
 import pytest
 
 import pulser
 from pulser.channels import Raman, Rydberg
+from pulser.waveforms import ConstantWaveform, BlackmanWaveform
 
 
 def test_device_channels():
@@ -64,10 +66,43 @@ def test_repr():
     )
     assert raman.__str__() == r1
 
-    ryd = Rydberg.Global(50, 2.5, phase_jump_time=300)
+    ryd = Rydberg.Global(50, 2.5, phase_jump_time=300, mod_bandwidth=4)
     r2 = (
         "Rydberg.Global(Max Absolute Detuning: 50 rad/µs, "
         "Max Amplitude: 2.5 rad/µs, Phase Jump Time: 300 ns, "
-        "Basis: 'ground-rydberg')"
+        "Basis: 'ground-rydberg', Modulation Bandwidth: 4 MHz)"
     )
     assert ryd.__str__() == r2
+
+
+def test_modulation():
+    rydberg_global = Rydberg.Global(2 * np.pi * 20, 2 * np.pi * 2.5)
+
+    raman_local = Raman.Local(
+        2 * np.pi * 20,
+        2 * np.pi * 10,
+        mod_bandwidth=4,  # MHz
+    )
+
+    wf = ConstantWaveform(100, 1)
+    assert rydberg_global.mod_bandwidth is None
+    with pytest.warns(UserWarning, match="No modulation bandwidth defined"):
+        out_samples = rydberg_global.modulate(wf.samples)
+    assert np.all(out_samples == wf.samples)
+
+    with pytest.raises(TypeError, match="doesn't have a modulation bandwidth"):
+        rydberg_global.calc_modulation_buffer(wf.samples, out_samples)
+
+    out_ = raman_local.modulate(wf.samples)
+    tr = raman_local.rise_time
+    assert len(out_) == wf.duration + 2 * tr
+    assert raman_local.calc_modulation_buffer(wf.samples, out_) == (tr, tr)
+
+    wf2 = BlackmanWaveform(800, np.pi)
+    side_buffer_len = 45
+    out_ = raman_local.modulate(wf2.samples)
+    assert len(out_) == wf2.duration + 2 * tr  # modulate() does not truncate
+    assert raman_local.calc_modulation_buffer(wf2.samples, out_) == (
+        side_buffer_len,
+        side_buffer_len,
+    )
diff --git a/pulser/tests/test_sequence.py b/pulser/tests/test_sequence.py
index a3d0bfa..1899b1e 100644
--- a/pulser/tests/test_sequence.py
+++ b/pulser/tests/test_sequence.py
@@ -19,6 +19,7 @@ import pytest
 
 import pulser
 from pulser import Sequence, Pulse, Register, Register3D
+from pulser.channels import Rydberg, Raman
 from pulser.devices import Chadoq2, MockDevice
 from pulser.devices._device_datacls import Device
 from pulser.sequence import _TimeSlot
@@ -554,3 +555,106 @@ def test_draw_register():
     seq3d.config_slm_mask([6, 15])
     with patch("matplotlib.pyplot.show"):
         seq3d.draw(draw_register=True)
+
+
+def test_hardware_constraints():
+    rydberg_global = Rydberg.Global(
+        2 * np.pi * 20,
+        2 * np.pi * 2.5,
+        phase_jump_time=120,  # ns
+        mod_bandwidth=4,  # MHz
+    )
+
+    raman_local = Raman.Local(
+        2 * np.pi * 20,
+        2 * np.pi * 10,
+        phase_jump_time=120,  # ns
+        fixed_retarget_t=200,  # ns
+        mod_bandwidth=7,  # MHz
+    )
+
+    ConstrainedChadoq2 = Device(
+        name="ConstrainedChadoq2",
+        dimensions=2,
+        rydberg_level=70,
+        max_atom_num=100,
+        max_radial_distance=50,
+        min_atom_distance=4,
+        _channels=(
+            ("rydberg_global", rydberg_global),
+            ("raman_local", raman_local),
+        ),
+    )
+    with pytest.warns(
+        UserWarning, match="should be imported from 'pulser.devices'"
+    ):
+        seq = Sequence(reg, ConstrainedChadoq2)
+    seq.declare_channel("ch0", "rydberg_global")
+    seq.declare_channel("ch1", "raman_local", initial_target="q1")
+
+    const_pls = Pulse.ConstantPulse(100, 1, 0, np.pi)
+    seq.add(const_pls, "ch0")
+    black_wf = BlackmanWaveform(500, np.pi)
+    black_pls = Pulse.ConstantDetuning(black_wf, 0, 0)
+    seq.add(black_pls, "ch1")
+    blackman_slot = seq._last("ch1")
+    # The pulse accounts for the modulation buffer
+    assert (
+        blackman_slot.ti == const_pls.duration + rydberg_global.rise_time * 2
+    )
+    seq.target("q0", "ch1")
+    target_slot = seq._last("ch1")
+    fall_time = black_pls.fall_time(raman_local)
+    assert (
+        fall_time
+        == raman_local.rise_time + black_wf.modulation_buffers(raman_local)[1]
+    )
+    fall_time += (
+        raman_local.clock_period - fall_time % raman_local.clock_period
+    )
+    assert target_slot.ti == blackman_slot.tf + fall_time
+    assert target_slot.tf == target_slot.ti + raman_local.fixed_retarget_t
+
+    assert raman_local.min_retarget_interval > raman_local.fixed_retarget_t
+    seq.target("q2", "ch1")
+    assert (
+        seq.get_duration("ch1")
+        == target_slot.tf + raman_local.min_retarget_interval
+    )
+
+    # Check for phase jump buffer
+    seq.add(black_pls, "ch0")  # Phase = 0
+    tf_ = seq.get_duration("ch0")
+    mid_delay = 40
+    seq.delay(mid_delay, "ch0")
+    seq.add(const_pls, "ch0")  # Phase = π
+    assert seq._last("ch0").ti - tf_ == rydberg_global.phase_jump_time
+    added_delay_slot = seq._schedule["ch0"][-2]
+    assert added_delay_slot.type == "delay"
+    assert (
+        added_delay_slot.tf - added_delay_slot.ti
+        == rydberg_global.phase_jump_time - mid_delay
+    )
+
+    tf_ = seq.get_duration("ch0")
+    seq.align("ch0", "ch1")
+    fall_time = const_pls.fall_time(rydberg_global)
+    assert seq.get_duration() == tf_ + fall_time
+
+    with pytest.raises(ValueError, match="'mode' must be one of"):
+        seq.draw(mode="all")
+
+    with patch("matplotlib.pyplot.show"):
+        with pytest.warns(
+            UserWarning,
+            match="'draw_phase_area' doesn't work in 'output' mode",
+        ):
+            seq.draw(
+                mode="output", draw_interp_pts=False, draw_phase_area=True
+            )
+        with pytest.warns(
+            UserWarning,
+            match="'draw_interp_pts' doesn't work in 'output' mode",
+        ):
+            seq.draw(mode="output")
+        seq.draw(mode="input+output")
diff --git a/pulser/tests/test_waveforms.py b/pulser/tests/test_waveforms.py
index c048736..fa10d51 100644
--- a/pulser/tests/test_waveforms.py
+++ b/pulser/tests/test_waveforms.py
@@ -21,6 +21,7 @@ import pytest
 from scipy.interpolate import interp1d, PchipInterpolator
 
 from pulser.json.coders import PulserEncoder, PulserDecoder
+from pulser.channels import Rydberg
 from pulser.parametrized import Variable, ParamObj
 from pulser.waveforms import (
     ConstantWaveform,
@@ -98,10 +99,16 @@ def test_integral():
 
 
 def test_draw():
+    rydberg_global = Rydberg.Global(
+        2 * np.pi * 20,
+        2 * np.pi * 2.5,
+        phase_jump_time=120,  # ns
+        mod_bandwidth=4,  # MHz
+    )
     with patch("matplotlib.pyplot.show"):
         composite.draw()
-        blackman.draw()
-        interp.draw()
+        blackman.draw(output_channel=rydberg_global)
+        interp.draw(output_channel=rydberg_global)
 
 
 def test_eq():
@@ -405,3 +412,21 @@ def test_get_item():
         assert wf[duration * 2 :].size == 0
         assert wf[duration * 2 : duration * 3].size == 0
         assert wf[-duration * 3 : -duration * 2].size == 0
+
+
+def test_modulation():
+    rydberg_global = Rydberg.Global(
+        2 * np.pi * 20,
+        2 * np.pi * 2.5,
+        phase_jump_time=120,  # ns
+        mod_bandwidth=4,  # MHz
+    )
+    mod_samples = constant.modulated_samples(rydberg_global)
+    assert np.all(mod_samples == rydberg_global.modulate(constant.samples))
+    assert constant.modulation_buffers(rydberg_global) == (
+        rydberg_global.rise_time,
+        rydberg_global.rise_time,
+    )
+    assert len(mod_samples) == constant.duration + 2 * rydberg_global.rise_time
+    assert np.isclose(np.sum(mod_samples) * 1e-3, constant.integral)
+    assert max(np.abs(mod_samples)) < np.abs(constant[0])