78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1function [bw,fc,spectrum,f,rfs,t]=calcRfBandwidth(rf, cutoff, df, dt)
2%calcRfBandwidth Compute the bandwidth and spectrum of an RF pulse.
3%
4% PURPOSE
5% Estimates the bandwidth of an RF pulse using a simple FFT under the
6% low-angle (small-tip) approximation. Also returns the pulse's center
7% frequency, complex spectrum, and the corresponding frequency and
8% resampled time axes. Typically used for plotting excitation profile
9% magnitude vs. frequency and for sanity-checking RF pulses produced by
10% mr.makeSincPulse, mr.makeGaussPulse, mr.makeBlockPulse,
11% mr.makeArbitraryRf, mr.makeSLRpulse, and mr.makeAdiabaticPulse.
12%
13% SIGNATURES
14% bw = mr.calcRfBandwidth(rf) % default cutoff=0.5 (FWHM), df=10 Hz, dt=1 us
15% [bw,fc] = mr.calcRfBandwidth(rf) % also return center frequency
16% [bw,fc,spectrum,f] = mr.calcRfBandwidth(rf) % also return spectrum and frequency axis
17% [bw,fc,spectrum,f,rfs,t] = mr.calcRfBandwidth(rf) % also return resampled RF waveform and time axis
18% [...] = mr.calcRfBandwidth(rf, cutoff) % override fractional cutoff
19% [...] = mr.calcRfBandwidth(rf, cutoff, df) % override spectral resolution
20% [...] = mr.calcRfBandwidth(rf, cutoff, df, dt) % override time sampling step
21%
22% Bandwidth is the frequency width at which |spectrum| drops to
23% cutoff*max(|spectrum|), with linear interpolation at the first and
24% last crossings. All four positional arguments must be given in order;
25% they cannot be passed as name/value pairs.
26%
27% INPUTS
28% rf [required] struct, RF event struct from mr.makeSincPulse,
29% mr.makeBlockPulse, mr.makeGaussPulse, mr.makeArbitraryRf,
30% mr.makeSLRpulse, or mr.makeAdiabaticPulse. Must have
31% fields .t (seconds), .signal (complex Hz),
32% .freqOffset (Hz), .freqPPM, .phaseOffset (radians),
33% and .center (seconds).
34% cutoff [optional] double, fractional threshold for bandwidth measurement,
35% dimensionless in (0,1]. 0.5 gives FWHM. Default: 0.5.
36% df [optional] double, spectral resolution in Hz. Default: 10.
37% dt [optional] double, time sampling step in seconds. Default: 1e-6.
38%
39% OUTPUT
40% bw double, Hz, bandwidth at the given cutoff threshold
41% fc double, Hz, center frequency (midpoint of the two threshold crossings)
42% spectrum complex vector, small-tip excitation response, scaled by
43% sin(2*pi*dt*s_ref)/s_ref with s_ref = |spectrum(fc)|
44% f double vector, Hz, frequency axis, step df, length round(1/df/dt)
45% rfs complex vector, RF waveform resampled onto the uniform
46% time grid t, with rf.freqOffset, rf.freqPPM, and
47% rf.phaseOffset folded into the phase modulation
48% t double vector, seconds, uniform time axis centered on
49% rf.center, step dt, span 1/df
50%
51% NOTES
52% - Computed from the RF envelope's Fourier transform (small-tip
53% equivalent). The Bloch-simulated bandwidth grows modestly at
54% large flip angles. For exact slice-profile analysis use mr.simRf.
55% - If rf.freqPPM is non-zero, the function calls mr.opts() to fetch
56% sys.gamma and sys.B0 for the PPM-to-Hz conversion, and emits a
57% warning reminding the caller to set system defaults via
58% mr.opts('setAsDefault', true). The system defaults are used
59% silently otherwise.
60% - The FFT length is round(1/df/dt) (default 1e5 points). Reducing df
61% or dt below their defaults increases compute time quadratically in
62% the product, not in either one alone.
63% - The spectrum is normalized so that |spectrum(fc)| approximates the
64% on-resonance small-tip excitation amplitude rather than an
65% unnormalized FFT magnitude.
66%
67% EXAMPLE
68% sys = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
69% 'MaxSlew', 170, 'SlewUnit', 'T/m/s');
70% rf = mr.makeSincPulse(pi/2, 'Duration', 3e-3, ...
71% 'SliceThickness', 3e-3, 'apodization', 0.5, ...
72% 'timeBwProduct', 4, 'system', sys);
73% [bw, fc, spectrum, f] = mr.calcRfBandwidth(rf);
74% fprintf('BW = %.1f Hz, center = %.2f Hz\n', bw, fc);
75% figure; plot(f, abs(spectrum)); xlim(3*[-bw bw]);
76% xlabel('Frequency, Hz'); title('Excitation pulse profile');
77%
78% SEE ALSO
79% mr.calcRfCenter, mr.calcRfPower, mr.simRf, mr.makeSincPulse,
80% mr.makeGaussPulse, mr.makeBlockPulse, mr.makeArbitraryRf,
81% mr.makeSLRpulse, mr.makeAdiabaticPulse
82%
84if nargin<2
85 cutoff=0.5;
86end
88if nargin<3
89 df=10; % spectral resolution in Hz
90end
92if nargin<4
93 dt=1e-6; % for now default sampling rate is 1Mhz, it's probably too high
94end
96if abs(rf.freqPPM)>eps
97 warning('mr.calcRfBandwidth((): relying on the system properties, like B0 and gamma, stored in the global environment by callimg mr.lims(''setAsDefault'',true)');
98 sys=mr.opts();
99 full_freqOffset=rf.freqOffset+rf.freqPPM*1e-6*sys.gamma*sys.B0;
100else
101 full_freqOffset=rf.freqOffset;
102end
104tc=rf.center;
106% resample the pulse to a resonable time array
107nn=round(1/df/dt);
108t=(-floor(nn/2):ceil(nn/2)-1)*dt;
110rfs=interp1(rf.t-tc,rf.signal.*exp(1i*(rf.phaseOffset+2*pi*full_freqOffset*rf.t)),t,'linear',0);
111spectrum=fftshift(fft(fftshift(rfs)));
112f=(-floor(nn/2):ceil(nn/2)-1)*df;
114w1=mr.aux.findFlank(f,spectrum,cutoff);
115w2=mr.aux.findFlank(f(end:-1:1),spectrum(end:-1:1),cutoff);
117bw=w2-w1;
118fc=(w2+w1)/2;
120% coarse STE scaling -- we normalize to the max of the spectrum, this works
121% better with frequency-shifted pulses than the abs(sum(shape)) -- the
122% 0-frequency response; yes, we could take the spectrum at fc but this
123% would have a problem for non-symmetric pulses...
124%s_ref=max(abs(spectrum));
125s_ref=interp1(f,abs(spectrum),fc);
126spectrum=sin(2*pi*dt*s_ref)*spectrum/s_ref;
128end