78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1function grad = makeExtendedTrapezoid(channel, varargin)
2%makeExtendedTrapezoid Create an extended trapezoid gradient event.
3%
4% g = makeExtendedTrapezoid(channel, lims, ...
5% 'times', times, ...
6% 'amplitudes', amplitudes)
7% Create a gradient by specifying a set of points (amplitudes) at
8% specified time points(times) at a given channel with given system
9% limits. This function returns an extended or arbitrary gradient object,
10% the latter if 'convert2arbitrary' is set to true.
11%
12% See also Sequence.addBlock mr.opts makeTrapezoid
13%
14% Stefan Kroboth <stefan.kroboth@uniklinik-freiburg.de>
16persistent parser
18if isempty(parser)
19 validChannels = {'x', 'y', 'z'};
20 parser = mr.aux.InputParserCompat;
21 parser.FunctionName = 'makeExtendedTrapezoid';
22 parser.addRequired('channel', ...
23 @(x) any(validatestring(x, validChannels)));
24 parser.addOptional('system',[],@isstruct);
25 parser.addParamValue('times', 0, @isnumeric);
26 parser.addParamValue('amplitudes', 0, @isnumeric);
27 parser.addParamValue('maxGrad', 0, @isnumeric);
28 parser.addParamValue('maxSlew', 0, @isnumeric);
29 parser.addParamValue('skip_check', false);
30 parser.addParamValue('convert2arbitrary', false);
32end
33parse(parser,channel,varargin{:});
34opt = parser.Results;
36if isempty(opt.system)
37 system=mr.opts();
38else
39 system=opt.system;
40end
42if any(size(opt.times) ~= size(opt.amplitudes))
43 error('Times and amplitudes must have the same length.');
44end
46if all(opt.times == 0)
47 error('At least one of the given times must be non-zero.');
48end
50if any(diff(opt.times)<=0)
51 error('Times must be in ascending order and all times must be distinct.');
52end
54if abs(round(opt.times(end)/system.gradRasterTime)*system.gradRasterTime-opt.times(end))>1e-8 % 10ns is an acceptable rounding error
55 error('The last time point must be on a gradient raster.');
56end
58%if all(opt.amplitudes == 0)
59% error('At least one of the given amplitudes must be non-zero.');
60%end
62if opt.skip_check == false && opt.times(1) > 0 && opt.amplitudes(1) ~= 0
63 error('If first amplitude of a gradient is nonzero, it must connect to previous block!');
64end
66maxSlew = system.maxSlew;
67maxGrad = system.maxGrad;
68if opt.maxGrad > 0
69 maxGrad = opt.maxGrad;
70end
71if opt.maxSlew > 0
72 maxSlew = opt.maxSlew;
73end
75if (opt.convert2arbitrary)
76 % represent the extended trapezoid on the regularly sampled time grid
77 waveform = mr.pts2waveform(opt.times, opt.amplitudes, system.gradRasterTime);
78 grad = mr.makeArbitraryGrad(channel, waveform, system, ...
79 'maxSlew', maxSlew,...
80 'maxGrad', maxGrad,...
81 'delay', opt.times(1));
82else
83 % keep the original possibly irregular sampling
84 if any(abs(round(opt.times/system.gradRasterTime)*system.gradRasterTime-opt.times)>1e-8) % 10ns is an acceptable rounding error
85 error('All time points must be on a gradient raster or "convert2arbitrary" option must be used.');
86 end
87 % check slew rate and gradient amplitude against the active limits
88 % the convert2arbitrary branch above gets these checks via mr.makeArbitraryGrad
89 slew = (opt.amplitudes(2:end)-opt.amplitudes(1:end-1)) ./ (opt.times(2:end)-opt.times(1:end-1));
90 if ~isempty(slew) && max(abs(slew))>maxSlew
91 error('Slew rate violation (%.0f%%)',max(abs(slew))/maxSlew*100);
92 end
93 if max(abs(opt.amplitudes))>maxGrad
94 error('Gradient amplitude violation (%.0f%%)',max(abs(opt.amplitudes))/maxGrad*100);
95 end
96 %
97 grad.type = 'grad';
98 grad.channel = opt.channel;
99 grad.waveform = opt.amplitudes(:);
100 grad.delay = round(opt.times(1)/system.gradRasterTime)*system.gradRasterTime;
101 grad.tt = opt.times(:) - grad.delay;
102 grad.shape_dur = round(grad.tt(end)/system.gradRasterTime)*system.gradRasterTime;
103 grad.area=0.5*sum((grad.tt(2:end)-grad.tt(1:end-1)).*(grad.waveform(2:end)+grad.waveform(1:end-1)));
104end
106% MZ: although makeArbitraryGrad sets the .first and .last for extended
107% trapezoids we can do it better
108grad.first=opt.amplitudes(1);
109grad.last=opt.amplitudes(end);
111end
115% figure; plot(waveform)
116% waveform = waveform(1:end-2);
117% opt.times = round(opt.times/system.gradRasterTime)*system.gradRasterTime; % round onto grid
118% times_diff = diff(opt.times);
119% amplitudes_diff = diff(opt.amplitudes);
120% waveform = [];
121% for ii = 1:length(opt.times)-1
122% % SK: there are no new points after the end, therefore we dont need to
123% % handle the overlap situation.
124% if ii == length(opt.times)-1
125% crop = 0;
126% else
127% crop = system.gradRasterTime;
128% end
129% y = amplitudes_diff(ii)/times_diff(ii)*...
130% (0:system.gradRasterTime:(opt.times(ii+1)-opt.times(ii)-crop))...
131% + opt.amplitudes(ii);
132% waveform = [waveform y(1:end)];
133% end