1function [total_energy, peak_pwr, rf_rms]=calcRfPower(rf, dt)
2%calcRfPower Compute relative energy, peak power, and RMS amplitude of an RF pulse.
3%
4% PURPOSE
5% Computes the relative energy, peak power, and RMS B1 amplitude of an
6% RF pulse by resampling its waveform onto a uniform time grid and
7% integrating |rf|^2. Returns 'relative' quantities in units of Hz^2
8% and Hz, not absolute SI power, because Pulseq RF amplitudes are
9% stored in Hz (gamma-scaled). Used for pulse-level sanity checks and
10% as a building block for the sequence-level mr.Sequence/calcRfPower.
11%
12% SIGNATURES
13% total_energy = mr.calcRfPower(rf) % default dt=1 us
14% [total_energy, peak_pwr] = mr.calcRfPower(rf)
15% [total_energy, peak_pwr, rf_rms] = mr.calcRfPower(rf)
16% [...] = mr.calcRfPower(rf, dt) % override sampling step
17%
18% The two arguments must be passed positionally and in order; they
19% cannot be given as name/value pairs. Outputs beyond nargout are not
20% computed.
21%
22% INPUTS
23% rf [required] struct, RF event struct from mr.makeSincPulse,
24% mr.makeBlockPulse, mr.makeGaussPulse, mr.makeArbitraryRf,
25% mr.makeSLRpulse, or mr.makeAdiabaticPulse. Must have
26% fields .t (seconds), .signal (complex Hz), and .shape_dur
27% (seconds).
28% dt [optional] double, resampling step in seconds. Default: 1e-6.
29%
30% OUTPUT
31% total_energy double, Hz (= Hz^2 * s), integral of |rf|^2 over the pulse duration
32% peak_pwr double, Hz^2, maximum of |rf|^2 over the resampled waveform
33% rf_rms double, Hz, RMS B1 amplitude, sqrt(total_energy/rf.shape_dur)
34%
35% NOTES
36% - Outputs are 'relative': amplitude is in Hz, so total_energy is in
37% Hz^2 * s and peak_pwr in Hz^2. To convert to SI: divide rf_rms by
38% sys.gamma to get Tesla; divide total_energy by sys.gamma^2 to get
39% Tesla^2 * s. Absolute SAR requires further scaling by tx-coil and
40% subject-specific factors (reference voltage, coil design) and is
41% not computed here.
42% - The pulse is resampled onto a uniform grid of step dt with bin
43% midpoints at ((0:nn-1)+0.5)*dt where nn = round(rf.shape_dur/dt).
44% Samples outside rf.t are filled with 0 by linear extrapolation.
45% For pulses already sampled on a uniform raster (sinc, Gauss, SLR,
46% arbitrary, adiabatic) the default dt=1e-6 oversamples mildly; for
47% block pulses (rf.t has only the two endpoints) resampling is
48% required.
49% - peak_pwr depends on dt because the resampled waveform may not hit
50% the original peak exactly; for narrow pulses, reduce dt to tighten
51% the estimate.
52% - rf_rms uses rf.shape_dur, which excludes any rf.delay or
53% post-pulse ringdown. For SAR-relevant duty cycle over a sequence
54% use mr.Sequence/calcRfPower instead.
55% - rf.freqOffset, rf.phaseOffset, and rf.freqPPM are ignored: the
56% calculation works on |rf.signal|^2, which is invariant under
57% frequency or phase modulation.
58%
59% EXAMPLE
60% sys = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
61% 'MaxSlew', 170, 'SlewUnit', 'T/m/s');
62% rf = mr.makeSincPulse(pi/2, 'Duration', 3e-3, ...
63% 'SliceThickness', 3e-3, 'apodization', 0.5, ...
64% 'timeBwProduct', 4, 'system', sys);
65% [total_energy, peak_pwr, rf_rms] = mr.calcRfPower(rf);
66% fprintf('energy = %.3g Hz, peak = %.3g Hz^2, rms = %.3g Hz\n', ...
67% total_energy, peak_pwr, rf_rms);
68%
69% SEE ALSO
70% mr.calcRfBandwidth, mr.calcRfCenter, mr.simRf, mr.makeSincPulse,
71% mr.makeBlockPulse, mr.makeGaussPulse, mr.makeArbitraryRf,
72% mr.makeSLRpulse, mr.makeAdiabaticPulse
74if nargin<2
75 dt=1e-6; % for now default sampling rate is 1Mhz
76end
78% resample the pulse to a resonable time array
79nn=round(rf.shape_dur/dt);
80t=((0:(nn-1))+0.5)*dt;
81rfs=interp1(rf.t,rf.signal,t,'linear',0);
82% TODO: avoid resampling of already sampled pulses (above)
84rfs_sq=rfs.*conj(rfs);
85total_energy=sum(rfs_sq)*dt;
86if nargout>1
87 peak_pwr=max(rfs_sq);
88 if nargout>2
89 rf_rms=sqrt(total_energy/rf.shape_dur);
90 end
91end