Convolution equation and program

What do we actually do in a computer when we multiply two Z-transforms together? The filter 2 + Z would be represented in a computer by the storage in memory of the coefficients (2, 1). Likewise, for 1 - Z, the numbers (1, -1) would be stored. The polynomial multiplication program should take these inputs and produce the sequence (2, -1, -1). Let us see how the computation proceeds in a general case, say

$$\begin{eqnarray} X(Z)\, B(Z) &\eq & Y(Z) \ (x_0 + x_1 Z + x_2 Z^2 + \cdots )\, (b_0 + b_1 Z + b_2 Z^2) &\eq & y_0 + y_1 Z + y_2 Z^2 + \cdots\end{eqnarray}$$ (5) (6)

Identifying coefficients of successive powers of Z, we get

$$\begin{eqnarray} y_0 &\eq & x_0 b_0 \nonumber \ y_1 &\eq & x_1 b_0 + x_0 b_1 \n... ...nonumber \ &\eq & \cdots\cdots\cdots\cdots\cdots\cdots \nonumber\end{eqnarray}$$ (7)

In matrix form this looks like  

$$\left[ \begin{array} {c} y_0 \ y_1 \ y_2 \ y_3... ... \begin{array} {c} b_0 \ b_1 \ b_2 \end{array} \right]$$ (8)

The following equation, called the ``convolution equation,'' carries the spirit of the group shown in (7):  

$$y_k \eq \sum^{N_b}_{i = 0} x_{k - i} b_i$$ (9)

To be correct in detail when we associate equation (9) with the group (7), we should also assert that either the input x_k vanishes before k=0 or N_b must be adjusted so that the sum does not extend before x₀. These end conditions are expressed more conveniently by defining j=k-i in equation (9) and eliminating k getting  

$$y_{j+i} \eq \sum^{N_b}_{i = 0} x_j b_i$$ (10)

A convolution program based on equation (10) including end effects on both ends, is convolve().  

#       convolution:    Y(Z) = X(Z) * B(Z)
#
subroutine convolve( nb, bb, nx, xx, yy )
integer nb      # number of coefficients in filter
integer nx      # number of coefficients in input
                # number of coefficients in output will be nx+nb-1
real    bb(nb)  # filter coefficients
real    xx(nx)  # input trace
real    yy(1)   # output trace
integer ib, ix, iy, ny
ny = nx + nb -1
call null( yy, ny)
do ib= 1, nb
        do ix= 1, nx
                yy( ix+ib-1) = yy( ix+ib-1) + xx(ix) * bb(ib)
return; end

Some details of the Ratfor programming language are given in an appendix, along with the subroutine zero() , which erases the space for the output.