%%%% Tutorial about the math of convolution.
%%%%
%%%% Written by Rajeev Raizada, Aug.19, 2002.
%%%% 
%%%% This program goes together with hrf_tutorial.m and
%%%% design_matrix_tutorial.m

% The math of convolution:
% So far we have just talked about a stimulus "kicking off" an HRF,
% and the HRFs just adding up on top of each other.
% The actual math of this is just multiplying and adding.
% The multiplying part is: we go along the stimulus-time-series vector,
% element by element, and we multiply each element by the HRF.
% e.g. if the vector starts off: [ 0 0 1 0 ... ]
% then we take multiply: 0*HRF, then 0*HRF, then 1*HRF, then 0*HRF etc.
%
% We then take those multiplied HRFs and add them all together,
% to make an output vector that is the result of the convolution.
%
% The only extra consideration is that each HRF gets "kicked off",
% i.e. multiplied and added, at a different position along the
% stimulus-time-series vector. 
% In our [ 0 0 1 0 ] example above, the 1 is in the third position
% of the stimulus-time-series vector, so we want our 1*HRF to get 
% added onto the same position of our output vector,
% i.e. starting at the third position.
%
% To make things easier, let's try this with a small HRF,
% that we talk more about in design_matrix_tutorial.m

hrf = [ 0  4  2  -1  0 ];

% And let's start with a simple stimulus-time-series vector, corresponding
% to a single stimulus coming on at the third time-point 
% (which will be the third TR, in the case of fMRI).

stim_time_series = [ 0 0 1 0 0 ];

% Now, we want to make an output vector that is the result of
% convolving the stimulus-time-series vector with the HRF.
% So, we go along the stimulus-time-series vector, element by element,
% multiply each element by the HRF, and add the result to
% the corresponding position of our output vector.
%
%        0*[ 0  4  2 -1  0 ]           +
%        
%           0*[ 0  4  2 -1  0 ]        +
%         
%              1*[ 0  4  2 -1  0 ]     +
%         
%                 0*[ 0  4  2 -1  0 ]  +
%
%                    0*[ 0  4  2 -1  0 ] 
%        
% Let's see how this would all add up for making 
% the fifth element in our output:
%
%        0*[ 0  4  2 -1  0 ]           +
%                        .        
%           0*[ 0  4  2 -1  0 ]        +
%                        .
%              1*[ 0  4  2 -1  0 ]     +
%                        .
%                 0*[ 0  4  2 -1  0 ]  +
%                        .
%                    0*[ 0  4  2 -1  0 ] 
%                        .  
% makes the result: 0*0 + 0*-1 + 1*2 +0*4 + 0*0 = 2
%
% In the example above, the 5th element of our output vector
% is equal to:
%
% 1st-element of the stimulus-time-series  *  5th-element of HRF  +
% 2nd-element of the stimulus-time-series  *  4th-element of HRF  +
% 3rd-element of the stimulus-time-series  *  3rd-element of HRF  +
% 4th-element of the stimulus-time-series  *  2nd-element of HRF  +
% 5th-element of the stimulus-time-series  *  1st-element of HRF 
%
% As an equation, this is:
% Output(pos_in_output) =  
%   Sum[ with pos_in_hrf going from 1 to length(hrf) ] of
%     stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
%
% Textbooks will often use a notation like this:
% Output(n) means the n-th element of the output of the result
% of our convolution.
% Typically people write the vector to be convolved as x
% In our case, x is stim_time_series.
% So, x(n) means the n-th element of the onsets-vector, x.
% People often call the vector that this gets convolved with the "filter".
% Another name for it is the "impulse response function",
% which is where the phrase "hemodynamic response function" comes from.
%
% In our case, the filter is the HRF.
% People often call our position in the HRF "k", and the length
% of the HRF "M".
% In that notation:
% Output(n) = Sum[ with k going from 0 to M ]  x(n-k) * hrf(k)
%
% In Matlab, the indices can't be quite as neat as this, because
% a vector doesn't have a zero-th entry. 
% The first entry is element number 1 in the vector,
% even though it corresponds to time=0 in the HRF.
% That's why in the equation
%
%   Output(pos_in_output) =  
%     Sum[ with pos_in_hrf going from 1 to length(hrf) ] of
%       stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
%
% we see that pos_in_hrf goes from 1 to length(hrf), rather than starting
% at zero, and also why there's a "+1" in 
% stim_time_series(pos_in_output - pos_in_hrf + 1)
%
% By the way, note that this means that the result of 
% convolving stim_time_series  with hrf is longer than 
% either stim_time_series or hrf.
% length(output) = length(stim_time_series) + length(hrf) - 1.
% This makes sense when you think of it in terms of a stimulus
% kicking off an HRF that runs to completion.
% Even if the stimulus is right at the very end of a scan and
% if we stop the scan immediately after showing the stimulus,
% it will still have kicked-off an HRF inside the subject's head.
% The result of the convolution is the overall HRF-response that
% the stimuli will evoke, and so if the last stimulus kicks off 
% a fresh HRF, that will continue by one HRF-length past the
% end of the stimulus-time-series.

