Usage

Running a Pulse Sequence

Pulse sequence programs are used programmatically via the Sequence class. Below is the code for a typical usage scenario. This code may be run directly in a jupyter notebook.

Load a pulse sequence by passing its file path:

from matipo import SEQUENCE_DIR, GLOBALS_DIR
from matipo.sequence import Sequence

# load system FID pulse sequence
fid_seq = Sequence(SEQUENCE_DIR+'FID.py')

Load parameter files by passing their file paths to the loadpar() method:

# load system frequency and pulse calibration files
fid_seq.loadpar(GLOBALS_DIR+'frequency.yaml')
fid_seq.loadpar(GLOBALS_DIR+'hardpulse_90.yaml')

Set parameters directly using the setpar() method:

fid_seq.setpar(
    n_scans=1,
    n_samples=10000,
    t_dw=10e-6,
    t_end=1
)

The current set of parameters can be read using the par property:

# print all parameters
print(fid_seq.par)
# use parameters to calculate total acquisition time: (dwell time)*(samples)
print('Acquisition time:', fid_seq.par.t_dw*fid_seq.par.n_samples)

Note

All parameters have a default value set in the pulse sequence program which will be used if not overriden by loadpar() or setpar().

Run the pulse sequence program using the run() method, which also returns the data. This is an async method, which means that await must be used to schedule it to run and obtain the result when finished:

data = await fid_seq.run()

The data from a sequence object’s last run can also be accessed later with the data property:

data = fid_seq.data

Data is returned as a numpy array of complex values. This data may be plotted with matplotlib:

import matplotlib.pyplot as plt
plt.plot(data.real, label='real')
plt.plot(data.imag, label='imag')
plt.legend()
plt.show()

Parameters may be saved to a file using the savepar() method:

fid_seq.savepar('FID_example.yaml')

This file can then be loaded with loadpar() to restore the same parameters at a later date.

Data may be saved using numpy:

import numpy as np
np.save('FID_example.npy', data)

Writing Custom Pulse Sequence Programs

Template

There is some boilerplate that all sequence programs must have. Here is a simple sequence which runs a single acquisition on a channel, and can be used as a template:

from matipo.sequence import ParDef
from collections import namedtuple
import numpy as np

# Define input parameters (must be named 'PARDEF')
PARDEF = [
    # args: name, dtype, default value
    ParDef('f', float, 1e6, unit='Hz'),
    ParDef('t_dw', float, 1e-6, min=0.1e-6, max=160e-6, unit='s'),
    ParDef('n_samples', int, 1000, min=2, unit=''),
    ParDef('rx_chan', int, 0, unit='')
]

# Create parameter set type (must be named 'ParameterSet')
ParameterSet = namedtuple('ParameterSet', [pd.name for pd in PARDEF])

# Main pulse sequence function, taking the sequencer object and a parameter set
def main(seq, par: ParameterSet):
    t_acq = par.n_samples * par.t_dw
    t_flatfilter_groupdelay = 10*par.t_dw

    yield seq.rx[par.rx_chan].freq(par.f)
    yield seq.rx[par.rx_chan].dwelltime(par.t_dw)
    yield seq.rx[par.rx_chan].mode(flatfilter=True)

    # Wait filter group delay after changing acquisition settings to avoid edge effects
    yield seq.wait(t_flatfilter_groupdelay)

    # Start acquisition with ID=0 and the desired number of samples
    yield seq.rx[par.rx_chan].acquire(0, par.n_samples)

    # Allow acquisition to finish completely before ending the sequence
    yield seq.wait(t_acq + t_flatfilter_groupdelay)

Input Parameter Definition

The input parameters are defined in the PARDEF list:

PARDEF = [
    # args: name, dtype, default value
    ParDef('f', float, 1e6, unit='Hz'),
    ParDef('t_dw', float, 1e-6, min=0.1e-6, max=160e-6, unit='s'),
    ParDef('n_samples', int, 1000, min=2, unit=''),
    ParDef('rx_chan', int, 0, unit='')
]

