Finished updates to examples

Finished current round of updates to the example scripts for clarity.
Added dynamic variable flip example as well

Author: agentmess

examples/demo_C13_multiband.m

@@ -25,8 +25,10 @@
 clc
 
 %% multiband pulse
-fprintf(1, '\nHere''s a C13 multiband excitation pulse example, for [1-13C]pyr+13C-urea\n');
-fprintf(1, 'dynamic MR spectroscopic or chemical shift imaging on a 3T clinical system\n\n');
+fprintf(1, '************************************************************\n')
+fprintf(1, 'Here''s a C13 multiband excitation pulse example, for [1-13C]pyr+13C-urea\n');
+fprintf(1, 'dynamic MR spectroscopic or chemical shift imaging on a 3T clinical system\n');
+fprintf(1, '************************************************************\n')
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
@@ -70,10 +72,12 @@
 pause;
 
 %% multiband pulse - minimum phase
-fprintf(1, '\nThis pulse is shortened by using a Minimum-phase spectral filter\n');
+fprintf(1, '************************************************************\n')
+fprintf(1, 'This pulse is shortened by using a Minimum-phase spectral filter\n');
 fprintf(1, 'This may add some phase offset between metabolites and may slightly distort lineshapes\n');
 fprintf(1, 'but under most circumstances won''t reduce SNR and allows for shorter TEs\n\n');
-fprintf(1, 'Here is the resulting pulse with a linear-phase (which is longer and may have higher peak power)\n\n');
+fprintf(1, 'Here is the resulting pulse with a linear-phase (which is longer and may have higher peak power)\n');
+fprintf(1, '************************************************************\n')
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
@@ -91,12 +95,14 @@
 s_ftype = 'min';
 
 %%  minimum phase - spectral correction
-fprintf(1, '\nThis pulse also is corrected for chemical-shift slice misregistration\n');
+fprintf(1, '************************************************************\n')
+fprintf(1, 'This pulse also is corrected for chemical-shift slice misregistration\n');
 fprintf(1, 'by using the ''Spect Correct'' option.\n');
 fprintf(1, 'If the frequency specification bandwidth is too large and/or not sparse, \n');
 fprintf(1, 'the spectral correction can fail and/or result in high RF pulse powers.\n\n');
 fprintf(1, 'Here is the resulting pulse without spectral correction\n');
-fprintf(1, '(the slices for different frequencies are slightly shifted and distorted).\n\n');
+fprintf(1, '(the slices for different frequencies are slightly shifted and distorted).\n');
+fprintf(1, '************************************************************\n')
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
@@ -113,8 +119,10 @@
 
 opt = ss_opt({'Spect Correct', 1}); 
 %% multiband pulse - thicker slice
-fprintf(1, '\nIncreasing slice thickness from 5mm to 1cm will allow for a higher\n');
-fprintf(1, 'spatial time-bandwidth and a sharper slice profile\n\n');
+fprintf(1, '************************************************************\n')
+fprintf(1, 'Increasing slice thickness from 5mm to 1cm will allow for a higher\n');
+fprintf(1, 'spatial time-bandwidth and a sharper slice profile\n');
+fprintf(1, '************************************************************\n')
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
@@ -132,9 +140,11 @@
 
 %% short duration pyruvate-lactate pulse
 
-fprintf(1, '\nHere''s a short C13 multiband excitation pulse, for [1-13C]pyr/lac imaging on a 3T clinical system\n');
+fprintf(1, '************************************************************\n')
+fprintf(1, 'Here''s a short C13 multiband excitation pulse, for [1-13C]pyr/lac imaging on a 3T clinical system\n');
 fprintf(1, 'with unspecified flip angles for other metabolites for cancer imaging applications.\n');
-fprintf(1, 'This type of design is used in UCSF clinical MRSI studies.\n\n');
+fprintf(1, 'This type of design is used in UCSF clinical MRSI studies.\n');
+fprintf(1, '************************************************************\n')
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
@@ -158,7 +168,7 @@
 % create vectors of angles, ripples, and band edges for input to pulse design
 [fspec, a_angs, d] = create_freq_specs(mets);
 fctr = 0;  % force pulse design to optimize for center of frequency specification
