Internal convolution

The output of the transient convolution program has the length larger than the length of the input. In many applications, it is convenient to set the two lengths equal to each other. Reformulating the convolution operator allows us to relate the output space to the input space and deal appropriately with the convolution end effects. The Fortran 90 implementation is in module icai1 ( Internal Convolution, Adjoint is the Input 1-D).

module icai1	
  use adj_mod

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

contains

subroutine icai1_init (filter, nd, lag)
    real, dimension (:), pointer   :: filter
    integer,           intent (in) :: nd
    integer, optional, intent (in) :: lag

integer                        :: na

b => filter ; nb = size (b) ; nx = nd

if (present (lag)) then
       lg = lag
    else
       lg = 1
    end if

y1 = nb + 1 - lg
    y2 = nx + 1 - lg
  end subroutine icai1_init

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

integer              :: i, x1, x2

call adjnull (adj, add, x, y)

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

stat = 0
  end function icai1_op

end module icai1

The lag parameter lag, saved at the initialization step to the private variable lg, defines how the output of the operator aligns with the input. The end effects of transient and internal convolutions are compared in Figure 1.

 

conv90 Figure 1 Example of convolution end effects. From top to bottom: (1) input; (2) filter; (3) output of tcai1(), (4) output of icai1() with ( lag=1).