[Rd] convolve: request for "usual" behaviour + some improvements + some fixes

[Herve Pages]

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.