[v1.5.x] Add sound()

examples/scripts/write_gre.py

@@ -5,7 +5,13 @@
 import pypulseq as pp
 
 
-def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_pypulseq.seq', paper_plot: bool = False):
+def main(
+    plot: bool = False,
+    write_seq: bool = False,
+    seq_filename: str = 'gre_pypulseq.seq',
+    paper_plot: bool = False,
+    play_sound: bool = False,
+):
     # ======
     # SETUP
     # ======
@@ -143,8 +149,11 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_p
 
         seq.write(seq_filename)
 
+    if play_sound:
+        seq.sound()
+
     return seq
 
 
 if __name__ == '__main__':
-    seq = main(plot=False, paper_plot=False, write_seq=False)
+    seq = main(plot=False, paper_plot=False, write_seq=False, play_sound=True)

pyproject.toml

@@ -32,6 +32,7 @@ dependencies = [
 [project.optional-dependencies]
 sigpy = ["sigpy>=0.1.26"]
 mplcursors = ["mplcursors"]
+sound = ["sounddevice"]
 Qt = ["PySide6"]
 test = [
   "coverage",

src/pypulseq/Sequence/sequence.py

@@ -33,6 +33,7 @@
 from pypulseq.utils.cumsum import cumsum
 from pypulseq.utils.paper_plot import paper_plot as ext_paper_plot
 from pypulseq.utils.seq_plot import SeqPlot
+from pypulseq.utils.seq_sound import seq_sound
 from pypulseq.utils.tracing import format_trace, trace, trace_enabled
 
 major, minor, revision = __version__.split('.')[:3]
@@ -863,6 +864,36 @@ def get_gradients(
             gw_pp.append(PPoly(np.stack((np.diff(gw[1]) / np.diff(gw[0]), gw[1][:-1])), gw[0], extrapolate=True))
         return gw_pp
 
+    def sound(
+        self,
+        time_range=(0.0, np.inf),
+        channel_weights: tuple = (1.0, 1.0, 1.0),
+        sample_rate: int = 44100,
+        only_produce_sound_data: bool = False,
+    ) -> np.ndarray:
+        """
+        Play out the sequence through the system speaker and return the sound data.
+
+        Parameters
+        ----------
+        time_range : iterable, default=(0, np.inf)
+            Time range for playing the sequence.
+            Default is 0 to infinity (entire sequence).
+        channel_weights : tuple, default=(1.0, 1.0, 1.0)
+            Weight for each gradient channel.
+        sample_rate : int, default=44100
+            Sampling rate for audio waveform
+        only_produce_sound_data : bool, default=False
+            If True, skip the playing part and only produce
+            the sound data
+
+        Returns
+        -------
+        sound_data : np.ndarray
+            Sound data of shape (2, N)
+        """
+        return seq_sound(self, time_range, channel_weights, sample_rate, only_produce_sound_data)
+
     def install(self, target: Union[str, None] = None, clear_cache: bool = False, **kwargs: Any) -> None:
         """Install a sequence to a target scanner.
 

src/pypulseq/utils/seq_sound.py

@@ -0,0 +1,122 @@
+from __future__ import annotations
+
+import typing
+
+import numpy as np
+from scipy.interpolate import interp1d
+from scipy.signal import convolve, windows
+
+try:
+    import sounddevice as sd
+
+    _HAS_SOUNDDEVICE = True
+except ImportError:
+    _HAS_SOUNDDEVICE = False
+
+if typing.TYPE_CHECKING:
+    from pypulseq.Sequence.sequence import Sequence
+
+
+def seq_sound(
+    seq: Sequence,
+    time_range=(0.0, np.inf),
+    channel_weights: tuple = (1.0, 1.0, 1.0),
+    sample_rate: int = 44100,
+    only_produce_sound_data: bool = False,
+) -> np.ndarray:
+    """
+    Play out the sequence through the system speaker and return the sound data.
+
+    Parameters
+    ----------
+    seq : Sequence
+        The Pulseq sequence object to played.
+    time_range : iterable, default=(0, np.inf)
+        Time range for playing the sequence.
+        Default is 0 to infinity (entire sequence).
+    channel_weights : tuple, default=(1.0, 1.0, 1.0)
+        Weight for each gradient channel.
+    sample_rate : int, default=44100
+        Sampling rate for audio waveform
+    only_produce_sound_data : bool, default=False
+        If True, skip the playing part and only produce
+        the sound data
+
+    Returns
+    -------
+    sound_data : np.ndarray
+        Sound data of shape (2, N)
+    """
+    # --- get waveform data ---
+    gw_data = seq.waveforms_and_times(False, time_range)
+    total_duration = np.sum(seq.block_durations)
+
+    dwell_time = 1.0 / sample_rate
+    sound_length = int(np.floor(total_duration / dwell_time)) + 1
+
+    sound_data = np.zeros((2, sound_length), dtype=float)
+
+    t_out = np.arange(sound_length) * dwell_time
+
+    # --- X channel → channel 1 ---
+    if gw_data[0] is not None and gw_data[0].size > 0:
+        f = interp1d(
+            gw_data[0][0],
+            gw_data[0][1] * channel_weights[0],
+            kind='linear',
+            bounds_error=False,
+            fill_value=0.0,
+        )
+        sound_data[0, :] = f(t_out)
+
+    # --- Y channel → channel 2 ---
+    if gw_data[1] is not None and gw_data[1].size > 0:
+        f = interp1d(
+            gw_data[1][0],
+            gw_data[1][1] * channel_weights[1],
+            kind='linear',
+            bounds_error=False,
+            fill_value=0.0,
+        )
+        sound_data[1, :] = f(t_out)
+
+    # --- Z channel → split equally ---
+    if gw_data[2] is not None and gw_data[2].size > 0:
+        f = interp1d(
+            gw_data[2][0],
+            0.5 * gw_data[2][1] * channel_weights[2],
+            kind='linear',
+            bounds_error=False,
+            fill_value=0.0,
+        )
+        tmp = f(t_out)
+        sound_data[0, :] += tmp
+        sound_data[1, :] += tmp
+
+    # --- Gaussian smoothing (same as MATLAB gausswin + conv) ---
+    gw_len = int(round(sample_rate / 6000) * 2 + 1)
+    gw = windows.gaussian(gw_len, std=gw_len / 6)
+    gw /= np.sum(gw)
+
+    sound_data[0, :] = convolve(sound_data[0, :], gw, mode='same')
+    sound_data[1, :] = convolve(sound_data[1, :], gw, mode='same')
+
+    # --- normalize ---
+    max_val = np.max(np.abs(sound_data))
+    if max_val > 0:
+        sound_data = 0.95 * sound_data / max_val
+
+    # --- play sound ---
+    if not only_produce_sound_data:
+        duration = sound_length * dwell_time
+        print(f'playing out the sequence waveform, duration {duration:.1f}s')
+
+        if not _HAS_SOUNDDEVICE:
+            raise RuntimeError('sounddevice not installed')
+
+        pad = np.zeros((2, sample_rate // 2))
+        play_data = np.hstack([pad, sound_data, pad]).T  # (N, 2)
+        sd.play(play_data, samplerate=sample_rate)
+        sd.wait()
+
+    return sound_data