Contents

% we define here a really crude True-FISP a.k.a. bSSFP sequence
% there is no user control for TR/TE, you just specify the ADC time and the
% RF parameters and the rest is calculated to find the fastest posible
% timing. The sequence intensively uses the extended trapezoid
% functionality to achieve near-optimal timing. Due to the requirement for
% splitting the sequence into blocks the TR is increased by approximately
% 40-60 us (rfDeadTime+adcDeadTime) in comparison to the true minimum TR

% set system limits
% had to slow down ramps and increase adc_duration to avoid stimulation
sys = mr.opts('MaxGrad',30,'GradUnit','mT/m',...
    'MaxSlew',140,'SlewUnit','T/m/s',...
    'rfRingdownTime', 20e-6, 'rfDeadTime', 100e-6, ...
    'adcDeadTime', 20e-6);

seq=mr.Sequence(sys);              % Create a new sequence object
fov=220e-3; Nx=256; Ny=256;     % Define FOV and resolution


% ADC duration (controls TR/TE)
adc_dur=2560; %us

% RF parameters
alpha=40; % deg
thick=4; %mm
rf_dur=600; % us
rf_apo=0.5;
rf_bwt=1.5;

% Create 'alpha' degree slice selection pulse and gradient
[rf, gz, gzReph] = mr.makeSincPulse(alpha*pi/180,'Duration',rf_dur*1e-6,...
    'SliceThickness',thick*1e-3,'apodization',rf_apo,'timeBwProduct',rf_bwt,'system',sys);

% Define other gradients and ADC events
deltak=1/fov;
gx = mr.makeTrapezoid('x','FlatArea',Nx*deltak,'FlatTime',adc_dur*1e-6,'system',sys);
adc = mr.makeAdc(Nx,'Duration',gx.flatTime,'Delay',gx.riseTime,'system',sys);
gxPre = mr.makeTrapezoid('x','Area',-gx.area/2,'system',sys);
phaseAreas = ((0:Ny-1)-Ny/2)*deltak;

% now we have to reshuffle gradients to achieve a half-way optimal timing
% new gz will consist of two parts:
% 1: slice refocusing from the previous TR followed by slice selection
%    including the plato an a small bit of the ramp-down
% 2: the remainder of the ramp-down and the slice refocusing for the next TR
gz_parts=mr.splitGradientAt(gz,mr.calcDuration(rf));
gz_parts(1).delay=mr.calcDuration(gzReph);
gz_1=mr.addGradients({gzReph,gz_parts(1)},'system',sys);
[rf,~]=mr.align('right',rf,gz_1);
gz_parts(2).delay=0;
gzReph.delay=mr.calcDuration(gz_parts(2));
gz_2=mr.addGradients({gz_parts(2),gzReph},'system',sys);

% new gr will consist of two parts:
% 1: prephaser followed by a part of the read gradient including the
%    beginning of the ramp-down
% 2: the remainer of the ramp-down and the second "prephaser"
gx_parts=mr.splitGradientAt(gx,ceil(mr.calcDuration(adc)/sys.gradRasterTime)*sys.gradRasterTime);
gx_parts(1).delay=mr.calcDuration(gxPre);
gx_1=mr.addGradients({gxPre,gx_parts(1)},'system',sys);
adc.delay=adc.delay+mr.calcDuration(gxPre); % we cannot use mr.align here because the adc duration maz be not aligneed to the grad raster
gx_parts(2).delay=0;
gxPre.delay=mr.calcDuration(gx_parts(2));
gx_2=mr.addGradients({gx_parts(2),gxPre},'system',sys);

% Calculate timing
%pe_dur=min(max(tp2(end),tpr1(end-2)),max(tp(end-3),tpr2(end)));
gxPre.delay=0; % otherwise duration below is misreported
%pe_dur=max(mr.calcDuration(gz_2),mr.calcDuration(gxPre)+gx.riseTime);
pe_dur=mr.calcDuration(gx_2);

% adjust delays to align objects
%gz_1.delay=max(mr.calcDuration(gx_2)-rf.delay,0);
gz_1.delay=max(mr.calcDuration(gx_2)-rf.delay+rf.ringdownTime,0); % this rf.ringdownTime is needed to center the ADC and the gradient echo in the center of RF-RF period

rf.delay=rf.delay+gz_1.delay;

% finish timing calculation
TR=mr.calcDuration(gz_1)+mr.calcDuration(gx_1);
TE=TR/2;

