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