Sequence.gradient_waveforms() does not interpolate trapezoids with gradient raster time.

[btasdelen]

Describe the bug
The issue arises when a trapezoid is defined by time points and amplitudes. gradient_waveforms() function directly puts gradient's waveform variable into the gradient axis without considering the timings, which results in incorrect rendering of the gradient. The responsible part seems to be in Sequence.py line 667:

grad_waveforms[j, l1:l2] = waveform

Ideally, there should be a code that interpolates the waveform using grad.tt to gradient raster time before the aforementioned line.

To Reproduce
Here is a minimal sequence that reproduces the bug:

import matplotlib.pyplot as plt
import numpy as np
from pypulseq.make_extended_trapezoid_area import make_extended_trapezoid_area
from pypulseq.opts import Opts
from pypulseq.Sequence.sequence import Sequence

# ## **SYSTEM LIMITS**
# Set the hardware limits and initialize sequence object

system = Opts(max_grad=26, grad_unit='mT/m', max_slew=45, slew_unit='T/m/s', 
              grad_raster_time=10e-6, rf_ringdown_time=10e-6, 
              rf_dead_time=100e-6)
seq = Sequence(system)

gz_spoil,gz_spoil_times,gz_spoil_amps = make_extended_trapezoid_area(channel='z', grad_start=0, grad_end=0, system=system, area=3536)

seq.add_block(gz_spoil)
            
## Waveform debug

gw = seq.gradient_waveforms()
gws = (gw[:, 1:] - gw[:, :-1]) / seq.system.grad_raster_time
gs = np.max(np.abs(gws), axis=1)/system.gamma

Expected behavior
The function should return correct gradient waveforms.

Screenshots
Here is a figure showing how the gradient is rendered. This also causes incorrect calculation of slew rate.

Desktop (please complete the following information):

• OS: macOS
• pypulseq version: dev branch

[btasdelen]

Workaround: Passing convert_to_arbitrary=True to make_extended_trapezoid() seems to work around the issue, naturally, but kind of defeats the purpose of extended trapezoids. Also, that parameter is not exposed to make_extended_trapezoid_area().

[sravan953]

Hi @bilal-tasdelen . Will take a look soon at the two issues you've opened, thanks!

[btasdelen]

If you are interested, I can send PRs?

[sravan953]

Yes PRs are most welcome! Definitely easier to maintain this framework with more helping hands :)

[btasdelen]

@sravan953 I see that both arbitrary gradients and extended trapezoids have the same type, 'grad'. Arbitrary grads assumed to be in GRT, whereas extended trapezoids are not. To fix the issue, I see two options.

1. We can define a separate type for extended trapezoids, e.g. 'ext_trap'. This way we can easily distinguish the types. An aside, type name 'grad' is not descriptive either.
2. We can check if the time axis of the gradients is on the GRT. I believe this is how it is done in 'waveforms_and_times()' function.
   Which one do you prefer?

[sravan953]

Thus far I have tried to keep PyPulseq identical to Pulseq in terms of code and structure so that users who use both platforms, or come from MATLAB to Python, don't find it jarring. Therefore, I would prefer method #2.

However, I believe waveforms_and_times() would be the more appropriate method to use here? It returns gradient values and their respective timepoints.

[btasdelen]

Got it. I agree that overall waveforms_and_times() is more useful, but there are a couple of cases where gradients directly returned on GRT is more convenient. One such example is detecting gradient resonance frequencies, where we take the FT of the entire waveform and check if there are peaks near the forbidden frequencies. Another case is gradient moment calculations.

I will check how it is implemented in the Pulseq, and make necessary changes accordingly.

Edit: I realized that function is dysfunctional in Pulseq. I believe the correct way to implement this function is to call waveforms_and_times(), and then call points_to_waveform() on the output to raster the waveforms. It is best this way because it will reuse the existing functions and reduce the risk of introducing bugs.

However, if you also want to keep the implementations the same, as well as the interface, I will understand. In that case, users can create their own convenience functions as I described.

[sravan953]

Thank you for sharing example use cases, I was not aware. Your suggestion sounds like a good idea -- and I agree with reusing code to reduce the risk of introducing bugs. Have you tried chaining the two methods and does it fix the issue?

[btasdelen]

I already did that for my fork. I did not observe any issues with my implementation yet. Creating a PR soon.