Shinnar-Le Roux Pulse design demo

The previous section showed how to simulate the profile of an RF pulse. In this section we show how to design much better RF pulses.

The SLR design starts by designing the beta polynomial, which is a lowpass function scalled to sine of half the flip angle. For an inversion pulse, this is sin(pi/2) = 1;

We start with a windowed sinc again

%use octave

try
    cd rf_tools_octave
catch
    try
        cd ../rf_tools_octave
    catch
        try
            cd ../rf_tools_octave
        catch
            cd ../../rf_tools_octave
        end
    end
end

pkg load signal

%use octave

b = msinc(256,2);

This is scaled to 1, which is what we want. We then design a consistent $\alpha$, and then solve for the RF pulse. $\alpha$is determined uniquely given that we want a minimum power pulse, so the result is

%use octave

rf = b2rf(b);

We will compare this to the windowed sinc scaled to $\pi$,

%use octave

rfw = b*pi;

Compute the two responses as inversions.

%use octave

x = [-64:64]/4;
mz = ab2inv(abr(rf,x));
mzw = ab2inv(abr(rfw,x));

%use sos
%get x --from Octave
%get mz --from Octave
%get mzw --from Octave

import matplotlib.pyplot as plt
import plotly.plotly as py
import plotly.graph_objs as go
import numpy as np
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
config={'showLink': False, 'displayModeBar': False}

init_notebook_mode(connected=True)

from IPython.core.display import display, HTML

re = go.Scatter(
    x = x,
    y = mz,
    name = 'slr',
    text = 'slr',
    hoverinfo = 'x+y+text', 
    line = dict(
        color = ('rgb(22, 96, 167)'),
        ),
)

im = go.Scatter(
    x = x,
    y = mzw,
    name = 'windowed sinc',
    text = 'windowed sinc',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(205, 12, 24)'),
        ),
)

data = [re, im]

layout = go.Layout(
    width=600,
    height=350,
    margin=go.layout.Margin(
        l=100,
        r=50,
        b=60,
        t=0,
    ),
    xaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    yaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    legend=dict(
        x=0.70,
        y=0.65,
        traceorder='normal',
        font=dict(
            family='Times New Roman',
            size=12,
            color='#000'
        ),
        bordercolor='#000000',
        borderwidth=2
    )
)

fig = dict(data=data, layout=layout)

plot(fig, filename = 'slr_fig_1.html', config = config)

display(HTML('slr_fig_1.html'))

(View plot code)

Note the sharper profile of the SLR pulse. This comes at a cost if we compare RF pulses

%use octave

t = [1:256]/32;

rfscale_rf = real(rfscale(rf,8));
rfscale_rfw = real(rfscale(rfw,8));

%use sos
%get t --from Octave
%get rfscale_rf --from Octave
%get rfscale_rfw --from Octave

re = go.Scatter(
    x = t,
    y = rfscale_rf,
    name = 'slr',
    text = 'slr',
    hoverinfo = 'x+y+text', 
    line = dict(
        color = ('rgb(22, 96, 167)'),
        ),
)

im = go.Scatter(
    x = t,
    y = rfscale_rfw,
    name = 'windowed sinc',
    text = 'windowed sinc',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(205, 12, 24)'),
        ),
)

data = [re, im]

layout = go.Layout(
    width=600,
    height=350,
    margin=go.layout.Margin(
        l=100,
        r=50,
        b=60,
        t=0,
    ),
    xaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    yaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2
    ),
    annotations=[
        dict(
            x=0.5004254919715793,
            y=-0.18,
            showarrow=False,
            text='time, ms',
            font=dict(
                family='Times New Roman',
                size=22
            ),
            xref='paper',
            yref='paper'
        ),
        dict(
            x=-0.15,
            y=0.5,
            showarrow=False,
            text='amplitude, kHz',
            font=dict(
                family='Times New Roman',
                size=22
            ),
            textangle=-90,
            xref='paper',
            yref='paper'
        ),
    ],
    legend=dict(
        x=0.70,
        y=0.65,
        traceorder='normal',
        font=dict(
            family='Times New Roman',
            size=12,
            color='#000'
        ),
        bordercolor='#000000',
        borderwidth=2
    )
)

fig = dict(data=data, layout=layout)


plot(fig, filename = 'slr_fig_2.html', config = config)

display(HTML('slr_fig_2.html'))

(View plot code)

rfscale(rf,t) takes an RF pulse that sums to its flip angle, and a pulse duration in ms, and scales it to an RF amplitude in kHz.

