Contents

This simulator demonstrates the interactions between slice thickness, field of view, te and bandwidth

You will find, while using the MR scanner, that some of the parameters that you choose for contrast, or resolution, interact with each other so that all combinations are not possible. This is because each of the events in the MR pulse sequence take some time, and because physical limitations of the equipment determine how fast these events can happen.

In this simulation, you will see how changing slice thickness, field of view, and field strength interact to determine the minimum te, and how they determine signal to noise ratio through changes in bandwidth.

The simulator will output the corresponding waveforms for slice selection and readout.

Enter your default parameters on lines 45 though 47, and the device parameters on lines 53 through 64, or use options on the command line.

Read the comments carefully to understand what is going on.

function SequenceSim(varargin)
if nargin < 1
    fprintf('USAGE: SequenceSim(OPTIONS)\n');
    fprintf('If no options are entered, default values are used\n');
    fprintf('OPTIONS\n');
    fprintf('     ''matrix''   - matrix size in pixels\n');
    fprintf('     ''fov''      - field of view in mm\n');
    fprintf('     ''slthick''  - Slice thickness in millimeters\n');
    fprintf('     ''field''    - Field strength in Tesla\n');
    fprintf('     ''slew''     - Gradient slew rate in milliTesla/m/msec\n');
    fprintf('     ''gmax''     - Maximum gradient amplitude in milliTesla/m\n');
    fprintf('     ''tr''       - Repetition time in milliseconds\n');
    fprintf('     ''nerd''     - Show calculated details\n');
    fprintf('For example:\n\tSequenceSim(''matrix'', 128, ''gmax'', 10, ''tr'', 1000);\n')
end

clf;
USAGE: SequenceSim(OPTIONS)
If no options are entered, default values are used
OPTIONS
     'matrix'   - matrix size in pixels
     'fov'      - field of view in mm
     'slthick'  - Slice thickness in millimeters
     'field'    - Field strength in Tesla
     'slew'     - Gradient slew rate in milliTesla/m/msec
     'gmax'     - Maximum gradient amplitude in milliTesla/m
     'tr'       - Repetition time in milliseconds
     'nerd'     - Show calculated details
For example:
	SequenceSim('matrix', 128, 'gmax', 10, 'tr', 1000);

User Parameters

To use this simulator, choose a matrix size, field of view, and slice thickness here:

If you wish to see how this plays out in technical details, set 'ShowNerdyDetails = true;'

ShowNerdyDetails = false;

Matrix = 256;   % Typically between 64 and 1024
FoVmm  = 192;   % Typically between 80 and 300
SliceThickness = 2e-3; % meters. Typically between 1e-3 and 5e-3
tr = 2000; % tr in milliseconds

Device Parameters

While you can't actually change these things once you have bought a scanner, these are key specs for the instruments you might buy. See their effects by changing them here.

FieldStrength = 3; % Tesla

Beyond field strength, each MR instrument has a series of specs that determine its ultimate performance. The gradient performance typically is measured in terms of maximum gradient strength, and the speed with which the gradient amplitudes can change (slew rate). The product of the gradient amplitude by time is the ultimate determinant of resolution, so larger, faster, gradients tend to outperform smaller, slower systems. We discuss the product of the gradient amplitude (in T/m) by time to be the gradient area

SlewRateMT = 20; % milliTesla/m/msec. Typically from 10 to 25
GradMaxMT  = 23; % milliTesla/m. Typically from 10 to 30

% Override with user-selected parameters
for j = 1:2:length(varargin)
    name = lower(varargin{j});  % make the options case insensitive
    value = varargin{j+1};
    switch name
        case 'matrix'
            Matrix = value;
        case 'fov'
            FoVmm = value;
        case 'slthick'
            SliceThickness = value/1000;
        case 'field'
            FieldStrength = value;
        case 'slew'
            SlewRateMT = value;
        case 'gmax'
            GradMaxMT = value;
        case 'tr'
            tr = value;
        case 'nerd'
            ShowNerdyDetails = true;
        otherwise
    end
