From fad099610b477d1f1f8bf4d939f732b27fe15329 Mon Sep 17 00:00:00 2001
From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com>
Date: Mon, 8 Jan 2024 14:48:28 +0100
Subject: [PATCH 1/5] Added gradient spectrum calculations
 (`Sequence.calculate_gradient_spectrum`) - Added `asc_to_acoustic_resonances`
 function to get acoustic resonance ranges from Siemens asc files

---
 pypulseq/Sequence/calc_grad_spectrum.py | 171 ++++++++++++++++++++++++
 pypulseq/Sequence/sequence.py           |  66 ++++++++-
 pypulseq/utils/siemens/asc_to_hw.py     |  25 ++++
 3 files changed, 261 insertions(+), 1 deletion(-)
 create mode 100644 pypulseq/Sequence/calc_grad_spectrum.py

diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
new file mode 100644
index 0000000..4aa095a
--- /dev/null
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -0,0 +1,171 @@
+from typing import Tuple, List
+
+import numpy as np
+from scipy.signal import spectrogram
+from matplotlib import pyplot as plt
+
+
+def calculate_gradient_spectrum(
+        obj, max_frequency: float = 2000,
+        window_width: float = 0.05,
+        frequency_oversampling: float = 3,
+        time_range: List[float] = None,
+        plot: bool = True,
+        combine_mode: str = 'max',
+        acoustic_resonances: List[dict] = []
+) -> Tuple[List[np.ndarray], np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Calculates the gradient spectrum of the sequence. Returns a spectrogram
+    for each gradient channel, as well as a root-sum-squares combined
+    spectrogram.
+    
+    Works by splitting the sequence into windows that are 'window_width'
+    long and calculating the fourier transform of each window. Windows
+    overlap 50% with the previous and next window. When 'combine_mode' is
+    not 'none', all windows are combined into one spectrogram.
+
+    Parameters
+    ----------
+    max_frequency : float, optional
+        Maximum frequency to include in spectrograms. The default is 2000.
+    window_width : float, optional
+        Window width (in seconds). The default is 0.05.
+    frequency_oversampling : float, optional
+        Oversampling in the frequency dimension, higher values make
+        smoother spectrograms. The default is 3.
+    time_range : List[float], optional
+        Time range over which to calculate the spectrograms as a list of
+        two timepoints (in seconds) (e.g. [1, 1.5])
+        The default is None.
+    plot : bool, optional
+        Whether to plot the spectograms. The default is True.
+    combine_mode : str, optional
+        How to combine all windows into one spectrogram, options:
+            'max', 'mean', 'sos' (root-sum-of-squares), 'none' (no combination)
+        The default is 'max'.
+    acoustic_resonances : List[dict], optional
+        Acoustic resonances as a list of dictionaries with 'frequency' and
+        'bandwidth' elements. Only used when plot==True. The default is [].
+
+    Returns
+    -------
+    spectrograms : List[np.ndarray]
+        List of spectrograms per gradient channel.
+    spectrogram_sos : np.ndarray
+        Root-sum-of-squares combined spectrogram over all gradient channels.
+    frequencies : np.ndarray
+        Frequency axis of the spectrograms.
+    times : np.ndarray
+        Time axis of the spectrograms (only relevant when combine_mode == 'none').
+
+    """
+    dt = obj.system.grad_raster_time # time raster
+    nwin = round(window_width / dt)
+    nfft = round(frequency_oversampling*nwin)
+
+    # Get gradients as piecewise-polynomials
+    gw_pp = obj.get_gradients(time_range=time_range)
+    ng = len(gw_pp)
+    max_t = max(g.x[-1] for g in gw_pp if g != None)
+    
+    # Determine sampling points
+    if time_range == None:
+        nt = int(np.ceil(max_t/dt))
+        t = (np.arange(nt) + 0.5)*dt
+    else:
+        tmax = min(time_range[1], max_t) - max(time_range[0], 0)
+        nt = int(np.ceil(tmax/dt))
+        t = max(time_range[0], 0) + (np.arange(nt) + 0.5)*dt
+    
+    # Sample gradients
+    gw = np.zeros((ng,t.shape[0]))
+    for i in range(ng):
+        if gw_pp[i] != None:
+            gw[i] = gw_pp[i](t)
+    
+    # Calculate spectrogram for each gradient channel
+    spectrograms = []
+    frequencies = None
+    spectrogram_sos = 0
+    
+    for i in range(ng):
+        # Use scipy to calculate the spectrograms
+        freq, times, Sxx = spectrogram(gw[i],
+                                       fs=1/dt,
+                                       mode='magnitude',
+                                       nperseg=nwin,
+                                       noverlap=nwin//2,
+                                       nfft=nfft,
+                                       detrend='constant',
+                                       window=('tukey', 1))
+        mask = freq<max_frequency
+        
+        # Accumulate spectrum for all gradient channels
+        spectrogram_sos += Sxx[mask]**2
+        
+        # Combine spectrogram over time axis
+        if combine_mode == 'max':
+            s = Sxx[mask].max(axis=1)
+        elif combine_mode == 'mean':
+            s = Sxx[mask].mean(axis=1)
+        elif combine_mode == 'sos':
+            s = np.sqrt((Sxx[mask]**2).sum(axis=1))
+        elif combine_mode == 'none':
+            s = Sxx[mask]
+        else:
+            raise ValueError(f'Unknown value for combine_mode: {combine_mode}, must be one of [max, mean, sos, none]')
+        
+        frequencies = freq[mask]
+        spectrograms.append(s)
+    
+    # Root-sum-of-squares combined spectrogram for all gradient channels
+    spectrogram_sos = np.sqrt(spectrogram_sos)
+    if combine_mode == 'max':
+        spectrogram_sos = spectrogram_sos.max(axis=1)
+    elif combine_mode == 'mean':
+        spectrogram_sos = spectrogram_sos.mean(axis=1)
+    elif combine_mode == 'sos':
+        spectrogram_sos = np.sqrt((spectrogram_sos**2).sum(axis=1))
+    
+    # Plot spectrograms and acoustic resonances if specified
+    if plot:
+        if combine_mode != 'none':
+            plt.figure()
+            plt.xlabel('Frequency (Hz)')
+            # According to spectrogram documentation y unit is (Hz/m)^2 / Hz = Hz/m^2, is this meaningful?
+            for s in spectrograms:
+                plt.plot(frequencies, s)
+            plt.plot(frequencies, spectrogram_sos)
+            plt.legend(['x', 'y', 'z', 'sos'])
+    
+            for res in acoustic_resonances:
+                plt.axvline(res['frequency'], color='k', linestyle='-')
+                plt.axvline(res['frequency'] - res['bandwidth']/2, color='k', linestyle='--')
+                plt.axvline(res['frequency'] + res['bandwidth']/2, color='k', linestyle='--')
+        else:
+            for s, c in zip(spectrograms, ['X', 'Y', 'Z']):
+                plt.figure()
+                plt.title(f'Spectrum {c}')
+                plt.xlabel('Time (s)')
+                plt.ylabel('Frequency (Hz)')
+                plt.imshow(abs(s[::-1]), extent=(times[0], times[-1], frequencies[0], frequencies[-1]),
+                           aspect=(times[-1]-times[0])/(frequencies[-1]-frequencies[0]))
+                
+                for res in acoustic_resonances:
+                    plt.axhline(res['frequency'], color='r', linestyle='-')
+                    plt.axhline(res['frequency'] - res['bandwidth']/2, color='r', linestyle='--')
+                    plt.axhline(res['frequency'] + res['bandwidth']/2, color='r', linestyle='--')
+            
+            plt.figure()
+            plt.title('Total spectrum')
+            plt.xlabel('Time (s)')
+            plt.ylabel('Frequency (Hz)')
+            plt.imshow(abs(spectrogram_sos[::-1]), extent=(times[0], times[-1], frequencies[0], frequencies[-1]),
+                       aspect=(times[-1]-times[0])/(frequencies[-1]-frequencies[0]))
+            
+            for res in acoustic_resonances:
+                plt.axhline(res['frequency'], color='r', linestyle='-')
+                plt.axhline(res['frequency'] - res['bandwidth']/2, color='r', linestyle='--')
+                plt.axhline(res['frequency'] + res['bandwidth']/2, color='r', linestyle='--')
+                
+    return spectrograms, spectrogram_sos, frequencies, times
diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py
index 258e736..90a80c6 100755
--- a/pypulseq/Sequence/sequence.py
+++ b/pypulseq/Sequence/sequence.py
@@ -17,6 +17,8 @@
 from pypulseq.Sequence.read_seq import read
 from pypulseq.Sequence.write_seq import write as write_seq
 from pypulseq.Sequence.calc_pns import calc_pns
+from pypulseq.Sequence.calc_grad_spectrum import calculate_gradient_spectrum
+
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.check_timing import check_timing as ext_check_timing
 from pypulseq.decompress_shape import decompress_shape
@@ -190,7 +192,69 @@ def add_block(self, *args: SimpleNamespace) -> None:
         """
         block.set_block(self, self.next_free_block_ID, *args)
         self.next_free_block_ID += 1
-            
+    
+    def calculate_gradient_spectrum(
+            self, max_frequency: float = 2000,
+            window_width: float = 0.05,
+            frequency_oversampling: float = 3,
+            time_range: List[float] = None,
+            plot: bool = True,
+            combine_mode: str = 'max',
+            acoustic_resonances: List[dict] = []
+    ) -> Tuple[List[np.ndarray], np.ndarray, np.ndarray, np.ndarray]:
+        """
+        Calculates the gradient spectrum of the sequence. Returns a spectrogram
+        for each gradient channel, as well as a root-sum-squares combined
+        spectrogram.
+        
+        Works by splitting the sequence into windows that are 'window_width'
+        long and calculating the fourier transform of each window. Windows
+        overlap 50% with the previous and next window. When 'combine_mode' is
+        not 'none', all windows are combined into one spectrogram.
+
+        Parameters
+        ----------
+        max_frequency : float, optional
+            Maximum frequency to include in spectrograms. The default is 2000.
+        window_width : float, optional
+            Window width (in seconds). The default is 0.05.
+        frequency_oversampling : float, optional
+            Oversampling in the frequency dimension, higher values make
+            smoother spectrograms. The default is 3.
+        time_range : List[float], optional
+            Time range over which to calculate the spectrograms as a list of
+            two timepoints (in seconds) (e.g. [1, 1.5])
+            The default is None.
+        plot : bool, optional
+            Whether to plot the spectograms. The default is True.
+        combine_mode : str, optional
+            How to combine all windows into one spectrogram, options:
+                'max', 'mean', 'sos' (root-sum-of-squares), 'none' (no combination)
+            The default is 'max'.
+        acoustic_resonances : List[dict], optional
+            Acoustic resonances as a list of dictionaries with 'frequency' and
+            'bandwidth' elements. Only used when plot==True. The default is [].
+
+        Returns
+        -------
+        spectrograms : List[np.ndarray]
+            List of spectrograms per gradient channel.
+        spectrogram_sos : np.ndarray
+            Root-sum-of-squares combined spectrogram over all gradient channels.
+        frequencies : np.ndarray
+            Frequency axis of the spectrograms.
+        times : np.ndarray
+            Time axis of the spectrograms (only relevant when combine_mode == 'none').
+
+        """
+        return calculate_gradient_spectrum(self, max_frequency=max_frequency,
+                                           window_width=window_width,
+                                           frequency_oversampling=frequency_oversampling,
+                                           time_range=time_range,
+                                           plot=plot,
+                                           combine_mode=combine_mode,
+                                           acoustic_resonances=acoustic_resonances)
+    
     def calculate_kspace(
         self,
         trajectory_delay: Union[float, List[float], np.ndarray] = 0,
diff --git a/pypulseq/utils/siemens/asc_to_hw.py b/pypulseq/utils/siemens/asc_to_hw.py
index 72568b6..a38a065 100644
--- a/pypulseq/utils/siemens/asc_to_hw.py
+++ b/pypulseq/utils/siemens/asc_to_hw.py
@@ -1,7 +1,32 @@
 from types import SimpleNamespace
+from typing import List
 import numpy as np
 
 
+def asc_to_acoustic_resonances(asc : dict) -> List[dict]:
+    """
+    Convert ASC dictionary from readasc to list of acoustic resonances
+
+    Parameters
+    ----------
+    asc : dict
+        ASC dictionary, see readasc
+
+    Returns
+    -------
+    List[dict]
+        List of acoustic resonances (specified by frequency and bandwidth fields).
+    """
+    
+    if 'aflGCAcousticResonanceFrequency' in asc:
+        freqs = asc['aflGCAcousticResonanceFrequency']
+        bw = asc['aflGCAcousticResonanceBandwidth']
+    else:
+        freqs = asc['asGPAParameters'][0]['sGCParameters']['aflAcousticResonanceFrequency']
+        bw = asc['asGPAParameters'][0]['sGCParameters']['aflAcousticResonanceBandwidth']
+    
+    return [dict(frequency=f, bandwidth=b) for f,b in zip(freqs.values(), bw.values()) if f != 0]
+
 def asc_to_hw(asc : dict, cardiac_model : bool = False) -> SimpleNamespace:
     """
     Convert ASC dictionary from readasc to SAFE hardware description.

From 57df1390a4650a09f095bd0e65d997542a6bf8e5 Mon Sep 17 00:00:00 2001
From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com>
Date: Thu, 15 Feb 2024 17:11:38 +0100
Subject: [PATCH 2/5] Added `use_derivative` option to
 `calculate_gradient_spectrum`

---
 pypulseq/Sequence/calc_grad_spectrum.py | 13 +++++++++++--
 pypulseq/Sequence/sequence.py           |  6 ++++++
 2 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
index 4aa095a..0b0bf69 100644
--- a/pypulseq/Sequence/calc_grad_spectrum.py
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -6,13 +6,15 @@
 
 
 def calculate_gradient_spectrum(
-        obj, max_frequency: float = 2000,
+        obj,
+        max_frequency: float = 2000,
         window_width: float = 0.05,
         frequency_oversampling: float = 3,
         time_range: List[float] = None,
         plot: bool = True,
         combine_mode: str = 'max',
-        acoustic_resonances: List[dict] = []
+        use_derivative: bool = False,
+        acoustic_resonances: List[dict] = [],
 ) -> Tuple[List[np.ndarray], np.ndarray, np.ndarray, np.ndarray]:
     """
     Calculates the gradient spectrum of the sequence. Returns a spectrogram
@@ -43,6 +45,10 @@ def calculate_gradient_spectrum(
         How to combine all windows into one spectrogram, options:
             'max', 'mean', 'sos' (root-sum-of-squares), 'none' (no combination)
         The default is 'max'.
+    use_derivative : bool, optional
+        Whether the use the derivative of the gradient waveforms instead of the
+        gradient waveforms for the gradient spectrum calculations. The default
+        is False
     acoustic_resonances : List[dict], optional
         Acoustic resonances as a list of dictionaries with 'frequency' and
         'bandwidth' elements. Only used when plot==True. The default is [].
@@ -83,6 +89,9 @@ def calculate_gradient_spectrum(
         if gw_pp[i] != None:
             gw[i] = gw_pp[i](t)
     
+    if use_derivative:
+        gw = np.diff(gw, axis=1)
+    
     # Calculate spectrogram for each gradient channel
     spectrograms = []
     frequencies = None
diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py
index 90a80c6..d7afbb5 100755
--- a/pypulseq/Sequence/sequence.py
+++ b/pypulseq/Sequence/sequence.py
@@ -200,6 +200,7 @@ def calculate_gradient_spectrum(
             time_range: List[float] = None,
             plot: bool = True,
             combine_mode: str = 'max',
+            use_derivative: bool = False,
             acoustic_resonances: List[dict] = []
     ) -> Tuple[List[np.ndarray], np.ndarray, np.ndarray, np.ndarray]:
         """
@@ -231,6 +232,10 @@ def calculate_gradient_spectrum(
             How to combine all windows into one spectrogram, options:
                 'max', 'mean', 'sos' (root-sum-of-squares), 'none' (no combination)
             The default is 'max'.
+        use_derivative : bool, optional
+            Whether the use the derivative of the gradient waveforms instead of the
+            gradient waveforms for the gradient spectrum calculations. The default
+            is False
         acoustic_resonances : List[dict], optional
             Acoustic resonances as a list of dictionaries with 'frequency' and
             'bandwidth' elements. Only used when plot==True. The default is [].
@@ -253,6 +258,7 @@ def calculate_gradient_spectrum(
                                            time_range=time_range,
                                            plot=plot,
                                            combine_mode=combine_mode,
+                                           use_derivative=use_derivative,
                                            acoustic_resonances=acoustic_resonances)
     
     def calculate_kspace(

From c4efecc55fe11501c22e5a19586c931f566f3115 Mon Sep 17 00:00:00 2001
From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com>
Date: Fri, 15 Mar 2024 16:07:29 +0100
Subject: [PATCH 3/5] Requested changes in `calculate_gradient_spectrum` -
 Minor beauty fixes and type hints - Renamed 'sos' to 'rss' - Renamed
 `spectrogram_sos` to `spectrogram_rss`

---
 pypulseq/Sequence/calc_grad_spectrum.py | 47 ++++++++++++-------------
 1 file changed, 23 insertions(+), 24 deletions(-)

diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
index 0b0bf69..66d7cad 100644
--- a/pypulseq/Sequence/calc_grad_spectrum.py
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -10,7 +10,7 @@ def calculate_gradient_spectrum(
         max_frequency: float = 2000,
         window_width: float = 0.05,
         frequency_oversampling: float = 3,
-        time_range: List[float] = None,
+        time_range: List[float] | None = None,
         plot: bool = True,
         combine_mode: str = 'max',
         use_derivative: bool = False,
@@ -43,7 +43,7 @@ def calculate_gradient_spectrum(
         Whether to plot the spectograms. The default is True.
     combine_mode : str, optional
         How to combine all windows into one spectrogram, options:
-            'max', 'mean', 'sos' (root-sum-of-squares), 'none' (no combination)
+            'max', 'mean', 'rss' (root-sum-of-squares), 'none' (no combination)
         The default is 'max'.
     use_derivative : bool, optional
         Whether the use the derivative of the gradient waveforms instead of the
@@ -57,7 +57,7 @@ def calculate_gradient_spectrum(
     -------
     spectrograms : List[np.ndarray]
         List of spectrograms per gradient channel.
-    spectrogram_sos : np.ndarray
+    spectrogram_rss : np.ndarray
         Root-sum-of-squares combined spectrogram over all gradient channels.
     frequencies : np.ndarray
         Frequency axis of the spectrograms.
@@ -72,7 +72,7 @@ def calculate_gradient_spectrum(
     # Get gradients as piecewise-polynomials
     gw_pp = obj.get_gradients(time_range=time_range)
     ng = len(gw_pp)
-    max_t = max(g.x[-1] for g in gw_pp if g != None)
+    max_t = max(g.x[-1] for g in gw_pp if g is not None)
     
     # Determine sampling points
     if time_range == None:
@@ -93,13 +93,12 @@ def calculate_gradient_spectrum(
         gw = np.diff(gw, axis=1)
     
     # Calculate spectrogram for each gradient channel
-    spectrograms = []
-    frequencies = None
-    spectrogram_sos = 0
+    spectrograms: List[np.ndarray] = []
+    spectrogram_rss = 0
     
     for i in range(ng):
         # Use scipy to calculate the spectrograms
-        freq, times, Sxx = spectrogram(gw[i],
+        freq, times, sxx = spectrogram(gw[i],
                                        fs=1/dt,
                                        mode='magnitude',
                                        nperseg=nwin,
@@ -110,31 +109,31 @@ def calculate_gradient_spectrum(
         mask = freq<max_frequency
         
         # Accumulate spectrum for all gradient channels
-        spectrogram_sos += Sxx[mask]**2
+        spectrogram_rss += sxx[mask]**2
         
         # Combine spectrogram over time axis
         if combine_mode == 'max':
-            s = Sxx[mask].max(axis=1)
+            s = sxx[mask].max(axis=1)
         elif combine_mode == 'mean':
-            s = Sxx[mask].mean(axis=1)
-        elif combine_mode == 'sos':
-            s = np.sqrt((Sxx[mask]**2).sum(axis=1))
+            s = sxx[mask].mean(axis=1)
+        elif combine_mode == 'rss':
+            s = np.sqrt((sxx[mask]**2).sum(axis=1))
         elif combine_mode == 'none':
-            s = Sxx[mask]
+            s = sxx[mask]
         else:
-            raise ValueError(f'Unknown value for combine_mode: {combine_mode}, must be one of [max, mean, sos, none]')
+            raise ValueError(f'Unknown value for combine_mode: {combine_mode}, must be one of [max, mean, rss, none]')
         
         frequencies = freq[mask]
         spectrograms.append(s)
     
     # Root-sum-of-squares combined spectrogram for all gradient channels
-    spectrogram_sos = np.sqrt(spectrogram_sos)
+    spectrogram_rss = np.sqrt(spectrogram_rss)
     if combine_mode == 'max':
-        spectrogram_sos = spectrogram_sos.max(axis=1)
+        spectrogram_rss = spectrogram_rss.max(axis=1)
     elif combine_mode == 'mean':
-        spectrogram_sos = spectrogram_sos.mean(axis=1)
-    elif combine_mode == 'sos':
-        spectrogram_sos = np.sqrt((spectrogram_sos**2).sum(axis=1))
+        spectrogram_rss = spectrogram_rss.mean(axis=1)
+    elif combine_mode == 'rss':
+        spectrogram_rss = np.sqrt((spectrogram_rss**2).sum(axis=1))
     
     # Plot spectrograms and acoustic resonances if specified
     if plot:
@@ -144,8 +143,8 @@ def calculate_gradient_spectrum(
             # According to spectrogram documentation y unit is (Hz/m)^2 / Hz = Hz/m^2, is this meaningful?
             for s in spectrograms:
                 plt.plot(frequencies, s)
-            plt.plot(frequencies, spectrogram_sos)
-            plt.legend(['x', 'y', 'z', 'sos'])
+            plt.plot(frequencies, spectrogram_rss)
+            plt.legend(['x', 'y', 'z', 'rss'])
     
             for res in acoustic_resonances:
                 plt.axvline(res['frequency'], color='k', linestyle='-')
@@ -169,7 +168,7 @@ def calculate_gradient_spectrum(
             plt.title('Total spectrum')
             plt.xlabel('Time (s)')
             plt.ylabel('Frequency (Hz)')
-            plt.imshow(abs(spectrogram_sos[::-1]), extent=(times[0], times[-1], frequencies[0], frequencies[-1]),
+            plt.imshow(abs(spectrogram_rss[::-1]), extent=(times[0], times[-1], frequencies[0], frequencies[-1]),
                        aspect=(times[-1]-times[0])/(frequencies[-1]-frequencies[0]))
             
             for res in acoustic_resonances:
@@ -177,4 +176,4 @@ def calculate_gradient_spectrum(
                 plt.axhline(res['frequency'] - res['bandwidth']/2, color='r', linestyle='--')
                 plt.axhline(res['frequency'] + res['bandwidth']/2, color='r', linestyle='--')
                 
-    return spectrograms, spectrogram_sos, frequencies, times
+    return spectrograms, spectrogram_rss, frequencies, times

From 26176251a83ff591971b9eead3bec7601ceced18 Mon Sep 17 00:00:00 2001
From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com>
Date: Fri, 15 Mar 2024 16:16:19 +0100
Subject: [PATCH 4/5] Minor changes in `calculate_gradient_spectrum` - Changed
 `time_range` type hint to `Union` to be compatible with older python versions
 - Changed initialization of `spectrogram_rss` from `0` to `np.array(0)`

---
 pypulseq/Sequence/calc_grad_spectrum.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
index 66d7cad..099d417 100644
--- a/pypulseq/Sequence/calc_grad_spectrum.py
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -1,4 +1,4 @@
-from typing import Tuple, List
+from typing import Tuple, List, Union
 
 import numpy as np
 from scipy.signal import spectrogram
@@ -10,7 +10,7 @@ def calculate_gradient_spectrum(
         max_frequency: float = 2000,
         window_width: float = 0.05,
         frequency_oversampling: float = 3,
-        time_range: List[float] | None = None,
+        time_range: Union[List[float], None] = None,
         plot: bool = True,
         combine_mode: str = 'max',
         use_derivative: bool = False,
@@ -94,7 +94,7 @@ def calculate_gradient_spectrum(
     
     # Calculate spectrogram for each gradient channel
     spectrograms: List[np.ndarray] = []
-    spectrogram_rss = 0
+    spectrogram_rss = np.array(0)
     
     for i in range(ng):
         # Use scipy to calculate the spectrograms

From fe8889203585bf740476c0eb36e163ba5e9e5f5f Mon Sep 17 00:00:00 2001
From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com>
Date: Fri, 15 Mar 2024 16:25:52 +0100
Subject: [PATCH 5/5] Reverted change in `spectrogram_rss` initialization in
 `calculate_gradient_spectrum`

---
 pypulseq/Sequence/calc_grad_spectrum.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
index 099d417..3172939 100644
--- a/pypulseq/Sequence/calc_grad_spectrum.py
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -94,7 +94,7 @@ def calculate_gradient_spectrum(
     
     # Calculate spectrogram for each gradient channel
     spectrograms: List[np.ndarray] = []
-    spectrogram_rss = np.array(0)
+    spectrogram_rss = 0
     
     for i in range(ng):
         # Use scipy to calculate the spectrograms