The SLR pulse has almost twice the peak amplitude. Note that the SLR pulse has closer zero spacing, even though it was based on the windowed sinc waveform.

%use octave

help dzrf

'dzrf' is a function from the file /home/jovyan/work/RF-Tools-in-Octave/rf_tools_octave/dzrf.m

   rf = dzrf(np,tb,ptype,ftype,d1,d2)

  Designs an rf pulse.  There are a lot of options, most of
  which have defaults.  For example, a reasonable 100 sample
  tb=4 spin-echo pulse can be designed with

   rf = dzrf(100,4,'se')

  Inputs are:
    np -- number of points.         (required)
    tb -- time-bandwidth product    (required)
    ptype -- pulse type.  Options are:
      st  -- small tip angle         (default)
      ex  -- pi/2 excitation pulse
      se  -- pi spin-echo pulse
      sat -- pi/2 saturation pulse
      inv -- inversion pulse
    ftype -- filter design method.  Options are:
      ms  -- Hamming windowed sinc (an msinc)
      pm  -- Parks-McClellan equal ripple
      ls  -- Least Squares           (default)
      min -- Minimum phase (factored pm)
      max -- Maximum phase (reversed min)
    d1 -- Passband ripple        (default = 0.01)
    d2 -- Stopband ripple        (default = 0.01)


Additional help for built-in functions and operators is
available in the online version of the manual.  Use the command
'doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at http://www.octave.org and via the help@octave.org
mailing list.

For example, to design a very low ripple inversion pulse using an equal ripple Parks-McClellan design,

%use octave

rf = dzrf(256,4,'inv','pm',0.001,0.001);

re_rf = real(rf);
im_rf = imag(rf);

ab2inv_rf = ab2inv(abr(rf,x));

re_ab2inv_rf = real(ab2inv_rf);
im_ab2inv_rf = imag(ab2inv_rf);

%use sos
%get re_rf --from Octave
%get im_rf --from Octave
%get re_ab2inv_rf --from Octave
%get im_ab2inv_rf --from Octave

sub1_1 = go.Scatter(
    x = [i for i in range(len(re_rf))],
    y = re_rf,
    name = 'real',
    text = 'real',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(22, 96, 167)'),
        ),
    showlegend=False,
)

sub1_2 = go.Scatter(
    x = [i for i in range(len(im_rf))],
    y = im_rf,
    name = 'imaginary',
    text = 'imaginary',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(205, 12, 24)'),
        ),
    showlegend=False,
)

sub2_1 = go.Scatter(
    x = [i for i in range(len(re_ab2inv_rf))],
    y = re_ab2inv_rf,
    name = 'real',
    text = 'real',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(22, 96, 167)'),
        ),
)

sub2_2 = go.Scatter(
    x = [i for i in range(len(im_ab2inv_rf))],
    y = im_ab2inv_rf,
    name = 'imaginary',
    text = 'imaginary',
    hoverinfo = 'x+y+text',
    line = dict(
        color = ('rgb(205, 12, 24)'),
        ),
)


from plotly import tools

fig = tools.make_subplots(rows=2, cols=1, print_grid=False)
fig.append_trace(sub1_1, 1, 1)
fig.append_trace(sub1_2, 1, 1)

fig.append_trace(sub2_1, 2, 1)
fig.append_trace(sub2_2, 2, 1)

layout_subplot = go.Layout(
    width=600,
    height=600,
    margin=go.layout.Margin(
        l=100,
        r=50,
        b=60,
        t=0,
    ),
    xaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2,
        domain=[0, 1]
    ),
    yaxis=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2,
        domain=[0.55, 1],
    ),
    xaxis2=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2,
        domain=[0, 1],
        anchor='y2'
    ),
    yaxis2=dict(
        showgrid=False,
        linecolor='black',
        linewidth=2,
        domain=[0, 0.45],
    ),
    legend=dict(
        x=0.75,
        y=0.9,
        traceorder='normal',
        font=dict(
            family='Times New Roman',
            size=12,
            color='#000'
        ),
        bordercolor='#000000',
        borderwidth=2
    )
)

fig["layout"]=layout_subplot

config={'showLink': False, 'displayModeBar': False}

plot(fig, filename = 'slr_fig_3.html', config = config)

display(HTML('slr_fig_3.html'))

(View plot code)

Something to try:

  • Look at how the RF pulses and profiles chance as the flip angle decreases
  • Scale $\beta$ to $\sin(\pi0.99/2)$, $\sin(\pi0.95/2)$, and $\sin(\pi*0.9/2)$ using the msinc filter from above.