1function [rf, gz, gzr, delay] = makeArbitraryRf(signal,flip,varargin)
2%makeArbitraryRf Create an RF pulse with the given pulse shape.
3% rf=makeArbitraryRf(singal, flip) Create RF pulse with complex signal
4% and given flip angle (in radians)
5%
6% rf=makeArbitraryRf(..., 'freqOffset', f,'phaseOffset',p)
7% Create arbitrary RF pulse with frequency offset and phase offset.
8%
9% rf=makeArbitraryRf(..., 'ppmOffset')
10% Create arbitrary RF pulse with frequency offset specified in PPM (e.g.
11% actual frequency offset proportional to the true Larmor frequency); can
12% be combined with the 'freqOffset' specified in Hz.
13%
14% [rf, gz]=makeArbitraryRf(..., 'Bandwidth', bw, 'SliceThickness', st)
15% Create RF pulse and corresponding slice select gradient. The bandwidth
16% of the pulse must be given for the specified shape. You can also check
17% yourself after creating the pulse object by calling mr.calcRfBandwidth()
18%
19% See also mr.calcRfBandwidth mr.makeSincPulse, Sequence.addBlock
21validPulseUses = mr.getSupportedRfUse();
23persistent parser
24if isempty(parser)
25 parser = mr.aux.InputParserCompat;
26 parser.FunctionName = 'makeArbitraryRf';
28 % RF params
29 addRequired(parser, 'signal', @isnumeric);
30 addRequired(parser, 'flipAngle', @isnumeric);
31 addOptional(parser, 'system', [], @isstruct);
32 addParamValue(parser, 'freqOffset', 0, @isnumeric);
33 addParamValue(parser, 'phaseOffset', 0, @isnumeric);
34 addParamValue(parser, 'freqPPM', 0, @isnumeric);
35 addParamValue(parser, 'phasePPM', 0, @isnumeric);
36 addParamValue(parser, 'timeBwProduct', 0, @isnumeric);
37 addParamValue(parser, 'bandwidth', 0, @isnumeric);
38 addParamValue(parser, 'center', NaN, @isnumeric);
39 % Slice params
40 addParamValue(parser, 'maxGrad', 0, @isnumeric);
41 addParamValue(parser, 'maxSlew', 0, @isnumeric);
42 addParamValue(parser, 'sliceThickness', 0, @isnumeric);
43 % Delay
44 addParamValue(parser, 'delay', 0, @isnumeric);
45 addParamValue(parser, 'dwell', 0, @isnumeric); % dummy default value
46 % whether it is a refocusing pulse (for k-space calculation)
47 addParamValue(parser, 'use', 'u', @(x) any(validatestring(x,validPulseUses)));
48end
49parse(parser, signal, flip,varargin{:});
50opt = parser.Results;
52if isempty(opt.system)
53 system=mr.opts();
54else
55 system=opt.system;
56end
58if opt.dwell==0
59 opt.dwell=system.rfRasterTime;
60end
62signal = signal./abs(sum(signal.*opt.dwell))*flip/(2*pi);
64if size(signal,1)>size(signal,2)
65 signal=signal.';
66end
68N= length(signal);
69duration = N*opt.dwell;
70t = ((1:N)'-0.5)*opt.dwell;
72rf.type = 'rf';
73rf.signal = signal(:);
74rf.t = t;
75rf.shape_dur=duration;
76rf.freqOffset = opt.freqOffset;
77rf.phaseOffset = opt.phaseOffset;
78rf.freqPPM = opt.freqPPM;
79rf.phasePPM = opt.phasePPM;
80rf.deadTime = system.rfDeadTime;
81rf.ringdownTime = system.rfRingdownTime;
82rf.delay = opt.delay;
83if ~isempty(opt.use)
84 rf.use=opt.use;
85end
86if rf.deadTime > rf.delay
87 rf.delay = rf.deadTime;
88end
90if isfinite(opt.center)
91 rf.center=opt.center;
92 if rf.center < 0, rf.center = 0; end
93 if rf.center > rf.shape_dur, rf.center = rf.shape_dur; end
94else
95 rf.center = mr.calcRfCenter(rf);
96end
98if opt.timeBwProduct>0
99 if opt.bandwidth > 0
100 error('Both ''bandwidth'' and ''timeBwProduct'' cannot be specified at the same time');
101 else
102 opt.bandwidth=opt.timeBwProduct/duration; % QL
103 end
104end
106if nargout>1
107 assert(opt.sliceThickness > 0, 'SliceThickness must be provided');
108 assert(opt.bandwidth > 0 || opt.timeBwProduct > 0, 'Bandwidth or BW-time-product of the pulse must be provided');
109 if opt.maxGrad > 0
110 system.maxGrad = opt.maxGrad;
111 end
112 if opt.maxSlew > 0
113 system.maxSlew = opt.maxSlew;
114 end
116 BW = opt.bandwidth;
117 if opt.timeBwProduct > 0
118 BW = opt.timeBwProduct/duration;
119 end
121 amplitude = BW/opt.sliceThickness;
122 area = amplitude*duration;
123 gz = mr.makeTrapezoid('z', system, 'flatTime', duration, ...
124 'flatArea', area);
126 if rf.delay > gz.riseTime
127 gz.delay = ceil((rf.delay - gz.riseTime)/system.gradRasterTime)*system.gradRasterTime; % round-up to gradient raster
128 end
129 if rf.delay < (gz.riseTime+gz.delay)
130 rf.delay = gz.riseTime+gz.delay; % these are on the grad raster already which is coarser
131 end
133 if nargout > 2
134 gzr= mr.makeTrapezoid('z', system, 'Area', -area*(1-rf.center)/rf.shape_dur-0.5*(gz.area-area));
135 end
136end
138% v1.4 finally eliminates RF zerofilling
139% if rf.ringdownTime > 0
140% tFill = (1:round(rf.ringdownTime/1e-6))*1e-6; % Round to microsecond
141% rf.t = [rf.t rf.t(end)+tFill];
142% rf.signal = [rf.signal, zeros(size(tFill))];
143% end
144if nargout > 3
145 delay=mr.makeDelay(mr.calcDuration(rf)); % calcDuration already includes the ringdown time
146end
148% RF amplitude check
149rf_amplitude=max(abs(rf.signal));
150if rf_amplitude>system.maxB1
151 warning('WARNING: system maximum RF amplitude exceeded (%.01f%%)', rf_amplitude/system.maxB1*100);
152end
154end