78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 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
5seq=mr.Sequence(); % Create a new sequence object
6fov=250e-3; Nx=64; Ny=64; % Define FOV and resolution
7thickness=3e-3; % slice thinckness
8Nslices=3;
9TE=40e-3;
11pe_enable=1; % a flag to quickly disable phase encoding (1/0) as needed for the delay calibration
12ro_os=1; % oversampling factor (in contrast to the product sequence we don't really need it)
13readoutTime=4.2e-4; % this controls the readout bandwidth
14partFourierFactor=0.75; % partial Fourier factor: 1: full sampling 0: start with ky=0
16tRFex=2e-3;
17tRFref=2e-3;
18spoilFactor=1.5; % spoiling gradient around the pi-pulse
20% Set system limits
21lims = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
22 'MaxSlew',130,'SlewUnit','T/m/s',...
23 'rfRingdownTime', 30e-6, 'rfDeadtime', 100e-6, 'adcDeadTime', 10e-6);
25% Create fat-sat pulse
26B0=2.89; % 1.5 2.89 3.0
27sat_ppm=-3.45;
28sat_freq=sat_ppm*1e-6*B0*lims.gamma;
29rf_fs = mr.makeGaussPulse(110*pi/180,'system',lims,'Duration',8e-3,...
30 'bandwidth',abs(sat_freq),'freqOffset',sat_freq, 'use', 'saturation');
31rf_fs.phaseOffset=-2*pi*rf_fs.freqOffset*mr.calcRfCenter(rf_fs); % compensate for the frequency-offset induced phase
32gz_fs = mr.makeTrapezoid('z',lims,'delay',mr.calcDuration(rf_fs),'Area',1/1e-4); % spoil up to 0.1mm
34% Create 90 degree slice selection pulse and gradient
35[rf, gz, gzReph] = mr.makeSincPulse(pi/2,'system',lims,'Duration',tRFex,...
36 'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4, 'use','excitation');
38% Create 90 degree slice refocusing pulse and gradients
39[rf180, gz180] = mr.makeSincPulse(pi,'system',lims,'Duration',tRFref,...
40 'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4,'PhaseOffset',pi/2,'use','refocusing');
41% [~, gzr_t, gzr_a]=mr.makeExtendedTrapezoidArea('z',gz180.amplitude,0,-gzReph.area+0.5*gz180.amplitude*gz180.fallTime,lims);
42% gz180n=mr.makeExtendedTrapezoid('z','system',lims,'times',[0 gz180.riseTime gz180.riseTime+gz180.flatTime+gzr_t]+gz180.delay, 'amplitudes', [0 gz180.amplitude gzr_a]);
43[~, gzr1_t, gzr1_a]=mr.makeExtendedTrapezoidArea('z',0,gz180.amplitude,spoilFactor*gz.area,lims);
44[~, gzr2_t, gzr2_a]=mr.makeExtendedTrapezoidArea('z',gz180.amplitude,0,-gzReph.area+spoilFactor*gz.area,lims);
45if gz180.delay>(gzr1_t(4)-gz180.riseTime)
46 gz180.delay=gz180.delay-(gzr1_t(4)-gz180.riseTime);
47else
48 rf180.delay=rf180.delay+(gzr1_t(4)-gz180.riseTime)-gz180.delay;
49 gz180.delay=0;
50end
51gz180n=mr.makeExtendedTrapezoid('z','system',lims,'times',[gzr1_t gzr1_t(4)+gz180.flatTime+gzr2_t]+gz180.delay, 'amplitudes', [gzr1_a gzr2_a]);
54% define the output trigger to play out with every slice excitatuion
55trig=mr.makeDigitalOutputPulse('osc0','duration', 100e-6); % possible channels: 'osc0','osc1','ext1'
57% Define other gradients and ADC events
58deltak=1/fov;
59kWidth = Nx*deltak;
61% Phase blip in shortest possible time
62blip_dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6/2)*10e-6*2; % we round-up the duration to 2x the gradient raster time
63% the split code below fails if this really makes a trpezoid instead of a triangle...
64gy = mr.makeTrapezoid('y',lims,'Area',-deltak,'Duration',blip_dur); % we use negative blips to save one k-space line on our way towards the k-space center
65%gy = mr.makeTrapezoid('y',lims,'amplitude',deltak/blip_dur*2,'riseTime',blip_dur/2, 'flatTime', 0);
67% readout gradient is a truncated trapezoid with dead times at the beginnig
68% and at the end each equal to a half of blip_dur
69% the area between the blips should be defined by kWidth
70% we do a two-step calculation: we first increase the area assuming maximum
71% slewrate and then scale down the amlitude to fix the area
72extra_area=blip_dur/2*blip_dur/2*lims.maxSlew; % check unit!;
73gx = mr.makeTrapezoid('x',lims,'Area',kWidth+extra_area,'duration',readoutTime+blip_dur);
74actual_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;
75gx.amplitude=gx.amplitude/actual_area*kWidth;
76gx.area = gx.amplitude*(gx.flatTime + gx.riseTime/2 + gx.fallTime/2);
77gx.flatArea = gx.amplitude*gx.flatTime;
79% calculate ADC
80% we use ramp sampling, so we have to calculate the dwell time and the
81% number of samples, which are will be qite different from Nx and
82% readoutTime/Nx, respectively.
83adcDwellNyquist=deltak/gx.amplitude/ro_os;
84% round-down dwell time to 100 ns
85adcDwell=floor(adcDwellNyquist*1e7)*1e-7;
86adcSamples=floor(readoutTime/adcDwell/4)*4; % on Siemens the number of ADC samples need to be divisible by 4
87% MZ: no idea, whether ceil,round or floor is better for the adcSamples...
88adc = mr.makeAdc(adcSamples,'Dwell',adcDwell,'Delay',blip_dur/2);
89% realign the ADC with respect to the gradient
90time_to_center=adc.dwell*((adcSamples-1)/2+0.5); % I've been told that Siemens samples in the center of the dwell period
91adc.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
92% this rounding actually makes the sampling points on odd and even readouts
93% to appear misalligned. However, on the real hardware this misalignment is
94% much stronger anyways due to the grdient delays
96% FOV positioning requires alignment to grad. raster... -> TODO
98% split the blip into two halves and produnce a combined synthetic gradient
99gy_parts = mr.splitGradientAt(gy, blip_dur/2, lims);
100[gy_blipup, gy_blipdown, ~]=mr.align('right',gy_parts(1),'left',gy_parts(2),gx);
101gy_blipdownup=mr.addGradients({gy_blipdown, gy_blipup}, lims);
103% pe_enable support
104gy_blipup.waveform=gy_blipup.waveform*pe_enable;
105gy_blipdown.waveform=gy_blipdown.waveform*pe_enable;
106gy_blipdownup.waveform=gy_blipdownup.waveform*pe_enable;
108% phase encoding and partial Fourier
110Ny_pre=round(partFourierFactor*Ny/2-1); % PE steps prior to ky=0, excluding the central line
111Ny_post=round(Ny/2+1); % PE lines after the k-space center including the central line
112Ny_meas=Ny_pre+Ny_post;
114% Pre-phasing gradients
115gxPre = mr.makeTrapezoid('x',lims,'Area',-gx.area/2);
116gyPre = mr.makeTrapezoid('y',lims,'Area',Ny_pre*deltak);
117[gxPre,gyPre]=mr.align('right',gxPre,'left',gyPre);
118% relax the PE prepahser to reduce stimulation
119gyPre = mr.makeTrapezoid('y',lims,'Area',gyPre.area,'Duration',mr.calcDuration(gxPre,gyPre));
120gyPre.amplitude=gyPre.amplitude*pe_enable;
122% Calculate delay times
123durationToCenter = (Ny_pre+0.5)*mr.calcDuration(gx);
124rfCenterInclDelay=rf.delay + mr.calcRfCenter(rf);
125rf180centerInclDelay=rf180.delay + mr.calcRfCenter(rf180);
126delayTE1=ceil((TE/2 - mr.calcDuration(rf,gz) + rfCenterInclDelay - rf180centerInclDelay)/lims.gradRasterTime)*lims.gradRasterTime;
127delayTE2=ceil((TE/2 - mr.calcDuration(rf180,gz180n) + rf180centerInclDelay - durationToCenter)/lims.gradRasterTime)*lims.gradRasterTime;
128assert(delayTE1>=0);
129%assert(delayTE2>=0);
130% now we merge slice refocusing, TE delay and pre-phasers into a single
131% block
132delayTE2=delayTE2+mr.calcDuration(rf180,gz180n);
133gxPre.delay=0;
134gxPre.delay=delayTE2-mr.calcDuration(gxPre);
135assert(gxPre.delay>=mr.calcDuration(rf180)); % gxPre may not overlap with the RF
136gyPre.delay=mr.calcDuration(rf180);
137assert(mr.calcDuration(gyPre)<=mr.calcDuration(gxPre)); % gyPre may not shift the timing
139% Define sequence blocks
141%seq.addBlock(mr.makeDelay(1)); % older scanners like Trio may need this
142 % dummy delay to keep up with timing
144for s=1:Nslices
145 seq.addBlock(rf_fs,gz_fs);
146 rf.freqOffset=gz.amplitude*thickness*(s-1-(Nslices-1)/2);
147 rf.phaseOffset=-2*pi*rf.freqOffset*mr.calcRfCenter(rf); % compensate for the slice-offset induced phase
148 rf180.freqOffset=gz180.amplitude*thickness*(s-1-(Nslices-1)/2);
149 rf180.phaseOffset=pi/2-2*pi*rf180.freqOffset*mr.calcRfCenter(rf180); % compensate for the slice-offset induced phase
150 seq.addBlock(rf,gz,trig);
151 seq.addBlock(mr.makeDelay(delayTE1));
152 seq.addBlock(rf180,gz180n,mr.makeDelay(delayTE2),gxPre,gyPre);
153 for i=1:Ny_meas
154 if i==1
155 seq.addBlock(gx,gy_blipup,adc); % Read the first line of k-space with a single half-blip at the end
156 elseif i==Ny_meas
157 seq.addBlock(gx,gy_blipdown,adc); % Read the last line of k-space with a single half-blip at the beginning
158 else
159 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
160 end
161 gx.amplitude = -gx.amplitude; % Reverse polarity of read gradient
162 end
163end
165%% check whether the timing of the sequence is correct
166[ok, error_report]=seq.checkTiming;
168if (ok)
169 fprintf('Timing check passed successfully\n');
170else
171 fprintf('Timing check failed! Error listing follows:\n');
172 fprintf([error_report{:}]);
173 fprintf('\n');
174end
176%% do some visualizations
179% trajectory calculation
181% plot k-spaces
183%% prepare the sequence output for the scanner
184seq.setDefinition('FOV', [fov fov thickness]);
185seq.setDefinition('Name', 'epi');
187seq.write('epi_se_rs.seq');