Where each parameter is defined with ParDef(<name>, <data type>, <default value>). min, max, and unit may also be defined, and will be read by experiments using BaseExperiment to set limits on the inputs. The Pint library is used to convert units if Pint Quantities are passed into the setpar() method.

The matipo.sequence.floatarray data type may be used to pass numpy arrays as parameters.

Main Function

The main() function is the actual pulse sequence program:

def main(seq, par: ParameterSet):
    t_acq = par.n_samples * par.t_dw
    t_flatfilter_groupdelay = 10*par.t_dw

    yield seq.rx[par.rx_chan].freq(par.f)
    yield seq.rx[par.rx_chan].dwelltime(par.t_dw)
    yield seq.rx[par.rx_chan].mode(flatfilter=True)

    # Wait filter group delay after changing acquisition settings to avoid edge effects
    yield seq.wait(t_flatfilter_groupdelay)

    # Start acquisition with ID=0 and the desired number of samples
    yield seq.rx[par.rx_chan].acquire(0, par.n_samples)

    # Allow acquisition to finish completely before ending the sequence
    yield seq.wait(t_acq + t_flatfilter_groupdelay)

This function is a python generator which returns pulse sequence commands using Python’s yield syntax.

Pulse Sequence Commands

Pulse sequence commands are methods of the Sequencer object and its children:

seq (SequencerTop):

wait(<time>)
join(<list of commands>)
tx[<channel index>] (RFTX)
- enable()
- disable()
- freq(<freq (Hz)>)
- phase(<phase (degrees)>)
- amp(<amplitude (ratio)>)
rx[<channel index>] (ACQ)
- acquire(<id>, <number of samples>)
- dwelltime(<sample dwell time (s)>)
- freq(<freq (Hz)>)
- phase(<phase (degrees)>)
- mode(flatfilter={True|False}, raw={True|False})
grad[<controller index>] (GRAD)
- vec(<x>, <y>, <z>)
- aux(<value>)
shim[<controller index>] (SHIM)
- set(<channel>, <value>)
gpo[<controller index>] (GPO)
- write(<value>)
- set(<bitmask>)
- clear(<bitmask>)

Composing Commands using Variables & Functions

A list of commands may be joined together with seq.join() and stored in a variable to be yielded multiple times for efficiency, e.g.:

pulse_90 = seq.join([
    seq.tx[0].freq(10e6),
    seq.tx[0].amp(1),
    seq.tx[0].enable(),
    seq.wait(10e-6),
    seq.tx[0].disable(),
    seq.tx[0].amp(0)
])
yield pulse_90
yield seq.wait(10e-6)
yield pulse_90
yield seq.wait(10e-6)
yield pulse_90
yield seq.wait(10e-6)

This code will result in three 10 microsecond full amplitude pulses at 10 MHz with 10 microseconds spacing.

Custom functions may also be defined using python syntax for convenience, e.g. to transmit an RF pulse with amplitude modulation:

def amplitude_modulated_pulse(seq, chan, amp_pts, duration):
    # calculate spacing of amplitude updates
    dt = duration/len(amp_pts)

    # start a list of commands with the TX enable command
    cmds = [seq.tx[chan].enable()]

    # for each amplitude value, append a command setting that amplitude and waiting for the spacing
    cmds += [
        seq.join([
            seq.tx[chan].amp(amp),
            seq.wait(dt)])
        for amp in amp_pts]

    # append TX disable command
    cmds += [seq.tx[chan].disable()]

    # join list of commands and return
    return seq.join(cmds)

# example usage, 1 millisecond 5 lobe sinc pulse:
def main(seq, par):
    ... snipped
    sinc_amp_pts = np.sinc(np.linspace(-3, 3, 100))
    yield amplitude_modulated_pulse(seq, 0, sinc_amp_pts, 1e-3)
    ... rest of sequence

Amplitude Ranges

Transmit amplitude is a value in the range -1.0 to 1.0, where 1 is the maximum possible amplitude of 0dBm.