-s_ftype = 'lin';  % linear-phase spectral filter
+s_ftype = 'min';  % linear-phase spectral filter
 
 % SPATIAL PULSE PARAMETERS
 z_thk = .5;  % thickness (cm)
@@ -171,3 +181,81 @@
     ss_design(z_thk, z_tb, [z_d1 z_d2], fspec, a_angs, d, ptype, ...
 	      z_ftype, s_ftype, ss_type, fctr);
 set(gcf,'Name', '[1-13C]pyr+lac Simple Multiband');
+
+
+%% short duration pyruvate-lactate pulse with variable flip angles
+
+fprintf(1, '************************************************************\n')
+fprintf(1, 'To improve the SNR for dynamic imaging, the flip angle should be varied over time.\n');
+fprintf(1, 'Here is an example to design a set of RF pulses, with identical gradients.\n\n');
+fprintf(1, 'and different variable flip angles for pyruvate and lactate.\n');
+fprintf(1, '************************************************************\n')
+fprintf(1,'Hit any key to continue: (design is slower since a set of pulses)\n');
+pause;
+
+clear all; ss_opt([]); ss_globals;				% Reset all options
+
+% GENERAL PULSE PARAMETERS
+ss_type = 'EP Whole';
+ptype = 'ex';  % excitation pulse
+opt = ss_opt({'Nucleus', 'Carbon', ...
+	      'Max Duration', 8e-3, ...
+	      'Max B1', 0.5, ...
+          'Spect Correct', 1});
+      
+% SPECTRAL PULSE PARAMETERS 
+B0 = 3e4; % G
+df = 0.5e-6 * B0 * SS_GAMMA; % 0.5 ppm = gamma_C13 * B0 * 0.5e-6
+% metabolite			frequency (Hz)		freq bandwidth (Hz)
+mets(1).name = 'pyr'; 	mets(1).f = -230; 	mets(1).df = 2*df;
+mets(2).name = 'lac'; 	mets(2).f = 165; 	mets(2).df = 2*df; 
+
+% define variable flip angles
+% (see hyperpolarized-mri-toolbox for more info and code on schemes)
+TR = 3; N = 16; 
+kPLest = .05;  R1est = [1/35 1/30]; % [T1pyr, T1lac], 1/s
+
+flips = zeros(2,N);
+% pyruvate flips for constant signal
+% vfa_const_amp(N, pi/2, exp(-TR * (R1(1)+kPL)))
+E1 = exp(-TR * (R1est(1)+kPLest));
+flips(1,N) = pi/2;
+for n = N-1:-1:1
+    flips(1,n) = atan(E1*sin(flips(1,n+1)));
+end
+
+% lactate flips for maximum total SNR
+% vfa_opt_signal(N, exp(-TR * R1(2)));
+E1 = exp(-TR * R1est(2));
+flips(2,:) = acos( sqrt((E1^2-E1.^(2.*(N-[1:N]+1))) ./ (1-E1.^(2*(N-[1:N]+1)))) );
+
+
+dr = .02; % here ripple is a fraction of the flip angle
+
+for n = 1:length(mets)
+    fspec(2*n-1) = mets(n).f - mets(n).df;
+    fspec(2*n) = mets(n).f + mets(n).df;
+end
+fmid = (mets(1).f+mets(end).f)/2;
+fspec = fspec - fmid;
+
+fctr = 0;  % force pulse design to optimize for center of frequency specification
+s_ftype = 'min';  % linear-phase spectral filter
+
+% SPATIAL PULSE PARAMETERS
+z_thk = .5;  % thickness (cm)
+z_tb = 3; % time-bandwidth, proportional to profile sharpness
+z_ftype='ls';  % least-squares filter design
+z_d1 = 0.01;  z_d2 = 0.01;  % slice profile pass and stop-band ripples, respectively
+
+% DESIGN THE PULSE!
+[g,rf,fs,z,f,mxy] = ...
+    ss_design_dyn(z_thk, z_tb, [z_d1 z_d2], fspec, flips, dr, ptype, ...
+	      z_ftype, s_ftype, ss_type, fctr);
+
+% % for saving dynamic pulses:
+% for t = 1:N
+%     ss_save_dyn(g(:,t),rf(:,t),max(flips(:,t)),z_thk, [], 'GE', fspec, flips(:,t), root_fname, t);
+% end
+
+