end

SlewRate = SlewRateMT/1000; % convert units to Tesla/m/msec
GradMax  = GradMaxMT/1000;  % convert units to Tesla/m
riseTime = GradMax/SlewRate; % Gradient rise time to max

fprintf('\nSEQUENCE PARAMETERS\n');
fprintf('\tMatrix Size (matrix):.............. %d\n', Matrix);
fprintf('\tField of View (fov):............... %d mm\n', FoVmm);
fprintf('\tSlice Thickness (slthick):......... %d mm\n', SliceThickness * 1000);
fprintf('\ttr (tr):........................... %0.2f mm\n', tr);
fprintf('\nINSTRUMENT SPECIFICATIONS\n');
fprintf('\tField Strength (field):............ %0.1f Tesla\n', FieldStrength);
fprintf('\tMax Gradient Strength (gmax):...... %0.1f milliTesla/meter\n', GradMaxMT);
fprintf('\tGradient Slew rate (slew):......... %0.1f milliTesla/meter/msec\n', GradMaxMT);
SEQUENCE PARAMETERS
	Matrix Size (matrix):.............. 256
	Field of View (fov):............... 192 mm
	Slice Thickness (slthick):......... 2 mm
	tr (tr):........................... 2000.00 mm

INSTRUMENT SPECIFICATIONS
	Field Strength (field):............ 3.0 Tesla
	Max Gradient Strength (gmax):...... 23.0 milliTesla/meter
	Gradient Slew rate (slew):......... 23.0 milliTesla/meter/msec

Other constraints

The final resolution is calculated as resolution = 1/(Gamma * GradAmp * t) where GradAmp is the amplitude of the gradient in use (T/m). Therefore:

FoV        = FoVmm/1000; % meters
Resolution = FoV/Matrix; % meters

% Physics
Gamma    = 42.5774E6;   % Gyromagnetic Ratio Hz/Tesla
GammaMT  = Gamma/1000; % Gyromagnetic Ratio in Hz/millTesla

Readout Waveform

The readout gradient shape is trapezoidal, reaching its maximum (GradMax) at a maximum in a time determined by slew rate (SlewMax). The area of the ramps becomes part of the total readout area.

%-
%
%          -------------      ......... GFinal
%         /             \
%        /               \
%       /  |           |  \
%      /   |           |   \
%     /    |           |    \
%    /     |           |     \
%   /      |           |      \
%---       |           |       -----
%   |<ramp>|<-flatTop->|<ramp>|

The total area of the gradient is GFinal * (flatTop + 2*ramp).

For this simulation, we assume that we are reading out the signal ONLY when the readout gradient has reached its plateau amplitude

An alternative involves "ramp sampling" in which the TOTAL gradient area is used for spatial encoding. To establish the gradient timing in that case, we find the minimum duration waveform that has our desired readout area.

% The total area of the gradient pulse (amplitude * time) determines the spatial resolution as below:
ReadoutArea = 1/(Resolution*GammaMT); % (T/m)*sec
flatTop = ReadoutArea/(GradMax);
BandWidth = 2000/flatTop;
ChemShift = FieldStrength * Gamma * 2.3e-6;

if ShowNerdyDetails
    fprintf('\n=== Readout pulse ===\n');

    fprintf('    Readout Area = %f (T/m)*sec', ReadoutArea);
    fprintf('         flatTop = %0.3f msec\n',flatTop);
    fprintf('        riseTime = %0.3f msec\n',riseTime);

The total duration of the readout period is the sum of the flatTop and the two ramp times:

    fprintf('Readout duration = %0.3f msec\n',2*riseTime + flatTop);

What the gradients do is to cause the magnetic field, and therefore the MR frequency, to vary with position. 'BANDWIDTH" is a measurement of the frequencies over which the MR signal is spread by these gradients. It is convenient to think about this in terms of frequency difference per pixel.

