% compute_convolution.m

%

% function [conv,iflag] = ...

%    compute_convolution(dt,t0,f,r,i_r_zero,iplot);

%

% This MATLAB function computes (and optionally plots)

% the convolution between a discretely-sampled input

% function and a response function.

%

% Input

% -----

% dt = the sampling interval

% f = the input signal

% r = the response function

% i_r_zero = the index of r corresponding to t = 0

% iplot = if non-zero, plot signal, response, and convolution

%

% K. Beers. MIT ChE. 8/14/2002

function [fCr,iflag] = ...

    compute_convolution(dt,t0,f,r,i_r_zero,iplot);

iflag = 0;

% Make sure both input signal and response are

% row vectors and not column vectors

if(size(f,1) ~= 1)

    f = f';

end

if(size(r,1) ~= 1)

    r = r';

end

% Next, determine the length of zero padding that

% must be added to the input signal.  To do so,

% we need to determine the number of non-zero values

% of the response function at postive and negative times.

% the number of non-zero r(t) values at t >= 0

num_r = length(r);

num_r_pos = num_r - i_r_zero + 1;

num_r_neg = i_r_zero;

num_pad = max(num_r_pos,num_r_neg);

% We now pad the input signal with this number of zeros.

% Also, at this time, we pad with zeros to bring the

% total number of points up to a power of 2.

N = length(f);  % original signal length

N_pad = N + num_pad;  % length with padding for r(t) length

two_power = log2(N_pad);

if(rem(two_power,1))  % two_power is not an integer

    N_pad = pow2(floor(two_power)+1);  % pad to get power of 2

end

num_pad = N_pad - length(f);  % total number of zeros to pad

f = [f zeros(1,num_pad)];  % padded input signal

% we now generate the response signal in wrap-around form

r_old = r;

r = zeros(1,N_pad);

r(1:num_r_pos) = r_old(1:num_r_pos);

for k=0:(num_r_neg-1)

    r(N_pad-k) = r_old(num_r-k-1);

end

% we now compute the fft of each signal

f_FT = fft(f);

r_FT = fft(r);

% the product of the Fourier transforms is calculated

fCr_FT = f_FT.*r_FT;

% the convolution is then obtained by an inverse

% Fourier transform

fCr = real(ifft(fCr_FT));

% if requested, plots of the signal, response, and

% convolution are calculated

if(exist('iplot'))

    if(iplot)

        figure;

        t_val = linspace(0,(N_pad-1)*dt,N_pad);

        t_plot_min = 0;

        t_plot_max = t_val(N_pad);

        % input signal

        subplot(2,2,1);

        plot(t_val,f);

        y_buffer = 0.05*(max(f)-min(f));

        axis([t_plot_min t_plot_max ...

                (min(f)-y_buffer) (max(f)+y_buffer)]);

        xlabel('t');

        ylabel('f(r)');

        title('Computing convolution');

        % response function

        subplot(2,2,2);

        plot(t_val,r);

        xlabel('t');

        ylabel('r(t)');

        y_buffer = 0.05*(max(r)-min(r));

        axis([t_plot_min t_plot_max ...

                (min(r)-y_buffer) (max(r)+y_buffer)]);

        % convolution

        subplot(2,2,3);

        plot(t_val,fCr);

        xlabel('t');

        ylabel('convolution f*r');

        y_buffer = 0.05*(max(fCr)-min(fCr));

        axis([t_plot_min t_plot_max ...

                (min(fCr)-y_buffer) (max(fCr)+y_buffer)]);

    end

end

iflag = 1;

return;