78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1function [rf, gz, gzr, delay] = makeAdiabaticPulse(type,varargin)
2%makeAdiabaticPulse make an aiabatic inversion pulse
3% a wrapper to a python function(see below). See supported params below
4% in the 'parser' section. Currently it will probably only work on
5% Linux. On my system I could install the required Python library by
6% executing "pip3 install sigpy"
7% type must be one of {'hypsec','wurst'}
8% BE CAREFUL, some parameters only affect certain pulse types and are
9% ignored for other; e.g. bandwidth is ignored if type='hypsec'.
10%
11% hypsec(n=512, beta=800, mu=4.9, dur=0.012)
12% Design a hyperbolic secant adiabatic pulse.
13%
14% mu * beta becomes the amplitude of the frequency sweep
15%
16% Args:
17% n (int): number of samples (should be a multiple of 4).
18% beta (float): AM waveform parameter.
19% mu (float): a constant, determines amplitude of frequency sweep.
20% dur (float): pulse time (s).
21%
22% Returns:
23% 2-element tuple containing
24%
25% - **a** (*array*): AM waveform.
26% - **om** (*array*): FM waveform (radians/s).
27%
28% References:
29% Baum, J., Tycko, R. and Pines, A. (1985). 'Broadband and adiabatic
30% inversion of a two-level system by phase-modulated pulses'.
31% Phys. Rev. A., 32:3435-3447.
32%
33% wurst(n=512, n_fac=40, bw=40000.0, dur=0.002)
34% Design a WURST (wideband, uniform rate, smooth truncation) adiabatic
35% inversion pulse
36%
37% Args:
38% n (int): number of samples (should be a multiple of 4).
39% n_fac (int): power to exponentiate to within AM term. ~20 or greater is
40% typical.
41% bw (float): pulse bandwidth.
42% dur (float): pulse time (s).
43%
44%
45% Returns:
46% 2-element tuple containing
47% - **a** (*array*): AM waveform.
48% - **om** (*array*): FM waveform (radians/s).
49%
50% References:
51% Kupce, E. and Freeman, R. (1995). 'Stretched Adiabatic Pulses for
52% Broadband Spin Inversion'.
53% J. Magn. Reson. Ser. A., 117:246-256.
56validPulseTypes = {'hypsec','wurst'};
57validPulseUses = mr.getSupportedRfUse();
59persistent parser
60if isempty(parser)
61 parser = mr.aux.InputParserCompat;
62 parser.FunctionName = 'makeAdiabaticPulse';
64 % RF params
65 addRequired(parser, 'type', @(x) any(validatestring(x,validPulseTypes)));
66 addOptional(parser, 'system', [], @isstruct);
67 addParamValue(parser, 'duration', 10e-3, @isnumeric);
68 addParamValue(parser, 'freqOffset', 0, @isnumeric);
69 addParamValue(parser, 'phaseOffset', 0, @isnumeric);
70 addParamValue(parser, 'freqPPM', 0, @isnumeric);
71 addParamValue(parser, 'phasePPM', 0, @isnumeric);
72 addParamValue(parser, 'beta', 800, @isnumeric);
73 addParamValue(parser, 'mu', 4.9, @isnumeric);
74 addParamValue(parser, 'n_fac', 40, @isnumeric);
75 addParamValue(parser, 'bandwidth', 40000, @isnumeric);
76 addParamValue(parser, 'adiabaticity', 4, @isnumeric);
77 % Slice params
78 addParamValue(parser, 'maxGrad', 0, @isnumeric);
79 addParamValue(parser, 'maxSlew', 0, @isnumeric);
80 addParamValue(parser, 'sliceThickness', 0, @isnumeric);
81 addParamValue(parser, 'delay', 0, @isnumeric);
82 addParamValue(parser, 'dwell', 0, @isnumeric); % dummy default value
83 % whether it is a refocusing pulse (for k-space calculation)
84 addParamValue(parser, 'use', 'u', @(x) any(validatestring(x,validPulseUses)));
85 % optional Python command
86 addParamValue(parser, 'pythonCmd', '', @(x)isstring(x)||ischar(x));
87end
89parse(parser, type, varargin{:});
90opt = parser.Results;
92if isempty(opt.system)
93 sys=mr.opts();
94else
95 sys=opt.system;
96end
98if opt.dwell==0
99 opt.dwell=sys.rfRasterTime;
100end
102% find/check python
103if ~isempty(opt.pythonCmd)
104 [status, result]=system([opt.pythonCmd ' --version']);
105 if status~=0
106 error(['provided python executable ''' opt.pythonCmd ''' returns an error on the version check']);
107 end
108 if ispc
109 [status, result] = system(sprintf('%s -c "import sigpy" 2>nul',opt.pythonCmd));
110 else
111 [status, result] = system(sprintf('%s -c "import sigpy" 2>/dev/null',opt.pythonCmd));
112 end
113 if status~=0
114 error(['provided python executable ''' opt.pythonCmd ''' returns an error on the sigPy check']);
115 end
116 python=opt.pythonCmd;
117else
118 [avail, python]=mr.aux.isSigPyAvailable();
119 if ~avail
120 error('python executable with installed sigPy not found, please check your system PATH settings and Python installation');
121 end
122end
123% add quotes in case Python install path contains spaces or alike characters
124if python(1)~='"'
125 python=['"' python '"'];
126end
128Nraw = round(opt.duration/opt.dwell+eps);
129N = floor(Nraw/4)*4; % number of points must be divisible by four -- this is a requirement of the underlying library
131if ispc()
132 switch type
133 case 'hypsec'
134 cmd=[python ' -c "import sigpy.mri.rf;pulse=sigpy.mri.rf.hypsec(' ... % hypsec(n=512, beta=800, mu=4.9, dur=0.012)
135 'n=' num2str(N) ',beta=' num2str(opt.beta) ',' ...
136 'mu=' num2str(opt.mu) ',dur=' num2str(opt.duration) ...
137 ');print(*pulse[0]);print(*pulse[1])"'];
138 case 'wurst'
139 cmd=[python ' -c "import sigpy.mri.rf;pulse=sigpy.mri.rf.wurst(' ... % wurst(n=512, n_fac=40, bw=40000.0, dur=0.002)
140 'n=' num2str(N) ',n_fac=' num2str(opt.n_fac) ',' ...
141 'bw=' num2str(opt.bandwidth) ',dur=' num2str(opt.duration) ...
142 ');print(*pulse[0]);print(*pulse[1])"'];
143 otherwise
144 error('unsupported adiabatic pulse type');
145 end
146else
147 switch type
148 case 'hypsec'
149 cmd=[python ' -c $''import sigpy.mri.rf\npulse=sigpy.mri.rf.hypsec(' ... % hypsec(n=512, beta=800, mu=4.9, dur=0.012)
150 'n=' num2str(N) ',beta=' num2str(opt.beta) ',' ...
151 'mu=' num2str(opt.mu) ',dur=' num2str(opt.duration) ...
152 ')\nprint(*pulse[0])\nprint(*pulse[1])'''];
153 case 'wurst'
154 cmd=[python ' -c $''import sigpy.mri.rf\npulse=sigpy.mri.rf.wurst(' ... % wurst(n=512, n_fac=40, bw=40000.0, dur=0.002)
155 'n=' num2str(N) ',n_fac=' num2str(opt.n_fac) ',' ...
156 'bw=' num2str(opt.bandwidth) ',dur=' num2str(opt.duration) ...
157 ')\nprint(*pulse[0])\nprint(*pulse[1])'''];
158 otherwise
159 error('unsupported adiabatic pulse type');
160 end
161end
162%fprintf('cmd=%s\n',cmd);
163[status, result]=system(cmd);
165if status~=0
166 error('executing python command failed');
167end
169lines = regexp(result,'\n','split'); % the response from the python call contains some garbage
170% look for two usable result vectors
171for i=1:length(lines)-1
172 try
173 am=str2num(lines{i});
174 fm=str2num(lines{i+1});
175 if length(am)==N && length(fm)==N
176 break;
177 end
178 catch
179 continue;
180 end
181end
182if length(am)~=N || length(fm)~=N
183 error('could not find usable data in the response of the Python command');
184end
186pm=cumsum(fm)*opt.dwell;
188[dfm,ifm]=min(abs(fm)); % find the center of the pulse
189% we will also use the ocasion to find the rate of change of the frequency
190% at the center of the pulse
191if dfm==0
192 pm0=pm(ifm);
193 am0=am(ifm);
194 roc_fm0=abs(fm(ifm+1)-fm(ifm-1))/2/opt.dwell;
195else
196 % we need to bracket the zero-crossing
197 if fm(ifm)*fm(ifm+1) < 0
198 b=1;
199 else
200 b=-1;
201 end
202 pm0=(pm(ifm)*fm(ifm+b)-pm(ifm+b)*fm(ifm))/(fm(ifm+b)-fm(ifm));
203 am0=(am(ifm)*fm(ifm+b)-am(ifm+b)*fm(ifm))/(fm(ifm+b)-fm(ifm));
204 roc_fm0=abs(fm(ifm)-fm(ifm+b))/opt.dwell;
205end
206pm=pm-pm0;
207a=(roc_fm0*opt.adiabaticity)^0.5/2/pi/am0;
209signal = a*am.*exp(1i*pm);
211if (N~=Nraw)
212 % we need to pad the signal vector
213 Npad=Nraw-N;
214 signal=[zeros(1,Npad-floor(Npad/2)) signal zeros(1,floor(Npad/2))];
215 N=Nraw;
216end
218%BW = opt.timeBwProduct/opt.duration;
219t = ((1:N)-0.5)*opt.dwell;
220%flip = abs(sum(signal))*opt.dwell*2*pi;
222rf.type = 'rf';
223rf.signal = signal;
224rf.t = t;
225rf.shape_dur=N*opt.dwell;
226rf.freqOffset = opt.freqOffset;
227rf.phaseOffset = opt.phaseOffset;
228rf.freqPPM = opt.freqPPM;
229rf.phasePPM = opt.phasePPM;
230rf.deadTime = sys.rfDeadTime;
231rf.ringdownTime = sys.rfRingdownTime;
232rf.delay = opt.delay;
233rf.center = mr.calcRfCenter(rf);
234if ~isempty(opt.use)
235 rf.use=opt.use;
236else
237 rf.use='inversion';
238end
239if rf.deadTime > rf.delay
240 rf.delay = rf.deadTime;
241end
243if nargout > 1
244 assert(opt.sliceThickness > 0,'SliceThickness must be provided');
245 if opt.maxGrad > 0
246 sys.maxGrad = opt.maxGrad;
247 end
248 if opt.maxSlew > 0
249 sys.maxSlew = opt.maxSlew;
250 end
252 switch type
253 case 'hypsec'
254 BW=mr.calcRfBandwidth(rf,0.1);
255 case 'wurst'
256 BW=opt.bandwidth;
257 otherwise
258 error('unsupported pulse type')
259 end
261 amplitude = BW/opt.sliceThickness;
262 area = amplitude*opt.duration;
263 gz = mr.makeTrapezoid('z', sys, 'flatTime', opt.duration, ...
264 'flatArea', area);
265 gzr= mr.makeTrapezoid('z', sys, 'Area', -area*(1-rf.center/rf.shape_dur)-0.5*(gz.area-area));
266 if rf.delay > gz.riseTime
267 gz.delay = ceil((rf.delay - gz.riseTime)/sys.gradRasterTime)*sys.gradRasterTime; % round-up to gradient raster
268 end
269 if rf.delay < (gz.riseTime+gz.delay)
270 rf.delay = gz.riseTime+gz.delay; % these are on the grad raster already which is coarser
271 end
272end
274% v1.4 finally eliminates RF zerofilling
275% if rf.ringdownTime > 0
276% tFill = (1:round(rf.ringdownTime/1e-6))*1e-6; % Round to microsecond
277% rf.t = [rf.t rf.t(end)+tFill];
278% rf.signal = [rf.signal, zeros(size(tFill))];
279% end
280if nargout > 3
281 delay=mr.makeDelay(mr.calcDuration(rf)); % calcDuration already includes the ringdown time
282end
284% RF amplitude check
285rf_amplitude=max(abs(rf.signal));
286if rf_amplitude>sys.maxB1
287 warning('WARNING: system maximum RF amplitude exceeded (%.01f%%)', rf_amplitude/sys.maxB1*100);
288end
290end