1function [rf, gz, gzr, delay] = makeSLRpulse(flip,varargin)
2%makeSLRpulse make an SLR 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%
8% sigpy.mri.rf.dzrf = dzrf(n=64, tb=4, ptype='st', ftype='ls', d1=0.01, d2=0.01, cancel_alpha_phs=False)
9% Primary function for design of pulses using the SLR algorithm.
10%
11% Args:
12% n (int): number of time points.
13% tb (int): pulse time bandwidth product.
14% ptype (string): pulse type, 'st' (small-tip excitation), 'ex' (pi/2
15% excitation pulse), 'se' (spin-echo pulse), 'inv' (inversion), or
16% 'sat' (pi/2 saturation pulse).
17% ftype (string): type of filter to use: 'ms' (sinc), 'pm'
18% (Parks-McClellan equal-ripple), 'min' (minphase using factored pm),
19% 'max' (maxphase using factored pm), 'ls' (least squares).
20% d1 (float): passband ripple level in :math:'M_0^{-1}'.
21% d2 (float): stopband ripple level in :math:'M_0^{-1}'.
22% filterType (str): filter type to use, e.g. sinc (ms),
23% least-squares (ls), etc. Refer to sigpy.rf documentation.
24%
25% Returns:
26% rf (array): designed RF pulse.
27%
28% References:
29% Pauly, J., Le Roux, Patrick., Nishimura, D., and Macovski, A.(1991).
30% Parameter Relations for the Shinnar-LeRoux Selective Excitation
31% Pulse Design Algorithm.
32% IEEE Transactions on Medical Imaging, Vol 10, No 1, 53-65.
33%
35validPulseUses = mr.getSupportedRfUse();
37persistent parser
38if isempty(parser)
39 parser = mr.aux.InputParserCompat;
40 parser.FunctionName = 'makeSLRpulse';
42 % RF params
43 addRequired(parser, 'flipAngle', @isnumeric);
44 addOptional(parser, 'system', [], @isstruct);
45 addParamValue(parser, 'duration', 1e-3, @isnumeric);
46 addParamValue(parser, 'freqOffset', 0, @isnumeric);
47 addParamValue(parser, 'phaseOffset', 0, @isnumeric);
48 addParamValue(parser, 'freqPPM', 0, @isnumeric);
49 addParamValue(parser, 'phasePPM', 0, @isnumeric);
50 addParamValue(parser, 'timeBwProduct', 4, @isnumeric);
51 addParamValue(parser, 'passbandRipple', 0.01, @isnumeric);
52 addParamValue(parser, 'stopbandRipple', 0.01, @isnumeric);
53 addParamValue(parser, 'filterType', 'ms', @isstr);
54 %addParamValue(parser, 'apodization', 0, @isnumeric);
55 %addParamValue(parser, 'centerpos', 0.5, @isnumeric);
56 % Slice params
57 addParamValue(parser, 'maxGrad', 0, @isnumeric);
58 addParamValue(parser, 'maxSlew', 0, @isnumeric);
59 addParamValue(parser, 'sliceThickness', 0, @isnumeric);
60 addParamValue(parser, 'delay', 0, @isnumeric);
61 addParamValue(parser, 'dwell', 0, @isnumeric); % dummy default value
62 % whether it is a refocusing pulse (for k-space calculation)
63 addParamValue(parser, 'use', 'excitation', @(x) any(validatestring(x,validPulseUses)));
64 addParamValue(parser, 'recenterOnSample', false, @islogical);
65 % optional Python command
66 addParamValue(parser, 'pythonCmd', '', @(x)isstring(x)||ischar(x));
67end
68parse(parser, flip, varargin{:});
69opt = parser.Results;
71if isempty(opt.system)
72 sys=mr.opts();
73else
74 sys=opt.system;
75end
77if opt.dwell==0
78 opt.dwell=sys.rfRasterTime;
79end
81% find/check python
82if ~isempty(opt.pythonCmd)
83 [status, result]=system([opt.pythonCmd ' --version']);
84 if status~=0
85 error(['provided python executable ''' opt.pythonCmd ''' returns an error on the version check']);
86 end
87 if ispc
88 [status, result] = system(sprintf('%s -c "import sigpy" 2>nul',opt.pythonCmd));
89 else
90 [status, result] = system(sprintf('%s -c "import sigpy" 2>/dev/null',opt.pythonCmd));
91 end
92 if status~=0
93 error(['provided python executable ''' opt.pythonCmd ''' returns an error on the sigPy check']);
94 end
95 python=opt.pythonCmd;
96else
97 [avail, python]=mr.aux.isSigPyAvailable();
98 if ~avail
99 error('python executable with installed sigPy not found, please check your system PATH settings and Python installation');
100 end
101end
102% add quotes in case Python install path contains spaces or alike characters
103if python(1)~='"'
104 python=['"' python '"'];
105end
107add_opt='';
109switch opt.use
110 case 'excitation'
111 if opt.flipAngle <= pi/6
112 ptype='st';
113 else
114 ptype='ex';
115 %add_opt=',cancel_alpha_phs=True';
116 end
117 case 'refocusing'
118 ptype='se';
119 case 'inversion'
120 ptype='inv';
121 case 'saturation'
122 ptype='sat';
123 otherwise
124 ptype='st';
125end
127N = round(opt.duration/opt.dwell);
128% on Windows it looks like the $ and '' are not needed and ; can be used in place of \n
129if ispc()
130 cmd=[python ' -c "import sigpy.mri.rf;pulse=sigpy.mri.rf.dzrf(' num2str(N) ...
131 ',' num2str(opt.timeBwProduct) ',ptype=''' ptype '''' ...
132 ',d1=' num2str(opt.passbandRipple) ',d2=' num2str(opt.stopbandRipple) ...
133 ',ftype=''' opt.filterType '''' add_opt ');print(*pulse)"'];
134else
135 cmd=[python ' -c $''import sigpy.mri.rf\npulse=sigpy.mri.rf.dzrf(' num2str(N) ...
136 ',' num2str(opt.timeBwProduct) ',ptype=\''' ptype '\''' ...
137 ',d1=' num2str(opt.passbandRipple) ',d2=' num2str(opt.stopbandRipple) ...
138 ',ftype=\''' opt.filterType '\''' add_opt ')\nprint(*pulse)'''];
139end
140%fprintf('cmd=%s\n',cmd);
141[status, result]=system(cmd);
143if status~=0
144 error('executing python command failed, error message was: %s', result);
145end
147lines = regexp(result,'\n','split'); % the response from the python call contains some garbage
148% look for a usable result vector
149for i=1:length(lines)
150 try
151 signal = str2num(lines{i});
152 if length(signal)==N
153 break;
154 end
155 catch
156 continue;
157 end
158end
159if length(signal)~=N
160 error('could not find usable data in the response of the Python command');
161end
163BW = opt.timeBwProduct/opt.duration;
164t = ((1:N)-0.5)*opt.dwell;
165flip = abs(sum(signal))*opt.dwell*2*pi;
166signal = signal*opt.flipAngle/flip;
168rf.type = 'rf';
169rf.signal = signal;
170rf.t = t;
171rf.shape_dur=N*opt.dwell;
172rf.freqOffset = opt.freqOffset;
173rf.phaseOffset = opt.phaseOffset;
174rf.freqPPM = opt.freqPPM;
175rf.phasePPM = opt.phasePPM;
176rf.deadTime = sys.rfDeadTime;
177rf.ringdownTime = sys.rfRingdownTime;
178rf.delay = opt.delay;
179rf.center = mr.calcRfCenter(rf);
180if ~isempty(opt.use)
181 rf.use=opt.use;
182end
183if rf.deadTime > rf.delay
184 rf.delay = rf.deadTime;
185end
187if nargout > 1
188 assert(opt.sliceThickness > 0,'SliceThickness must be provided');
189 if opt.maxGrad > 0
190 sys.maxGrad = opt.maxGrad;
191 end
192 if opt.maxSlew > 0
193 sys.maxSlew = opt.maxSlew;
194 end
196 amplitude = BW/opt.sliceThickness;
197 area = amplitude*opt.duration;
198 gz = mr.makeTrapezoid('z', sys, 'flatTime', opt.duration, ...
199 'flatArea', area);
200 gzr= mr.makeTrapezoid('z', sys, 'Area', -area*(1-rf.center/rf.shape_dur)-0.5*(gz.area-area));
201 if rf.delay > gz.riseTime
202 gz.delay = ceil((rf.delay - gz.riseTime)/sys.gradRasterTime)*sys.gradRasterTime; % round-up to gradient raster
203 end
204 if rf.delay < (gz.riseTime+gz.delay)
205 rf.delay = gz.riseTime+gz.delay; % these are on the grad raster already which is coarser
206 end
207end
209% v1.4 finally eliminates RF zerofilling
210% if rf.ringdownTime > 0
211% tFill = (1:round(rf.ringdownTime/1e-6))*1e-6; % Round to microsecond
212% rf.t = [rf.t rf.t(end)+tFill];
213% rf.signal = [rf.signal, zeros(size(tFill))];
214% end
215if nargout > 3
216 delay=mr.makeDelay(mr.calcDuration(rf)); % calcDuration already includes the ringdown time
217end
219% RF amplitude check
220rf_amplitude=max(abs(rf.signal));
221if rf_amplitude>sys.maxB1
222 warning('WARNING: system maximum RF amplitude exceeded (%.01f%%)', rf_amplitude/sys.maxB1*100);
223end
225end