Transient convolution

When the matrix has a special form, such as  

$$\left[ \begin{array} {c} y_1 \ y_2 \ y_3 \ y_4... ...x_1 \ x_2 \ x_3 \ x_4 \ x_5 \end{array} \right]$$ (1)

then the matrix multiplication and transpose multiplication still fit easily in the same computational framework. The operation $\bold B\bold x$ convolves b_t with x_t, whereas the operation $\bold B'\bold y$crosscorrelates b_t with y_t. I will leave it to you to verify the pseudocode

		do ix = 1, nx {
		do ib = 1, nb {
		 		iy = ib + ix - 1
		 		if operator itself  (convolution)
		 		 		y(iy) = y(iy) + b(ib) $\times$ x(ix)
		 		if adjoint    (correlation)
		 		 		x(ix) = x(ix) + b(ib) $\times$ y(iy)
		      }}

Again, notice that the ``bottom line'' in the program is that x and y are simply interchanged.

Equation (1) could be rewritten as  

$$\left[ \begin{array} {c} y_1 \ y_2 \ y_3 \ y_4... ... \begin{array} {c} b_1 \ b_2 \ b_3 \end{array} \right]$$ (2)

which we abbreviate by $\bold y = \bold X \bold b$.So we can choose between $\bold y = \bold X \bold b$and $\bold y=\bold B\bold x$.In one case the output $\bold y$is dual to the filter $\bold b$,and in the other case the output $\bold y$is dual to the input $\bold x$.In applications, sometimes we will solve for $\bold b$ and sometimes for $\bold x$;so sometimes we will use equation (2) and sometimes (1).

The program contran() can be used with either equation (1) or equation (2), because the calling program can swap the xx and bb variables. The name contran() denotes convolution with ``transpose'' and with ``transient'' end effects.  

#       Convolve and correlate (adjoint convolve).
#
subroutine contran( adj, add, nx, xx, nb, bb, yy)
integer ix, ib, ny, adj, add, nx,     nb
real    xx(nx)                  # input signal
real    bb(nb)                  # filter      (or output crosscorrelation)
real    yy(nx+nb-1)             # filtered signal (or second input signal)
ny = nx + nb - 1                # length of filtered signal
call adjnull(       adj, add,         bb, nb, yy, ny)
do ib= 1, nb {
do ix= 1, nx {
        if( adj == 0 )
                yy( ib+ix-1) = yy( ib+ix-1) + xx( ix) * bb(ib)
        else
                bb( ib)      = bb( ib)      + xx( ix) * yy( ib+ix-1)
        }}
return; end