examples/demo_general.m

@@ -83,7 +83,7 @@
 fprintf(1, '\n************************************************************\n')
 fprintf(1, 'Here''s an example of a water/fat spectral spatial pulse for 1.5T\n\n');
 fprintf(1, 'The pulse is a ''Flyback Half'' type, meaning that it\n');
-fprintf(1, 'uses a flyback trajectory, and has a asymmetric frequency\n');
+fprintf(1, 'uses a flyback trajectory, and has a symmetric frequency\n');
 fprintf(1, 'response.\n');
 fprintf(1, 'The frequency spec includes a passband at [%3f,%3f]\n',fspec(3),fspec(4));
 fprintf(1, 'and stopbands at [%3f,%3f]\n',fspec(1),fspec(2));
@@ -226,7 +226,8 @@
 fprintf(1,'Hit any key to continue:\n');
 pause;
 
-opt = ss_opt({'B1 Verse', 0});
+opt = ss_opt({'B1 Verse', 0, ...
+	      'Max B1', 0.2});
 
 %%
 
@@ -252,17 +253,17 @@
 fprintf(1, 'The package also allows for echo-planar (''EP'') trajectory designs\n');
 fprintf(1, 'which can reduce the pulse duration and allow for thinner slices\n');
 fprintf(1, 'but are more sensitive to eddy currents and timing errors.\n');
-fprintf(1, 'Here is an example of the original design using a EP trajectory\n');
-fprintf(1, 'There is some more ripple in the stopband, but it still meets the spec\n');
+fprintf(1, 'Here is an example of a fat-water EP trajectory design\n');
+fprintf(1, 'The slice is thinner with less stopband ripple, which the EP trajectory can achieve.\n');
+fprintf(1, 'However, with this initial design the ripple specification is not met in the stopband...\n');
 fprintf(1, '************************************************************\n');
 
 ss_opt(default_opt);
-%opt = ss_opt({'Max B1', 0.4});
 ss_type = 'EP Whole';
 f_ctr = [];
-z_thk = .5;
 
-%d = [0.01 0.005];
+z_thk = .5;
+d = [0.005 0.005];
 
 [g_ew,rf_ew,fs,z,f,mxy] = ...
     ss_design(z_thk, z_tb, [z_d1 z_d2], fspec, a*ang, d, ptype, ...
@@ -298,17 +299,19 @@
 fprintf(1, '\n************************************************************\n');
 fprintf(1, 'There is a more sophisticated solution....\n');
 fprintf(1, 'The ripple arises from the nonuniform sampling in time\n');
-fprintf(1, 'that occurs for EP trajectories\n');
+fprintf(1, 'that occurs for EP trajectories.\n');
+fprintf(1, 'This can be corrected for by using the ''Spect Correct'' option.\n');
 fprintf(1, '************************************************************\n');
 
 ss_type = 'EP Whole';
 f_ctr = [];
-ss_opt({'Spect Correct', 1});
+ss_opt({'Spect Correct', 1, ...
+    'Spect Correct Reg', 0.001});  % small amount of regularization added to better condition spectral correction inversion
+
 [g_fw,rf_fw,fs,z,f,mxy] = ...
     ss_design(z_thk, z_tb, [z_d1 z_d2], fspec, a*ang, d, ptype, ...
 	      z_ftype, s_ftype, ss_type, f_ctr, dbg);
 set(gcf,'Name', 'Water/Fat EP Whole - Spect Correct');
 
-fprintf(1,'Hit any key to continue:\n');
-pause;
+
 

ss_design_dyn.m

@@ -543,6 +543,7 @@
     for t = [1, round(Nt/2) Nt]
         [f_plot,z_plot,m_plot] = ss_plot(g(:,t),rf(:,t),SS_TS, ptype,z_thk*2,bw,...
             SS_GAMMA, fmid);
+        set(gcf, 'Name', ['Pulse #' int2str(t)])
     end
 end