NOTE: In our case we are trying to use the gradients at their maximum amplitude all the time, maxing out the bandwidth (at a cost in SNR) The factor of 2000 in the formula below converts milliseconds to frequency, and includes factor of two to account for the fact the frequencies are above and below the center frequency.

    fprintf('Bandwidth at this readout setting = %0.f Hz/pixel\n',BandWidth);
    fprintf('Typical readout flatTop: %0.f ms\n', 2000/(FieldStrength * Gamma * 2.3e-6));
end

Usually, the scanner will suggest a more reasonable bandwidth, on the order of Field_Strength * Gamma * Fat-Water_shift

Readout pre-encode

Generally, we center the signal in the readout period by adding a negative going pulse whose area is established to equal the area of the readout itself at the middle of the flatTop period. This is usually set equal to te.

%-
%
%            ------
%           /      \
% ---  pre /  read  \-----
%    \    /
%     \__/
%
%..............|..........
%              te
%
APreEncode = GradMax*(riseTime/2) + ReadoutArea/2; % Pre-encode area. You might want to check this

if ShowNerdyDetails
    fprintf('\n=== Pre-encode calculation ===\n');
    fprintf('pre-encode area = %f (T/m)*sec\n',APreEncode);
%
    fprintf('\n --> You might notice that the pre-encode area is LARGER than the readout area. Why is this?\n');
end

We are not collecting data during the pre-encode period, so any waveshape is possible. If the pre-encode area is smaller than the area of two complete ramps, we will use a triangular pulse shape (no flat top) whose maximum amplitude is less than Gmax.

%-
%
%   -----      -----
%        \    /
%         \  /
%          \/
%
%     -->| d  |<--

Letting d indicate the duration, A denote the area, and S indicate the slew rate, a little algebra shows that:

d = 2sqrt(A/s)
minAreaGmax = riseTime * GradMax;
if minAreaGmax < APreEncode
    preFlatTop = (APreEncode - minAreaGmax)/GradMax;
    preEncodeDur = 2*riseTime + preFlatTop;
    if ShowNerdyDetails
        fprintf('\nThe pre-encode shape will be a trapezoid.\n');
    end
else
    preEncodeDur =  2*sqrt(APreEncode/SlewRate);
    GPreEncode   = SlewRate * preEncodeDur/2; % The peak amplitude of the pre-encode gradient
    if ShowNerdyDetails
        fprintf('\nThe pre-encode shape will be triangular.\n');
    end
end
% we will plot this waveform below

Slice Selection

Slice selection is performed by transmitting an RF whose spectral width (the range of frequencies it contains) is equal to the range of frequencies covered in a slice. The latter is created by applying a gradient pulse. For example, if a gradient of 1 gauss/cm is applied (10 mT/m) the frequency differs by 4258 Hz/cm. Selecting a 1 mm slice would require an RF pulse whose spectral width is 425.8 Hz.

However, narrowing the width of an RF pulse necessitates lengthening it in time. We consider an RF pulse to have a dimensionless parameter, TBW (the time-bandwidth product) that allows one to trade time against spectral width. For a typical windowed sinc pulse, TBW might be 4.

Therefore, a pulse of 3.2 msec duration would have a bandwidth of 4/3.2e-3

TBW=4;
NominalRFPulseDur = 3.2e-3;
RFBandwidth = TBW/NominalRFPulseDur; % typical spectral width of an RF pulse = 1250Hz

It is a simple matter to calculate the maximum change in frequency with position in the magnet, as it is just the maximum gradient amplitude multipled by the gyromagnetic ratio.

MaxFreqGradient = GradMax * Gamma; % Hertz/m
MinSliceWidth = RFBandwidth / MaxFreqGradient;

if ShowNerdyDetails
    fprintf('\n                  Maximum Slice Gradient = %0.2f Hz/millimeter\n',MaxFreqGradient/1000);
    fprintf('Minimum slice thickness for 3.2 ms pulse = %0.3f mm\n',MinSliceWidth*1000);
end

The only way to make a narrower slice, given a physically limited maximum gradient amplitude, is to increase the RF duration to lower its bandwidth

