CAUSAL INTEGRATION RECURSION

Causal integration is defined as

$$y(t) \eq \int_{-\infty}^t \ x(t)\ dt$$ (23)

Sampling the time axis gives a matrix equation which we should call causal summation, but we often call it causal integration.  

$$\left[ \begin{array} {c} y_0 \ y_1 \ y_2 \ y_3 \ y... ...\ x_5 \ x_6 \ x_7 \ x_8 \ x_9 \ \end{array} \right]$$ (24)

(In some applications the 1 on the diagonal is replaced by 1/2.) Causal integration is the simplest prototype of a recursive operator. The coding is trickier than operators we considered earlier. Notice when you compute y₅ that it is the sum of 6 terms, but that this sum is more quickly computed as y₅=y₄ + x₅. Thus equation (24) is more efficiently thought of as the recursion  

$$y_t \quad = \quad y_{t-1} + x_t \quad \quad \quad {\rm for\ increasing\ } t$$ (25)

(which may also be regarded as a numerical representation of the differential equation dy/dt=x.)

When it comes time to think about the adjoint, however, it is easier to think of equation (24) than of (25). Let the matrix of equation (24) be called $\bold C$.Transposing to get $\bold C '$ and applying it to $\bold y$gives us something back in the space of $\bold x$,namely $\tilde{\bold x} = \bold C' \bold y$.From it we see that the adjoint calculation, if done recursively, needs to be done backwards like  

$$\tilde x_{t-1} \quad = \quad \tilde x_{t} + y_{t-1} \quad \quad \quad {\rm for\ decreasing\ } t$$ (26)

We can sum up by saying that the adjoint of causal integration is anticausal integration.

A subroutine to do these jobs is causint() . The code for anticausal integration is not obvious from the code for integration and the adjoint coding tricks we learned earlier. To understand the adjoint, you need to inspect the detailed form of the expression $\tilde{\bold x} = \bold C' \bold y$and take care to get the ends correct.  

# causal integration  (1's on diagonal)
#
subroutine causint( adj, add,        n,xx,  yy   )
integer i,       n, adj, add;   real xx(n), yy(n )
temporary real tt( n)
call adjnull(       adj, add,        xx,n,  yy,n )
if( adj == 0){                                  tt(1) = xx(1) 
                                        do i= 2, n
                                                tt(i) = tt(i-1) + xx(i)
        do i= 1, n
                yy(i) = yy(i) + tt(i)
        }
else {                                          tt(n) = yy(n)
                                        do i= n, 2, -1 
                                                tt(i-1) = tt(i) + yy(i-1)
        do i= 1, n
                xx(i) = xx(i) + tt(i)
        }
return; end

 

causint Figure 6 in1 is an input pulse. C in1 is its causal integral. C' in1 is the anticausal integral of the pulse. in2 is a separated doublet. Its causal integration is a box and its anticausal integration is the negative. CC in2 is the double causal integral of in2. How can a triangle be built?

Later we will consider equations to march wavefields up towards the earth's surface, a layer at a time, an operator for each layer. Then the adjoint will start from the earth's surface and march down, a layer at a time, into the earth.

EXERCISES:

1. Modify the calculation in Figure 7 to make a triangle waveform on the bottom row.