1function [Mz_z,Mz_xy,F,ref_eff,Mx_xy,My_xy]=simRf(rf,rephase_factor,prephase_factor)
2%simRf Simulate an RF pulse with the given pulse shape.
3% [Mz_z,Mz_xy,F,ref_eff,Mx_xy,My_xy]=simRf(pulse,prephase_factor,rephase_factor)
4% Performs a rapid RF pulse simulation based on the rotation formalism.
5% The algorithm is optimized by using quaternions to represent rotations.
6% The compulsory parameter 'rf' is the Pulseq RF pulse. Optional
7% parameter 'rephase_factor' is needed in several cases e.g. to correclty
8% visualize the phase of the magnetization for slice-selective
9% excitation. Another optional parameter 'prephase_factor' is an
10% experimental parameter useful for simulating refocusing pulses or
11% spoiling needed.
12% Return values:
13% Mz_z,Mz_xy: z and xy comnponents of the magnetisation after the pulse
14% assuming the unit magnetization was aligned with z before
15% the pulse. Useful for assessing excitation RF pulses.
16% F: frequency axis in Hz
17% ref_eff: Refocusing efficiency of the pulse as a complex value.
18% Magnitude of ref_eff seems to closely follow Mz_z. Phase
19% of ref_eff is related to the effective phase of the RF
20% pulse, e.g. the axis of the planar flip.
21% Mx_xy,My_xy: xy magnetizations after the RF pulse assuming the unit
22% magnetization was aligned with x or y axis prior to the
23% pulse, respectively. Useful for detailed analyses of
24% refocusing pulses.
25%
26% The implementation was inspired by the example by Dr. Tony Stoecker
27% (https://github.com/stoeckert/mr-simu-example-ismrm19)
28% The algorithm was rewritten to quaternions and vectorized for
29% performance by MZ
30%
32bw_mul=4; % simulation bandwidth (multiplier of the pulse bandwidth)
33df=1; % spectral resolution [Hz]
34dt=10e-6; % (re-)sampling interval
36if nargin < 2
37 if isfield(rf,'use') && strcmp(rf.use,'refocusing')
38 rephase_factor = 0;
39 else
40 rephase_factor = -(rf.shape_dur-rf.center)/rf.shape_dur;
41 end
42end
44if nargin < 3
45 prephase_factor = 0;
46end
48[bw,f0]=mr.calcRfBandwidth(rf,0.5,df*10,dt);
50% our bandwidth here is the total bandwidth relative to 0, so we have to add the center frequency
51bw = abs(bw) + abs(f0);
53% adapt time stepping -- just some compromizes -- we stick to dt~1/bw/50
54if bw>4e3
55 dt=5e-6;
56 if bw>1e4
57 dt=2e-6;
58 if bw>20000
59 dt=1e-6;
60 end
61 end
62end
64T = (1:round(rf.shape_dur/dt))*dt-0.5*dt; % timesteps axis [s]
65F = 2*pi*linspace(f0-bw_mul*bw/2,f0+bw_mul*bw/2,bw/df)'; % offset frequencies [rad/s]
67if abs(rf.freqPPM)>eps || abs(rf.phasePPM)>eps
68 warning('relying on the system properties, like B0 and gamma, stored in the global environment by callimg mr.lims(''setAsDefault'',true)');
69 sys=mr.opts();
70 full_freqOffset=rf.freqOffset+rf.freqPPM*1e-6*sys.gamma*sys.B0;
71 full_phaseOffset=rf.phaseOffset+rf.phasePPM*1e-6*sys.gamma*sys.B0;
72else
73 full_freqOffset=rf.freqOffset;
74 full_phaseOffset=rf.phaseOffset;
75end
77shapea = interp1(rf.t, 2*pi*rf.signal.*exp(1i*(full_phaseOffset+2*pi*full_freqOffset*rf.t)),T,'linear',0);
79% intialize result vectors
80M_ROT=zeros(size(F));
81Z_ROT=zeros(size(F));
82sf=size(F);
83q=zeros(sf(1),4);
84q(:,1)=1; % init rotation quaternions
86% prephaser / left spoiler
87W = -F*dt*length(T)*prephase_factor; % effective field rotation angle
88Q = [cos(W/2) zeros(sf) zeros(sf) sin(W/2)];
89q=quat_multiply(q,Q);
91% RF pulse simulation
92for j=1:length(T)
93 W = -dt*sqrt(abs(shapea(j))^2+F.^2); % effective field rotation angles
94 n = dt * [real(shapea(j))*ones(sf) imag(shapea(j))*ones(sf) F]./abs(W); % effective field rotation axes
95 Q = [cos(W/2) sin(W/2).*n];
96 q=quat_multiply(q,Q);
97end
99% rephaser / right spoiler / refocusing pulse
100W = -F*dt*length(T)*rephase_factor; % effective field rotation angle
101Q = [cos(W/2) zeros(sf) zeros(sf) sin(W/2)];
102q=quat_multiply(q,Q);
104% export results
105F=F/(2*pi);
106m=zeros(sf(1),4);
108% excitation: start with M0=M_z
109m(:,4)=1;
110m0rf=quat_multiply(quat_conj(q),quat_multiply(m,q));
111Mz_z=m0rf(:,4);
112Mz_xy=m0rf(:,2)+1i*m0rf(:,3);
114% refocusing: start both with M0=M_x and them M0=M_y
115m=zeros(sf(1),4);
116m(:,2)=1;
117Mx_xy=quat_multiply(quat_conj(q),quat_multiply(m,q));
118Mx_xy=Mx_xy(:,2)+1i*Mx_xy(:,3);
119m=zeros(sf(1),4);
120m(:,3)=1;
121My_xy=quat_multiply(quat_conj(q),quat_multiply(m,q));
122My_xy=My_xy(:,2)+1i*My_xy(:,3);
123ref_eff=(Mx_xy+My_xy*1i)/2;
124end
126function qout = quat_multiply( q, r )
127% quat_multiply: Calculate the product of two quaternions.
129% Calculate vector portion of quaternion product
130% vec = s1*v2 + s2*v1 + cross(v1,v2)
131vec = [q(:,1).*r(:,2) q(:,1).*r(:,3) q(:,1).*r(:,4)] + ...
132 [r(:,1).*q(:,2) r(:,1).*q(:,3) r(:,1).*q(:,4)]+...
133 [ q(:,3).*r(:,4)-q(:,4).*r(:,3) ...
134 q(:,4).*r(:,2)-q(:,2).*r(:,4) ...
135 q(:,2).*r(:,3)-q(:,3).*r(:,2)];
137% Calculate scalar portion of quaternion product
138% scalar = s1*s2 - dot(v1,v2)
139scalar = q(:,1).*r(:,1) - q(:,2).*r(:,2) - ...
140 q(:,3).*r(:,3) - q(:,4).*r(:,4);
142qout = [scalar vec];
143end
145function q = quat_conj( q )
146% quat_conj Calculate the conjugate of a quaternion.
147q(:,2:4) = -q(:,2:4);
148end