1function [tc, ic, fi]=calcRfCenter(rf)
2%calcRfCenter Calculate the effective center time of an RF pulse.
3%
4% PURPOSE
5% Returns the time point of the effective rotation of an RF pulse,
6% used to align echo timing, set off-resonance phase compensation, and
7% compute TE/TI delays. For shaped pulses this is the peak of the RF
8% amplitude; for block pulses it is the midpoint of the constant
9% plateau. Zero-padding in rf.signal is treated as part of the shape.
10% The rf.delay field is NOT included in the returned time, so callers
11% that need an absolute time within a block typically compute
12% rf.delay + mr.calcRfCenter(rf).
13%
14% SIGNATURES
15% tc = mr.calcRfCenter(rf) % time of center, seconds
16% [tc, ic] = mr.calcRfCenter(rf) % also integer sample index
17% [tc, ic, fi] = mr.calcRfCenter(rf) % also fractional offset
18%
19% If rf has a .center field (set by the modern RF constructors
20% mr.makeSincPulse, mr.makeBlockPulse, mr.makeGaussPulse,
21% mr.makeArbitraryRf, mr.makeSLRpulse, mr.makeAdiabaticPulse), tc is
22% taken directly from rf.center and ic is the nearest index in rf.t.
23% Otherwise tc is computed from the amplitude peak of rf.signal.
24%
25% INPUTS
26% rf [required] struct, RF event struct. Must have fields .t (seconds, time
27% axis on the RF raster) and .signal (complex Hz, waveform).
28% If field .center (seconds) is present it is used directly;
29% otherwise the function detects the peak of abs(rf.signal).
30% Typically produced by mr.makeSincPulse, mr.makeBlockPulse,
31% mr.makeGaussPulse, mr.makeArbitraryRf, mr.makeSLRpulse,
32% or mr.makeAdiabaticPulse.
33%
34% OUTPUT
35% tc double, seconds, time of the RF center relative to the start
36% of the RF shape (rf.delay not included)
37% ic integer, 1-based index into rf.t / rf.signal of the sample
38% nearest to tc
39% fi double, dimensionless, fractional offset in [-0.5, 0.5] from
40% rf.t(ic) toward the previous (negative) or next (positive)
41% sample, normalized by the local raster step. 0 when tc lies
42% exactly on rf.t(ic) (within 1 ns)
43%
44% NOTES
45% - When rf.center is absent, the peak detector treats samples within
46% 0.001% of max(abs(rf.signal)) as part of the same plateau and
47% returns the midpoint. This is what makes block pulses (constant
48% amplitude) yield a center at the middle of the pulse rather than
49% at the first sample.
50% - rf.delay is intentionally excluded from tc. If a sequence places
51% an RF event after a delay, the absolute time within the block is
52% rf.delay + tc.
53% - Returned tc is calculated as a floating-point value in seconds
54% relative to the beginning of the RF shape (the leading edge of the
55% first RF raster cell of the shape); ic is the index of the raster
56% cell (in Matlab indexing convention), which center is the closest to
57% tc; fi captures the sub-raster offset when rf.center does not
58% coincide with a sample (center of the RF raster cell) and is
59% constrained to the range [-0.5 0.5).
60%
61% EXAMPLE
62% sys = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
63% 'MaxSlew', 170, 'SlewUnit', 'T/m/s', ...
64% 'rfDeadTime', 100e-6); % bumps rf.delay to 100 us
65% rf = mr.makeSincPulse(pi/2, 'Duration', 3e-3, ...
66% 'SliceThickness', 3e-3, 'system', sys);
67% % off-resonance phase compensation for an off-center slice:
68% rf.freqOffset = sys.gamma * 1e-3 * 5e-3; % 5 mm offset
69% rf.phaseOffset = -2*pi*rf.freqOffset*mr.calcRfCenter(rf);
70% % absolute time of the RF center within its block:
71% tCenter = rf.delay + mr.calcRfCenter(rf); % 100 us + 1.5 ms = 1.6 ms
72%
73% SEE ALSO
74% mr.calcRfBandwidth, mr.calcRfPower, mr.makeSincPulse,
75% mr.makeBlockPulse, mr.makeGaussPulse, mr.makeArbitraryRf,
76% mr.makeSLRpulse, mr.makeAdiabaticPulse
77%
79% % detect zero-padding
80% last=length(rf.signal);
81% for first=1:last
82% if abs(rf.signal(first))>eps
83% break;
84% end
85% end
86% for last=last:-1:first
87% if abs(rf.signal(last))>eps
88% break;
89% end
90% end
92% rfmax=max(abs(rf.signal(first:last)));
93% ipeak=find(abs(rf.signal(first:last))>=rfmax-eps);
95 if isfield(rf,'center')
96 tc=rf.center;
97 [~,ic]=min(abs(rf.t-tc));
98 else
100 % we detect the excitation peak and if it is a plato we take its center
101 rfmax=max(abs(rf.signal));
102 ipeak=find(abs(rf.signal)>=rfmax*0.99999);
103 tc=(rf.t(ipeak(1))+rf.t(ipeak(end)))/2;
104 ic=ipeak(round(end/2));
105 end
107 ft=tc-rf.t(ic);
108 if ic<length(rf.t) && ft>1e-9 % 1 ns
109 fi=ft/(rf.t(ic+1)-rf.t(ic));
110 elseif ic>1 && ft<-1e-9 % -1 ns
111 fi=ft/(rf.t(ic)-rf.t(ic-1));
112 else
113 fi=0;
114 end
116% % detect the excitation peak (this code is far from being ideal...)
117% rfmin=min(abs(rf.signal(first:last))); % pure max check fails for the block pulse!!!
118% [rfmax,ic]=max(abs(rf.signal(first:last)));
119% if (rfmax-rfmin)<=eps
120% ic=round((last-first+1)/2); % we take the center of the pulse for block pulses
121% tc=(rf.t(first)+rf.t(last))/2;
122% else
123% tc=rf.t(first-1+ic);
124% end
125end