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

[Herve Pages]

Last but not least: convolve2 can be made 100 times or 1000 times faster
than convolve by choosing a power of 2 for the length of the fft-buffer
(a length of 2^n is the best case for the fft, the worst case being when
the length is a prime number):

> x <- 1:100003
> y <- 1:1
> system.time(cc <- convolve(x, y, type="o")) # uses buffer length of 100003
   user  system elapsed
 76.428   0.016  76.445
> system.time(cc <- convolve2(x, y, type="o")) # uses buffer length of 2^17
   user  system elapsed
  0.164   0.012   0.179

Here is the modified 'convolve2':

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 {
        nz0 <- nx + ny - 1
        nz <- 2^ceiling(log2(nz0))
    }
    if (nz > nx)
        x[(nx+1):nz] <- as.integer(0)
    if (nz > ny)
        y[(ny+1):nz] <- as.integer(0)
    fz <- fft(x) * fft(y)
    z <- fft(fz, inverse=TRUE) / nz
    if (type == "open") {
        z <- z[1:nz0]
    } else {
        if (type == "filter")
            z <- z[1:nx]
    }
    if (is.numeric(x) && is.numeric(y))
        z <- Re(z)
    if (is.integer(x) && is.integer(y))
        z <- as.integer(round(z))
    z
}

In fact, it should try to be smarter than that and not use the fft at all
when one of the 2 input sequences is very short (less than 3 or 4) or
e.g. when one is 10000 times shorter than the other one.

Cheers,
H.


Herve Pages wrote:
> 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.
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>