Simulating Spin-Echo Metabolite NMR / MR Spectra with PyGamma (VeSPA)

PyGamma is a fantastic python library that provides a wrapper around the C++ library Gamma (General Approach to Magnetic resonance Mathematical Analysis). Simulating NMR or MRS spectra with PyGamma is easy, however, the examples that come bundled with PyGamma a little limited if you just want a quick simulation. So let us have a look at a modified version of the example FID pulse sequence in PyGamma to get our acquisition table mx (if you just want the code without the explanation, you can go here).

GABA simulation

Simulation output


First, lets set up a python virtual environment to run this, based on the requirements found with the code here. A range of systems are included under ‘sys/’, generally, these systems have the Govindaraju et al. values. There’s a range of 23 of metabolites in there, including NAA, Creatine, Glutamine, Glutamate, I’d double check these if you need specific or exact values.

virtualenv venv && source venv/bin/activate && pip install -r requirements.txt
import pygamma as pg

def fid(sys, b0=123.23):
    H = pg.Hcs(sys) + pg.HJ(sys)
    D = pg.Fm(sys)
    ACQ = pg.acquire1D(pg.gen_op(D), H, 0.1)
    sigma = pg.sigma_eq(sys)
    sigma0 = pg.Ixpuls(sys, sigma, 90.0)
    mx = pg.TTable1D(ACQ.table(sigma0))
    return mx

This is a modified version of the example that can be found here. There’s a couple of additional parts added here, first we’ve wrapped the ACQ.table in a pg.TTable1D() to stop the C table object being garbage collected when we exit the fid function. This means we can pass the mx table object around between methods. Secondly, I’ve added the ability to adjust the b0 frequency (in MHz) of the scanner via Omega adjust.

Next, we need to take the mx object and get it to return our FID, FFT and nu. If you’re interested in the methods invoked here, check out the L2 section of the Gamma user guide here.

def binning(mx, b0=123.23, high_ppm=-3.5, low_ppm=-1.5, dt=5e-4, npts=4096, linewidth=1):
    mx.broaden(linewidth)       # add some broadening
    mx.resolution(0.5)          # Combine transitions within 0.5 rad/sec

    ADC = []
    acq = mx.T(npts, dt)
    for ii in range(0, acq.pts()):
        ADC.extend([acq.get(ii).real() + (1j * acq.get(ii).imag())])

    FFT = []
    acq = mx.F(npts, b0 * high_ppm, b0 * low_ppm)
    for ii in range(0, acq.pts()):
        FFT.extend([acq.get(ii).real() + (1j * acq.get(ii).imag())])

    ADC = normalise_complex_signal(np.array(ADC))
    FFT = normalise_complex_signal(np.array(FFT))
    nu = np.linspace(high_ppm, low_ppm, npts)

    return ADC, FFT, nu

The FFT returned by PyGamma is not the same as FFT(ADC), rather than having a bandwidth determined by acquisition parameters, it has the exact frequency range specified by low_ppm and high_ppm. There’s also small function included for normalizing the complex signals to be within the range [-1:1]. This makes plotting our data much nicer:

def normalise_complex_signal(signal):
    return signal / np.max([np.abs(np.real(signal)), np.abs(np.imag(signal))])

Finally, let’s build some plots for this data.

def plot(sys, ADC, FFT, nu):
    fig, axes = plt.subplots(2, 2, sharex='col', sharey=True,  figsize=(19.2, 10.8), dpi=85)
    # pulse sequence name, system name and magnetic field strength
    plt.suptitle('FID: %s B0: %.2fT' %(, hz_to_tesla(sys.Omega())))

    # Plot the reported FFT from PyGamma
    axes[0, 0].set_title('FFT Magnitude spectrum')
    axes[0, 0].plot(nu, np.abs(FFT)/max(np.abs(FFT)))
    axes[0, 0].set_xlim([min(nu), max(nu)])

    axes[1, 0].set_title('FFT Real & Imag. spectrum')
    axes[1, 0].plot(nu, np.real(FFT), label='Real')
    axes[1, 0].plot(nu, np.imag(FFT), label='Imag')
    axes[1, 0].set_xlabel('PPM')
    axes[1, 0].legend(['Real', 'Imaginary'])

    # Plot the ADC
    axes[0, 1].set_title('ADC Real')
    axes[0, 1].plot(np.linspace(0, len(ADC), len(ADC)), np.real(ADC))
    axes[0, 1].set_xlim([0, len(ADC)])

    axes[1, 1].set_title('ADC Imaginary')
    axes[1, 1].plot(np.linspace(0, len(ADC), len(ADC)),np.imag(ADC))
    axes[1, 1].set_xlabel('T')

Loading, simulation and binning for GABA a 6 spin system takes 0.03 seconds on my dual-core laptop.

Code for this and all my PyGamma examples can be found on GitHub.


PfSense DNS Resolver with PiHole DNS forwarder for network-wide ad-blocking

Configuring Unbound on PfSense as a DNS Resolver to register DHCP hostnames on localdomain and using PiHole as a DNS Forwarder to perform network-wide ad-blocking.

Setting up a Greylog server for central logging with Ansible

Learning how to use Ansible through setting up a Graylog server for centralised syslogging.

Proxmox Setup v2 - moving root onto an SSD and away from the ZFS array

Moving away from ZFS on root, and using it for the HDD array.

Proxmox host migration; new Home Server day

A workbook of migrating Proxmox to a new host for the first time.

Migrating Docker from Digital Ocean to home

Moving my docker containers to a new host at home.

Setting up a VLAN with PfSence

A quick introduction to setting up VLANs with PfSence for a guest network.

Converting a BBC Micro keyboard to USB with an Arduino Micro (ATmega 32U4)

Using an Arduino to turn a BBC Micro keyboard to be USB compatible.