diff --git a/pypulseq/Sequence/calc_grad_spectrum.py b/pypulseq/Sequence/calc_grad_spectrum.py
new file mode 100644
index 00000000..3172939d
--- /dev/null
+++ b/pypulseq/Sequence/calc_grad_spectrum.py
@@ -0,0 +1,179 @@
+from typing import Tuple, List, Union
+
+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: Union[List[float], None] = 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]:
+    """
+    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', '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
+        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 [].
+
+    Returns
+    -------
+    spectrograms : List[np.ndarray]
+        List of spectrograms per gradient channel.
+    spectrogram_rss : 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 is not 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)
+    
+    if use_derivative:
+        gw = np.diff(gw, axis=1)
+    
+    # Calculate spectrogram for each gradient channel
+    spectrograms: List[np.ndarray] = []
+    spectrogram_rss = 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_rss += 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 == 'rss':
+            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, rss, none]')
+        
+        frequencies = freq[mask]
+        spectrograms.append(s)
+    
+    # Root-sum-of-squares combined spectrogram for all gradient channels
+    spectrogram_rss = np.sqrt(spectrogram_rss)
+    if combine_mode == 'max':
+        spectrogram_rss = spectrogram_rss.max(axis=1)
+    elif combine_mode == 'mean':
+        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:
+        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_rss)
+            plt.legend(['x', 'y', 'z', 'rss'])
+    
+            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_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:
+                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_rss, frequencies, times
diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py
index 258e736e..d7afbb55 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,75 @@ 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',
+            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
+        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'.
+        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 [].
+
+        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,
+                                           use_derivative=use_derivative,
+                                           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 72568b64..a38a0655 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.