Gradient and shim amplitudes are also defined in the range -1.0 to 1.0, which corresponds to an output range of -10V to 10V for the gradient channels and -5V to 5V for the shim channels.

RF pulse amplitude may be negative for convenience, which results in a 180 degree phase shift. This makes defining shaped pulses simpler as the phase can be kept constant and the amplitude varied with a waveform that contains negative values.

Timing Constraints

The minimum wait time is 100 nanoseconds (100e-9), but the timing precision is 10 nanoseconds. This means you may have a delay of exactly 110 nanoseconds, for example, but a wait time of 116 nanoseconds would be rounded to 120 nanoseconds.

There must be at least a 10e-6 (10 microseconds) delay between successive updates on a single gradient/shim controller. Note that you may still have multiple commands before a time delay as they will be combined together into a single update to the controller by the driver.

If a timing contraint violation occurs, they driver will terminate the sequence and return an error code.

Acquisition Modes/Filters

The signal processing pipeline is a Digital-Down-Converter (DDC) with three stages:

Frequency Mixer: Shifts the reference frame to the desired frequency and outputs quadrature data. The frequency and phase are set with the freq() and phase() commands.
CIC Decimator: Decimates the quadrature data to the desired dwell time, set by the dwelltime() command. The antialiasing filter has a fast settling time of 6 samples but does not have a flat frequency response.
FIR Compensation Filter: compensates for the response of the CIC filter to produce a flat frequency domain with a sharp rolloff. This filter has a settling time of approx. 20 samples.

The DDC may be disabled entirely by setting raw=True using the mode() command. In raw mode the output data will be the unprocessed 16 bit samples from the ADC.

Note

In raw mode, all acquisitions must have a number of samples that is a multiple of 4.

The flat filter may be disabled by setting flatfilter=False using the mode() command. This may be desirable in applications that require fast settling time.

General Purpose Outputs

The GPO state can be controlled using the gpo_set() and gpo_clear() commands. Here is an example of using a GPO to control an active T/R switch while transmitting an RF pulse:

TR_SWITCH_ENABLE = 0x00000001 # Port 3 pin 1 bitmask
yield seq.gpo[0].set(TR_SWITCH_ENABLE) # set high, T/R switch transmit mode
yield seq.tx[0].phase(0) # 0 degrees phase
yield seq.tx[0].amp(0.5) # 50% amplitude
yield seq.tx[0].enable() # enable RF output
yield seq.wait(100e-6) # wait 100 microseconds
yield seq.tx[0].disable() # disable RF output
yield seq.gpo[0].clear(TR_SWITCH_ENABLE) # clear low, T/R switch receive mode

GPO Pin Mapping

	GPO Port
Pin	0	1	2	3
1	0x01000000	0x00010000	0x00000100	0x00000001
2	0x02000000	0x00020000	0x00000200	0x00000002
3	0x04000000	0x00040000	0x00000400	0x00000004
4	0x08000000	0x00080000	0x00000800	0x00000008
5	0x10000000	0x00100000	0x00001000	0x00000010
6	0x20000000	0x00200000	0x00002000	0x00000020
7	0x40000000	0x00400000	0x00004000	0x00000040
8	0x80000000	0x00800000	0x00008000	0x00000080

Error Codes

If the pulse sequence execution driver detects an error, the pulse sequence will be terminated and the Sequence.run() method will raise an exception with an error code:

-131: RFTX frequency out of bounds
-132: RFTX amplitude out of bounds
-141: ACQ frequency out of bounds
-142: ACQ raw mode acquisition length not divisible by 4
-144: ACQ dwell time out of bounds
-145: ACQ number of samples out of bounds
-146: ACQ data too large for DMA buffer
-162: Sequencer wait time out of bounds (too small, must be >=100e-9 seconds)
-171: SHIM channel out of range
-1106: Input buffer ran out of pulse sequence commands before the end of the sequence. Usually may be remedied by optimising the sequence code, contact support for guidance.