In practice, the scanners make THICKER slices, by lowering hte gradient amplitude while attempting to keep the RF slice selection pulse bandwidth constant.

We consider te to be equal to the time between the center of the RF pulse and the center of the readout period. Because dephasing occurs whenever the gradients are on, there must be a compensating gradient applied after the RF pulse. Its area is about 0.6 times the area under the total RF.

%-
%
%       ----
%  RF   |  |
%    ---    --------------
%
%       ----
%  GZ  /    \
%     /      \
%    -        \    /-----
%              \  /
%               --

Suppose that the user desires a given slice thickness (SliceThickness defined above):

RequiredRFBW = MaxFreqGradient * SliceThickness; % Required RF bandwidth
RFPulseDur = TBW/RequiredRFBW; % RF pulse duration
% Let's build the gradient waveforms, changing to units of msec:
RFPulseDur = RFPulseDur * 1000;

if ShowNerdyDetails
    fprintf('              Required RF pulse duration = %0.2f ms\n', RFPulseDur);
end

% Create the slice select gradient waveform by making a table of (time, amplitude)
% pairs. Here I use pt to indicate time and k as an array index.
pt = 0; k = 1;
RFGradientPulseWaveform(k,:) = [0,0];
k = k+1; pt = pt + riseTime;
RFCenter = pt + RFPulseDur/2; % te is calculated from the center of the RF pulse.
RFGradientPulseWaveform(k,:) = [pt,GradMax]; % Begin RF pulse
k = k+1; pt = pt + RFPulseDur;
RFGradientPulseWaveform(k,:) = [pt,GradMax];% End RF pulse
ReadoutStartTime = pt; % At this point, we can start the readout pre-encode
k = k+1; pt = pt + riseTime; % Rise time and fall time are equal
RFGradientPulseWaveform(k,:) = [pt,0]; % At the end of the fall time, the gradient returns to zero

Slice select refocus

A_RFRefocus = 0.6 * RFPulseDur * GradMax;

Similar to the readout pre-encode pulses, the slice selection refocus area may result in either a triangular or a trapezoidal pulse shape depending on the total area required. Using formulae from the above on readout:

if minAreaGmax < A_RFRefocus
    RefocusFlatTop = (A_RFRefocus - minAreaGmax)/GradMax;
    RefocusDur = 2*riseTime + RefocusFlatTop;
    if ShowNerdyDetails
        fprintf('\nThe slice refocus pulse will be a trapezoid\n');
    end
else
    RefocusDur = 2*sqrt(A_RFRefocus/SlewRate);
    GRefocus   = SlewRate * RefocusDur/2; % The peak amplitude of the pre-encode gradient
    if ShowNerdyDetails
        fprintf('\nThe slice refocus pulse will be triangular\n');
    end
end

k= k+1;
if minAreaGmax < A_RFRefocus % trapezoidal
    pt = pt + riseTime;
    RFGradientPulseWaveform(k,:) = [pt, -GradMax];
    k = k+1; pt = pt+RefocusFlatTop;
    RFGradientPulseWaveform(k,:) = [pt, -GradMax];
    k = k+1; pt = pt+riseTime;
    RFGradientPulseWaveform(k,:) = [pt, 0];
else  % triangular
    pt = pt + RefocusDur/2;
    RFGradientPulseWaveform(k,:) = [pt,-GRefocus];
    k = k+1; pt = pt + RefocusDur/2;
    RFGradientPulseWaveform(k,:) = [pt,0];
end

% We cannot start the actual readout until slice refocus is complete.
MinReadoutStartTime = pt;

Plotted outputs

Create the readout gradient waveform in the same manner.

clf;
ReadoutPulseWaveform(1,:) = [0,0];
pt = ReadoutStartTime;
ReadoutPulseWaveform(2,:) = [pt,0];
k = 3;
if minAreaGmax < APreEncode % trapezoidal
    pt = pt + riseTime;
    ReadoutPulseWaveform(k,:) = [pt, -GradMax];
    k = k+1; pt = pt+preFlatTop;
    ReadoutPulseWaveform(k,:) = [pt, -GradMax];
    k = k+1; pt = pt+riseTime;
    ReadoutPulseWaveform(k,:) = [pt, 0];
