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 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: .
module tcai1 use adj_modreal, 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: . Module tcaf1 ( Transient Convolution, Adjoint is the Filter 1-D) does the job.
module tcaf1 use adj_modreal, 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