[Rd] convolve: request for "usual" behaviour + some improvements + some fixes
Herve Pages
hpages at fhcrc.org
Sat Feb 3 02:03:26 CET 2007
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.
More information about the R-devel
mailing list