Up: Annotated Program Listings
findintconvo (``find integrals of convolutions'') computes the
at many different values of k2, for up to three weighting functions
wi. Here, T1 and T2 represent the start and end time any
particular frame; in concrete terms, the inner integrand is just the
across each frame, and is
computed once for each frame. The outer integral then integrates across
all frames to arrive at a single value for each weighting function
wi and each value of k2.
- The function declaration line.
function [int1, int2, int3] = ...
findintconvo (Ca_even, ts_even, k2_lookup,...
midftimes, flengths, w1, w2, w3, progress)
The input arguments are:
- the blood data resampled at an evenly-spaced time
domain. Usually, this is as returned by resampleblood, and
possibly delay- and dispersion-corrected by correctblood.
- the time domain at which Ca_even was resampled.
- the list of k2 values at which to evaluate equation
. The output will consist of one value per
weighting function per value of k2.
- the mid-frame times (used to perform the
- the length of each frame (also used in the
- the three weighting functions wi. w1 must be
supplied, but if it is an empty matrix, then a weighting function of
1 is assumed. w2 and w3 are not required, and if they
are not given, then int2 and int3 (respectively) will
not be defined in the output. Generally, the weighting funtions are
simple functions of t such as 1, t, or t2.
- an optional argument indicating whether progress should
be reported (in the form of printing a dot when the weighted
integrals for each value of k2 have been computed). This is
quite useful, as findintconvo is one of the most
time-consuming steps of RCBF analysis.
The output arguments int1, int2, and int3 are
simply vectors containing the value of equation 17 for
the three different weighting functions; each element of these vectors
corresponds to one value of k2 as supplied in k2_lookup.
- Do some initial setup. We need to get the sizes of various
vectors, and initialize vectors that we will fill element by element
later. Initializing vectors to all zero before filling them allows
better memory management by MATLAB.
NumEvenTimes = length(ts_even);
NumFrames = length(midftimes);
fstart = midftimes - (flengths / 2);
TableSize = length (k2_lookup);
integrand = zeros (NumFrames, 1);
if (nargin >= 6); int1 = zeros (1, TableSize); end;
if (nargin >= 7); int2 = zeros (1, TableSize); end;
if (nargin == 8); int3 = zeros (1, TableSize); end;
% if w1 is empty, assume that it should be all ones
if isempty (w1)
w1 = ones (size(NumFrames));
- Calculate each element of the integrals, one at a time.
Unfortunately, there does not seem to be any way to vectorize this
operation, and it must therefore be performed within an inefficient
for i = 1:TableSize
- First calculate the innermost terms of equation
17. exp_fun is
e-k2 t, and convo
We then perform the inner
integration, i.e. average across all frames. This is performed by
nframeint, a special CMEX routine that has been carefully
optimised for this very common operation in RCBF analysis. Type
help nframeints in MATLAB for more information.
exp_fun = exp(-k2_lookup(i) * ts_even);
convo = nconv(Ca_even, exp_fun, ts_even(2) - ts_even(1));
integrand = nframeint (ts_even, convo(1:length(ts_even)),...
- Find any frames that are not completely covered by the
ts_even time domain. nframeint actually does this
for us, returning NaN (not-a-number) in any such element of its
output vector. Thus, we only use data from frames that are not equal to NaN, selected using the built-in isnan
select = ~isnan(integrand);
- Now we calculate the outer integral (the weighted integral
across all frames) for each weighting function with the current
value of k2. Again, we use a CMEX routine to enhance
performance; ntrapz is a replacement for MATLAB's
trapz function that is much faster, as well as having an
optional argument to allow weighting of the integrand. Note the use
of select to make sure that we only use frames actually
spanned by the blood data.
int1 (i) = ntrapz(midftimes(select), integrand(select), w1(select));
int2 (i) = ntrapz(midftimes(select), integrand(select), w2(select));
int3 (i) = ntrapz(midftimes(select), integrand(select), w3(select));
(We have omitted the code that makes w2 and w3
optional in favour of concentrating on the numerical analysis.)
- End the for loop.
Up: Annotated Program Listings