1% this is an experimental high-performance EPI sequence
2% which uses split gradients to overlap blips with the readout
3% gradients combined with ramp-samping
5% Set system limits
6sys = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
7 'MaxSlew',130,'SlewUnit','T/m/s',...
8 'rfRingdownTime', 30e-6, 'rfDeadtime', 100e-6,...
9 'adcDeadTime', 10e-6, 'B0', 2.89, ... % this is Siemens' 3T
10 'setAsDefault', true ... % new way of handing over the system parameters implicitly
11);
13seq=mr.Sequence(sys); % Create a new sequence object
14fov=256e-3; Nx=64; Ny=Nx; % Define FOV and resolution
15thickness=4e-3; % slice thinckness in mm
16sliceGap=1e-3; % slice gap im mm
17Nslices=1;
19pe_enable=1; % a flag to quickly disable phase encoding (1/0) as needed for the delay calibration
20ro_os=1; % oversampling factor (in contrast to the product sequence we don't really need it)
21readoutTime=4.2e-4; % this controls the readout bandwidth
22partFourierFactor=1; % partial Fourier factor: 1: full sampling 0: start with ky=0
24% Create fat-sat pulse
25sat_ppm=-3.45;
26rf_fs = mr.makeGaussPulse(110*pi/180,'system',sys,'Duration',8e-3,...
27 'bandwidth',abs(sat_ppm*1e-6*sys.B0*sys.gamma),'freqPPM',sat_ppm,'use','saturation');
28rf_fs.phasePPM=-2*pi*rf_fs.freqPPM*rf_fs.center; % compensate for the frequency-offset induced phase
29rf_fs.name='fat-sat'; % useful for debugging, can be seen in seq.plot
30gz_fs = mr.makeTrapezoid('z','delay',mr.calcDuration(rf_fs),'Area',0.1/1e-4); % spoil up to 0.1mm
31% Create 90 degree slice selection pulse and gradient
32[rf, gz, gzReph] = mr.makeSincPulse(pi/2,'Duration',2e-3,...
33 'SliceThickness',thickness,'apodization',0.42,'timeBwProduct',4,'use','excitation');
34rf.name='rf90'; % useful for debugging, can be seen in seq.plot
36% define the output trigger to play out with every slice excitatuion
37trig=mr.makeDigitalOutputPulse('osc0','duration', 100e-6); % possible channels: 'osc0','osc1','ext1'
39% Define other gradients and ADC events
40deltak=1/fov;
41kWidth = Nx*deltak;
43% Phase blip in shortest possible time
44blip_dur = ceil(2*sqrt(deltak/sys.maxSlew)/10e-6/2)*10e-6*2; % we round-up the duration to 2x the gradient raster time
45% the split code below fails if this really makes a trpezoid instead of a triangle...
46gy = mr.makeTrapezoid('y','Area',-deltak,'Duration',blip_dur); % we use negative blips to save one k-space line on our way towards the k-space center
47%gy = mr.makeTrapezoid('y',lims,'amplitude',deltak/blip_dur*2,'riseTime',blip_dur/2, 'flatTime', 0);
49% readout gradient is a truncated trapezoid with dead times at the beginnig
50% and at the end each equal to a half of blip_dur
51% the area between the blips should be defined by kWidth
52% we do a two-step calculation: we first increase the area assuming maximum
53% slewrate and then scale down the amlitude to fix the area
54extra_area=blip_dur/2*blip_dur/2*sys.maxSlew; % check unit!;
55gx = mr.makeTrapezoid('x','Area',kWidth+extra_area,'duration',readoutTime+blip_dur);
56actual_area=gx.area-gx.amplitude/gx.riseTime*blip_dur/2*blip_dur/2/2-gx.amplitude/gx.fallTime*blip_dur/2*blip_dur/2/2;
57gx.amplitude=gx.amplitude/actual_area*kWidth;
58gx.area = gx.amplitude*(gx.flatTime + gx.riseTime/2 + gx.fallTime/2);
59gx.flatArea = gx.amplitude*gx.flatTime;
60gx.name='Gro'; % useful for debugging, can be seen in seq.plot
62% calculate ADC
63% we use ramp sampling, so we have to calculate the dwell time and the
64% number of samples, which are will be qite different from Nx and
65% readoutTime/Nx, respectively.
66adcDwellNyquist=deltak/gx.amplitude/ro_os;
67% round-down dwell time to 100 ns
68adcDwell=floor(adcDwellNyquist*1e7)*1e-7;
69adcSamples=floor(readoutTime/adcDwell/4)*4; % on Siemens the number of ADC samples need to be divisible by 4
70% MZ: no idea, whether ceil,round or floor is better for the adcSamples...
71adc = mr.makeAdc(adcSamples,'Dwell',adcDwell,'Delay',blip_dur/2);
72% realign the ADC with respect to the gradient
73time_to_center=adc.dwell*((adcSamples-1)/2+0.5); % I've been told that Siemens samples in the center of the dwell period
74adc.delay=round((gx.riseTime+gx.flatTime/2-time_to_center)*1e6)*1e-6; % we adjust the delay to align the trajectory with the gradient. We have to aligh the delay to 1us
75% this rounding actually makes the sampling points on odd and even readouts
76% to appear misalligned. However, on the real hardware this misalignment is
77% much stronger anyways due to the grdient delays
79% FOV positioning requires alignment to grad. raster... -> TODO
81% split the blip into two halves and produce a combined synthetic gradient
82gy_parts = mr.splitGradientAt(gy, blip_dur/2, sys);
83[gy_blipup, gy_blipdown,~]=mr.align('right',gy_parts(1),'left',gy_parts(2),gx);
84gy_blipdownup=mr.addGradients({gy_blipdown, gy_blipup}, sys);
86% pe_enable support
87gy_blipup.waveform=gy_blipup.waveform*pe_enable;
88gy_blipdown.waveform=gy_blipdown.waveform*pe_enable;
89gy_blipdownup.waveform=gy_blipdownup.waveform*pe_enable;
91% phase encoding and partial Fourier
93Ny_pre=round(partFourierFactor*Ny/2-1); % PE steps prior to ky=0, excluding the central line
94Ny_post=round(Ny/2+1); % PE lines after the k-space center including the central line
95Ny_meas=Ny_pre+Ny_post;
97% Pre-phasing gradients
98gxPre = mr.makeTrapezoid('x','Area',-gx.area/2);
99gyPre = mr.makeTrapezoid('y','Area',Ny_pre*deltak);
100[gxPre,gyPre,gzReph]=mr.align('right',gxPre,'left',gyPre,gzReph);
101% relax the PE prepahser to reduce stimulation
102gyPre = mr.makeTrapezoid('y','Area',gyPre.area,'Duration',mr.calcDuration(gxPre,gyPre,gzReph));
103gyPre.amplitude=gyPre.amplitude*pe_enable;
105% slice positions
106slicePositions=(thickness+sliceGap)*((0:(Nslices-1)) - (Nslices-1)/2);
107slicePositions=slicePositions([1:2:Nslices 2:2:Nslices]); % reorder slices for an interleaved acquisition (optional)
109% Define sequence blocks
110%seq.addBlock(mr.makeDelay(1)); % older scanners like Trio may need this
111 % dummy delay to keep up with timing
112for s=1:Nslices
113 seq.addBlock(rf_fs,gz_fs);
114 rf.freqOffset=gz.amplitude*slicePositions(s);
115 rf.phaseOffset=-2*pi*rf.freqOffset*mr.calcRfCenter(rf); % compensate for the slice-offset induced phase
116 seq.addBlock(rf,gz,trig);
117 seq.addBlock(gxPre,gyPre,gzReph);
118 for i=1:Ny_meas
119 if i==1
120 seq.addBlock(gx,gy_blipup,adc); % Read the first line of k-space with a single half-blip at the end
121 elseif i==Ny_meas
122 seq.addBlock(gx,gy_blipdown,adc); % Read the last line of k-space with a single half-blip at the beginning
123 else
124 seq.addBlock(gx,gy_blipdownup,adc); % Read an intermediate line of k-space with a half-blip at the beginning and a half-blip at the end
125 end
126 gx.amplitude = -gx.amplitude; % Reverse polarity of read gradient
127 end
128end
130%% check whether the timing of the sequence is correct
131[ok, error_report]=seq.checkTiming;
133if (ok)
134 fprintf('Timing check passed successfully\n');
135else
136 fprintf('Timing check failed! Error listing follows:\n');
137 fprintf([error_report{:}]);
138 fprintf('\n');
139end
141%% do some visualizations
145rf.freqOffset=0;
146rf.phaseOffset=0;
147[rf_bw,rf_f0,rf_spectrum,rf_w]=mr.calcRfBandwidth(rf);
148xlim(3*[-rf_bw rf_bw]);
150%% trajectory calculation
152% plot k-spaces
155%axis off;
157%% prepare the sequence output for the scanner
158seq.setDefinition('Name', 'epi');
159seq.setDefinition('FOV', [fov fov max(slicePositions)-min(slicePositions)+thickness]);
160seq.setDefinition('ReceiverGainHigh',1);
161% the following definitions only have effect in conjunction with LABELs
162%seq.setDefinition('SlicePositions', slicePositions);
163%seq.setDefinition('SliceThickness', thickness);
164%seq.setDefinition('SliceGap', sliceGap);
166seq.write('epi_rs.seq');