else  % triangular
    pt = pt + preEncodeDur/2;
    ReadoutPulseWaveform(k,:) = [pt, -GPreEncode];
    k = k+1; pt = pt + preEncodeDur/2;
    ReadoutPulseWaveform(k,:) = [pt,0];
end
pt = pt + riseTime;
ReadoutStartTime = pt; % This is where the readout flat top begins

% If the slice refocus is not complete, we have to delay the readout
Delay = MinReadoutStartTime - ReadoutStartTime;
if Delay < 0 % we introduce a delay only if this occurs
    Delay = 0;
end

k = k+1;
ReadoutPulseWaveform(k,:) = [pt, GradMax];
te = pt + flatTop/2 - RFCenter;
pt = pt + flatTop;
k = k+1;
ReadoutPulseWaveform(k,:) = [pt, GradMax];
pt = pt + riseTime;
k = k+1;
ReadoutPulseWaveform(k,:) = [pt, 0];

min_slice_time = pt

fprintf('\nWith this tr of %0.2f msec you could acquire %d slices using these parameters.\n', tr, fix(tr/min_slice_time));
fprintf('However, in practice there are additional delays. A more likely number is %d slices.\n', fix(tr/(10 + min_slice_time)));

% For plotting purposes we have to keep track of the length of the waveform
% description vectors
lenRead = size(ReadoutPulseWaveform, 1);
lenRF   = size(RFGradientPulseWaveform, 1);

RFGradientPulseWaveform(lenRF+1,:) = [ReadoutPulseWaveform(lenRead)+Delay, 0];

plot(RFGradientPulseWaveform(1:lenRF+1,1), RFGradientPulseWaveform(1:lenRF+1,2), 'linewidth',2);
offset = 2.5 * GradMax; % for plotting
ylim([-GradMax - offset*1.1, GradMax*1.3]);
set(gcf,'color','white');
hold all;

ReadoutPulseWaveform(1,1) = -Delay;
plot(ReadoutPulseWaveform(1:lenRead,1)+Delay, ReadoutPulseWaveform(1:lenRead,2)-offset,'linewidth',2);
ylabel('Gradient Amplitude','FontSize', 16);
xlabel('Time (msec)','FontSize', 16);
text(RFCenter+Delay,-1.2 * GradMax,'\uparrow Center of RF pulse','FontSize', 16);
text(te+RFCenter+Delay,-1.2 * GradMax,'\downarrow Minimum te','FontSize', 16);
legend(' Slice select',' Readout');
min_slice_time =

    9.2596


With this tr of 2000.00 msec you could acquire 215 slices using these parameters.
However, in practice there are additional delays. A more likely number is 103 slices.

Image quality

fprintf('\nIMAGE QUALITY FACTORS\n')
fprintf('\tIn-plane resolution:............... %0.3f mm\n', Resolution*1000);
fprintf('\tFat-Water shift at %0.1f Tesla:...... %0.f Hz\n',FieldStrength, ChemShift);
fprintf('\tThis is a spatial shift of:........ %0.3f pixels\n', ChemShift/BandWidth);
IMAGE QUALITY FACTORS
	In-plane resolution:............... 0.750 mm
	Fat-Water shift at 3.0 Tesla:...... 294 Hz
	This is a spatial shift of:........ 0.200 pixels

The wider the bandwidth the lower the SNR. This is because noise appears at all frequencies. If you limit the pickup to signal only in a narrow frequency range, you will also limit the amount of noise that is brought in. This is a relative calculation balancing chemical shift and bandwidth.

SNRLoss = 100*(1 - sqrt(ChemShift/BandWidth));
fprintf('\tSNR loss from bandwidth: %0.3f %%\n', SNRLoss);
fprintf('\nMinimum te = %0.1f msec\n\n', te+Delay);
	SNR loss from bandwidth: 55.279 %

Minimum te = 5.3 msec