1function [rf, gz, gzr, delay] = makeSincPulse(flip,varargin)
2%makeSincPulse Create a slice selective since pulse.
3% rf=makeSincPulse(flip, 'Duration', dur) Create sinc pulse
4% with given flip angle (rad) and duration (s).
5%
6% rf=makeSincPulse(..., 'freqOffset', f,'phaseOffset',p)
7% Create sinc pulse with frequency offset (Hz) and phase offset (rad).
8%
9% rf=makeSincPulse(..., 'freqPPM',-3.3)
10% Create arbitrary RF pulse with frequency offset specified in PPM (e.g.
11% actual frequency offset proportional to the true Larmor frequency), in
12% this example -3.3 ppm as often used for fat saturation; can be combined
13% with the 'freqOffset' specified in Hz.
14%
15% [rf, gz]=makeSincPulse(...,'SliceThickness',st) Return the
16% slice select gradient corresponding to given slice thickness (m).
17%
18% [rf, gz]=makeSincPulse(flip,lims,...) Create slice selection gradient
19% with the specificed gradient limits (e.g. amplitude, slew).
20%
21% [rf, gz, gzr]=makeSincPulse(flip,lims,...) Create slice selection and
22% slice refocusing gradients with the specificed gradient limits
23% (e.g. amplitude, slew) and taking into account 'centerpos' parameter
24%
25% See also Sequence.addBlock
27validPulseUses = mr.getSupportedRfUse();
29persistent parser
30if isempty(parser)
31 parser = mr.aux.InputParserCompat;
32 parser.FunctionName = 'makeSincPulse';
34 % RF params
35 addRequired(parser, 'flipAngle', @isnumeric);
36 addOptional(parser, 'system', [], @isstruct);
37 %addParamValue(parser, 'system', [], @isstruct);
38 addParamValue(parser, 'duration', 0, @isnumeric);
39 addParamValue(parser, 'freqOffset', 0, @isnumeric);
40 addParamValue(parser, 'phaseOffset', 0, @isnumeric);
41 addParamValue(parser, 'freqPPM', 0, @isnumeric);
42 addParamValue(parser, 'phasePPM', 0, @isnumeric);
43 addParamValue(parser, 'timeBwProduct', 4, @isnumeric);
44 addParamValue(parser, 'apodization', 0, @isnumeric);
45 addParamValue(parser, 'centerpos', 0.5, @isnumeric);
46 % Slice params
47 addParamValue(parser, 'maxGrad', 0, @isnumeric);
48 addParamValue(parser, 'maxSlew', 0, @isnumeric);
49 addParamValue(parser, 'sliceThickness', 0, @isnumeric);
50 addParamValue(parser, 'delay', 0, @isnumeric);
51 addParamValue(parser, 'dwell', 0, @isnumeric); % dummy default value
52 % whether it is a refocusing pulse (for k-space calculation)
53 addParamValue(parser, 'use', 'u', @(x) any(validatestring(x,validPulseUses)));
54end
55parse(parser, flip, varargin{:});
56opt = parser.Results;
58if isempty(opt.system)
59 system=mr.opts();
60else
61 system=opt.system;
62end
64if opt.dwell==0
65 opt.dwell=system.rfRasterTime;
66end
68if opt.duration<=0
69 error('rf pulse duration must be positive');
70end
72BW = opt.timeBwProduct/opt.duration;
73alpha = opt.apodization;
74N = round(opt.duration/opt.dwell);
75t = ((1:N)-0.5)'*opt.dwell;
76tt = t - opt.duration*opt.centerpos;
77window = (1.0-alpha+alpha*cos(2*pi*tt/opt.duration));
78signal = window.*sinc(BW*tt);
79flip = sum(signal)*opt.dwell*2*pi;
80signal = signal*opt.flipAngle/flip;
82rf.type = 'rf';
83rf.signal = signal;
84rf.t = t;
85rf.shape_dur=N*opt.dwell;
86rf.freqOffset = opt.freqOffset;
87rf.phaseOffset = opt.phaseOffset;
88rf.freqPPM = opt.freqPPM;
89rf.phasePPM = opt.phasePPM;
90rf.deadTime = system.rfDeadTime;
91rf.ringdownTime = system.rfRingdownTime;
92rf.delay = opt.delay;
93if ~isempty(opt.use)
94 rf.use=opt.use;
95end
96if rf.deadTime > rf.delay
97 rf.delay = rf.deadTime;
98end
99rf.center=opt.duration*opt.centerpos;
101if nargout > 1
102 assert(opt.sliceThickness > 0,'SliceThickness must be provided');
103 if opt.maxGrad > 0
104 system.maxGrad = opt.maxGrad;
105 end
106 if opt.maxSlew > 0
107 system.maxSlew = opt.maxSlew;
108 end
110 amplitude = BW/opt.sliceThickness;
111 area = amplitude*opt.duration;
112 gz = mr.makeTrapezoid('z', system, 'flatTime', opt.duration, ...
113 'flatArea', area);
114 gzr= mr.makeTrapezoid('z', system, 'Area', -area*(1-opt.centerpos)-0.5*(gz.area-area));
115 if rf.delay > gz.riseTime
116 gz.delay = ceil((rf.delay - gz.riseTime)/system.gradRasterTime)*system.gradRasterTime; % round-up to gradient raster
117 end
118 if rf.delay < (gz.riseTime+gz.delay)
119 rf.delay = gz.riseTime+gz.delay; % these are on the grad raster already which is coarser
120 end
121end
123% v1.4 finally eliminates RF zerofilling
124% if rf.ringdownTime > 0
125% tFill = (1:round(rf.ringdownTime/1e-6))*1e-6; % Round to microsecond
126% rf.t = [rf.t rf.t(end)+tFill];
127% rf.signal = [rf.signal, zeros(size(tFill))];
128% end
129if nargout > 3
130 delay=mr.makeDelay(mr.calcDuration(rf)); % calcDuration already includes the ringdown time
131end
133% RF amplitude check
134rf_amplitude=max(abs(rf.signal));
135if rf_amplitude>system.maxB1
136 warning('WARNING: system maximum RF amplitude exceeded (%.01f%%)', rf_amplitude/system.maxB1*100);
137end
139function y = sinc(x)
140 % sinc Calculate the sinc function:
141 % sinc(x) = sin(pi*x)/(pi*x)
142 %
143 % This is a useful helper function for those without the signal
144 % processing toolbox
146 i = find(x == 0);
147 x(i) = 1;
148 y = sin(pi*x)./(pi*x);
149 y(i) = 1;
150end
152end