R devel - Feb 2007 - Inaccuracy in ?convolve

Hi,

Man page for 'convolve' says:

    conj: logical; if 'TRUE', take the complex _conjugate_ before
          back-transforming (default, and used for usual convolution).

The complex conjugate of 'x', of 'y', of both?

In fact it seems that it takes the complex conjugate of 'y' only which
is OK but might be worth mentioning because (1) conj=TRUE is the default
and (2) with this default then convolve(x,y) is not the same as convolve(y,x).
Note that 'convolve' is commutative with conj=FALSE and would be too
if conj=TRUE was taking the complex conjugate of x _and_ y...

Cheers,
H.

Hi again,

There are many problems with current 'convolve' function.
The author of the man page seems to be aware that 'convolve' does _not_
the
usual thing:

  Note that the usual definition of convolution of two sequences 'x'
  and 'y' is given by 'convolve(x, rev(y), type =
"o")'.

and indeed, it does not:

  > x <- 1:3
  > y <- 3:1
  > convolve(x, y, type="o")
  [1]  1  4 10 12  9

The "usual" convolution would rather give:

  > convolve(x, rev(y), type="o")
  [1]  3  8 14  8  3

Also the "usual" convolution is commutative:

  > convolve(y, rev(x), type="o")
  [1]  3  8 14  8  3

but convolve is not:

  > convolve(y, x, type="o")
  [1]  9 12 10  4  1

Of course I could write the following wrapper:

  usual_convolve <- function(x, y, ...) convolve(x, rev(y))

to work around those issues but 'convolve' has other problems:

  (1) The input sequences shouldn't need to have the same length when
      type = "circular" (the shortest can be right-padded with 0s up
      to the length of the longest).
  (2) If the input sequences are both integer vectors, then the result
      should be an integer vector too.
  (3) The "filter" feature seems to be broken (it's not even clear
      what it should do or why we need it though):
        > x <- 1:9
        > y <- 1
        > convolve(x, y, type="f")
        Error in convolve(x, y, type = "f") : subscript out of bounds
        > convolve(y, x, type="f")
        numeric(0)
  (4) If you look at the source code, you'll see that 'x' is first
left-padded
      with '0's. The "usual" convolution doesn't do that:
it always padd
      sequences on the _same_ side (generally on the right).
  (5) It's not clear why we need a 'conj' arg. All what it does is
      take the conjugate of fft(y) before it does the product with fft(x).
      But this has the "non-usual" effect of reverting the expected
result:
        > round(convolve(as.integer(c(0,0,0,1)), 1:7, type="o"))
        [1] 0 0 0 7 6 5 4 3 2 1

Here below is my version of 'convolve' just in case. It does the
"usual"
convolution plus:
  - no need to have 'x' and 'y' of the same length when
'type' is "circular",
  - when 'x' and 'y' are integer vectors, the output is an
integer vector,
  - no more 'conj' arg (not needed, only leads to confusion),
  - when type is "filter", the output sequence is the same as with
    type="open" but is truncated to the length of 'x' (the
original signal)
    It can be seen has the result of 'x' filtered by 'y' (the
filter).

convolve2 <- function(x, y, type = c("circular", "open",
"filter"))
{
    type <- match.arg(type)
    nx <- length(x)
    ny <- length(y)
    if (type == "circular")
        nz <- max(nx, ny)
    else
        nz <- nx + ny - 1
    if (nz > nx)
        x[(nx+1):nz] <- as.integer(0)
    if (nz > ny)
        y[(ny+1):nz] <- as.integer(0)
    fx <- fft(x)
    fy <- fft(y)
    fz <- fx * fy
    z <- fft(fz, inverse=TRUE) / nz
    if (is.numeric(x) && is.numeric(y))
        z <- Re(z)
    if (is.integer(x) && is.integer(y))
        z <- as.integer(round(z))
    if (type == "filter")
        z[1:nx]
    else
        z
}

Cheers,
H.