Transient convolution

As with any linear operator, the convolution operator has a representation in the form of matrix multiplication. However, it is more convenient to implement it in a different form. Compare the mathematical definition of convolution $y_{k} = \sum_{i} b_{i} x_{k-i}$ with its implementation in module tcai1. The name of this module abbreviates Transient Convolution, Adjoint is the Input 1-D, which refers to the fact that the adjoint operator cross-correlates the output with the filter to return to the input space: $x_j = \sum_{i} b_{i} y_{i+j}$.

module tcai1	
  use adj_mod

real, dimension (:), pointer, private :: b
  integer, private                      :: nb, nx

contains

subroutine tcai1_init (filter, nd)
    real, dimension (:), pointer :: filter
    integer, intent (in)         :: nd

b => filter
    nb = size (b)
    nx = nd
  end subroutine tcai1_init

function tcai1_op (adj, add, x, y) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: x, y

integer              :: i

call adjnull (adj, add, x, y)

do i = 1, nb
       if (adj) then
          x (1 :     nx    ) = x (1 :     nx    ) + b (i) * y (i : i + nx - 1)
       else
          y (i : i + nx - 1) = y (i : i + nx - 1) + b (i) * x (1 :     nx    )
       end if
    end do

stat = 0

end function tcai1_op

end module tcai1

Analogously to the case of matrix multiplication, we initialize the operator with a pointer to the filter. The pointer is stored in a private variable b to be referenced by function tcai1_op.

In some applications, we need to consider the filter b as the output of the adjoint process. In this case, the adjoint operation corresponds to cross-correlation of a fixed portion of the filter input x with a variable portion the output y: $b_i = \sum_{j} x_j y_{i+j}$. Module tcaf1 ( Transient Convolution, Adjoint is the Filter 1-D) does the job.

module tcaf1	
  use adj_mod

real,    dimension (:), pointer, private :: x
  logical, dimension (:), pointer, private :: m
  integer, private                         :: nx

contains

subroutine tcaf1_init (input)
    real,    dimension (:), pointer :: input

x => input ; nx = size (x)
  end subroutine tcaf1_init

function tcaf1_op (adj, add, b, y) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: b, y

integer              :: i

call adjnull (adj, add, b, y)

do i = 1, size (b)
       if (adj) then
          b (i) = b (i) + sum (y (i : i + nx - 1) * x (1 : nx))
       else
          y (i : i + nx - 1) = y (i : i + nx - 1) + x (1 : nx) * b (i)
       end if
    end do

stat = 0
  end function tcaf1_op

end module tcaf1