% alpha / 2 preparation: shorter TR and half-angle, no PE, no RO
% create 0.5*alpha prep pulse
rf05=rf;
rf05.signal=0.5*rf.signal;
seq.addBlock(rf05,gz_1,mr.makeLabel('SET','ONCE', 1)); % we also label the few following blocks as preparing to exclude them if the sequence is repeated
% the following delay calculation fails for agressive sequence timing
prepDelay=mr.makeDelay(round((TR/2-mr.calcDuration(gz_1))/sys.gradRasterTime)*sys.gradRasterTime); % I know this round() is not 100% correct
gx_1_1=mr.makeExtendedTrapezoidArea('x',0,gx_2.first,-gx_2.area,sys);
gyPre_2 = mr.makeTrapezoid('y','Area',phaseAreas(end),'Duration',pe_dur,'system',sys); % last PE step (in case of repetitions)
seq.addBlock(mr.align('left',prepDelay,gz_2,gyPre_2,'right',gx_1_1));
seq.addBlock(mr.makeLabel('SET','ONCE', 0)); % remove preparing block label

% Loop over phase encodes and define sequence blocks
for i=1:Ny
    rf.phaseOffset=pi*mod(i,2);
    adc.phaseOffset=pi*mod(i,2);

    gyPre_1 = mr.scaleGrad(gyPre_2,-1); % undo previous PE step
    gyPre_2 = mr.makeTrapezoid('y','Area',phaseAreas(i),'Duration',pe_dur,'system',sys); % current PE step

    seq.addBlock(rf,gz_1, gyPre_1, gx_2);
    seq.addBlock(gx_1,gyPre_2, gz_2,adc);
end
% finish the x-grad shape
seq.addBlock(gx_2,mr.makeLabel('SET','ONCE', 2)); % we also label this block as the exit block, which excludes it from all but last repetitions if the sequence is repeated

% check that the calculated TR was reached
% alpha / 2 prep takes 4 blocks
assert(TR==(mr.calcDuration(seq.getBlock(5))+mr.calcDuration(seq.getBlock(6))));

fprintf('Sequence ready\n');
fprintf('TR=%03f ms  TE=%03f ms\n', TR*1e3, TE*1e3);
Sequence ready
TR=4.740000 ms  TE=2.370000 ms

check whether the timing of the sequence is correct

[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end
Timing check passed successfully

prepare export

seq.setDefinition('FOV', [fov fov thick*1e-3]);
seq.setDefinition('Name', 'trufi');

seq.write('trufi.seq')       % Write to pulseq file

%seq.install('siemens');

plots and displays

seq.plot('timeDisp','us');

plot entire interpolated wave forms -- good for debugging of gaps, jumps,

etc, but is relatively slow

%gw=seq.gradient_waveforms();
%figure; plot(gw'); % plot the entire gradient waveform
gw_data=seq.waveforms_and_times();
figure; plot(gw_data{1}(1,:),gw_data{1}(2,:),gw_data{2}(1,:),gw_data{2}(2,:),gw_data{3}(1,:),gw_data{3}(2,:)); % plot the entire gradient waveform
title('gradient waveforms');

trajectory calculation

[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

% plot k-spaces

figure; plot(t_ktraj,ktraj'); title('k-space components as functions of time'); % plot the entire k-space trajectory
figure; plot(ktraj(1,:),ktraj(2,:),'b',...
             ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % a 2D plot
title('2D k-space');

very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within slewrate limits

rep = seq.testReport;
fprintf([rep{:}]);
Number of blocks: 516
Number of events:
   RF:      257
   Gx:      514
   Gy:      513
   Gz:      514
   ADC:     256
Sequence duration: 1.216560s
TE: 0.002370s
TR: 0.004740s
Flip angle: 20.00°
Flip angle: 40.00°
Unique k-space positions (a.k.a. columns, rows, etc): 256
Unique k-space positions (a.k.a. columns, rows, etc): 256
Dimensions: 2
   Spatial resolution: 0.86 mm
   Spatial resolution: 0.86 mm
Repetitions/slices/contrasts: 1  range: [1 1]
   65536 k-space position(s) repeated 1 times
Cartesian encoding trajectory detected
Block timing check passed successfully
Max. Gradient: 1276596 Hz/m == 29.98 mT/m
Max. Gradient: 1003135 Hz/m == 23.56 mT/m
Max. Gradient: 1109375 Hz/m == 26.06 mT/m
Max. Slew Rate: 5.80271e+09 Hz/m/s == 136.29 T/m/s
Max. Slew Rate: 5.9596e+09 Hz/m/s == 139.98 T/m/s
Max. Slew Rate: 5.68182e+09 Hz/m/s == 133.45 T/m/s
Max. Absolute Gradient: 1966390 Hz/m == 46.19 mT/m
Max. Absolute Slew Rate: 1.00733e+10 Hz/m/s == 236.60 T/m/s