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