%%%%%%%%% Using a for-loop to implement the convolution %%%%%%%%%%%%%%
%
% So, now we have the basic equation for the convolution:
%
%   Output(pos_in_output) =  
%     Sum[ with pos_in_hrf going from 1 to length(hrf) ] of
%       stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
%
% and we have to write a computer program that will plug the actual numbers
% into that equation and see what they add up to.
%
% In order to do that, we want to go along the each position in the output,
% one position at a time, and then go along the positions in the HRF,
% carrying out the multiplication:
%    stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
% for each position in the HRF, and then adding up the results of
% all of these multiplications.
%
% In Matlab, and in almost any programming language, the easiest way
% to go along a series of positions in a list, one step at a time,
% is using a for-loop.
% (A for-loop isn't always the most elegant or most efficient method, but
%  it's often the simplest way, so that's how we'll do it here).
%
% Here's how we'd write a for-loop to do the convolution.
%
% Start off with an output vector all zeros, 
% with one row and (length(stim_time_series) + length(hrf) - 1) columns
output = zeros(1,length(stim_time_series) + length(hrf) - 1);           

%%%% Now build two loops.
%%%% The main loop goes along positions in the output vector.
%%%% Nested inside that, we have a smaller loop that goes along
%%%% positions in the HRF.
for pos_in_output = 1:length(output),   
              %%% Go around length(output) times, with the value of
              %%% pos_in_output starting at 1, and increasing by 1 
              %%% each time we go around the loop. 
    
   for pos_in_hrf = 1:length(hrf), 
              %%% This loop is completely inside the pos_in_output loop.
              %%% So, for any given pos_in_output value, we loop
              %%% through all the positions in the HRF, going from 1
              %%% to length(hrf), which is the last element in the HRF.
       
      % Now we want to do summing as we move along the positions in the HRF
      % This involves looking back to the earler positions of 
      % stim_time_series, back to position (pos_in_output - pos_in_hrf + 1).
      % However, if we're at the beginning of the stim_time_series,
      % then these earlier positions don't exist, so they don't
      % play any role in the convolution. 
      % In Matlab, trying to access non-existent parts of the stim_time_series
      % vector would cause an error, so first we check that we're far enough
      % along. Similarly, we can't try to access elements in stim_time_series
      % that happen after the end of that vector.
      
      % Output(pos_in_output) =  
      %   Sum[ with pos_in_hrf going from 1 to length(hrf) ] of
      %     stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
      pos_in_stim_time_series = pos_in_output - pos_in_hrf + 1;
      
      if (pos_in_stim_time_series > 0) & ...  % If it's far enough along and...
         (pos_in_stim_time_series <= length(stim_time_series)),
                                              % it's not too far along.
                                              % <= means less than or equal to
         output(pos_in_output) = output(pos_in_output) + ...
            stim_time_series(pos_in_stim_time_series) * hrf(pos_in_hrf);

      end;  % End of the if-statement
        
   end;     % End of loop along the HRF

end;        % End of loop along the positions in the output vector


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Let's display the output in the Matlab command window, by
% typing it without a semi-colon at the end:

output

% We can check that this is correct by using Matlab's built-in
% convolution function. This should be the same as the output above!

output_from_matlab_conv_function = conv(stim_time_series,hrf)

% You can also express convolution as a matrix multiplication,
% but I'm not going to get into that here.
% If your Matlab has the Signal Processing Toolbox, type
% help convmtx
% for more info.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Ok, that seemed to work.
% Now let's try the same code with a slightly less trivial
% convolution, and let's plot the results.
%
% Instead of just having one stimulus, at time t=3, 
% we'll have several stimuli, and randomly jittered times

stim_time_series = [ 0 0 1 0 1 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 ];

output = zeros(1,length(stim_time_series) + length(hrf) - 1);           

%%%%% Below is exactly the same code as above, but with
%%%%% some of the lengthier comments taken out, to save repetition.
for pos_in_output = 1:length(output),   
    
   for pos_in_hrf = 1:length(hrf), 
 
      % Output(pos_in_output) =  
      %   Sum[ with pos_in_hrf going from 1 to length(hrf) ] of
      %     stim_time_series(pos_in_output - pos_in_hrf + 1) * hrf(pos_in_hrf)
      pos_in_stim_time_series = pos_in_output - pos_in_hrf + 1;
      
      if (pos_in_stim_time_series > 0) & ...  % If it's far enough along and...
         (pos_in_stim_time_series <= length(stim_time_series)),
                                              % it's not too far along.
                                              % <= means less than or equal to
         output(pos_in_output) = output(pos_in_output) + ...
            stim_time_series(pos_in_stim_time_series) * hrf(pos_in_hrf);

      end;  % End of the if-statement        
   end;     % End of loop along the HRF
end;        % End of loop along the positions in the output vector

% We can check that this is correct by using Matlab's built-in
% convolution function. This should be the same as the output above!
output_from_matlab_conv_function = conv(stim_time_series,hrf);

%%%%%%%%%%%%% Let's plot the results 

figure(1);
clf;
subplot(3,1,1);    % This lines up three subplots on this figure window
stem(stim_time_series,'bo-'); 
grid on;
legend('Stimulus time-series to be convolved');
axis([0 23 0 1.8]);
ylabel('Stimulus present / absent');

subplot(3,1,2);
plot(output,'rx-'); % Red crosses, solid line 
grid on;
legend('Convolved with HRF using for-loop');
axis([0 23 -2 12]);
ylabel('Predicted fMRI signal');

subplot(3,1,3);
plot(output_from_matlab_conv_function,'k*-'); % Black stars, solid line 
grid on;
legend('Convolved with HRF using built-in Matlab function');
axis([0 23 -2 12]);
ylabel('Predicted fMRI signal');