[R] convolve bug?

[Ben Bolker]

On Thu, 18 Nov 1999, Bill Simpson wrote:


> 1. First example, from the classic Bracewell The Fourier transform and its
> applications, chap 3 (p.32 in 2nd edition):
> {2 2 3 3 4} * {1 1 2} = {2 4 9 10 13 10 8}
> where * denotes convolution (discrete in this case). See bottom of this 
> email for his definition of convolution.
> In R I get:
> x<-c(2,2,3,3,4)
> h<-c(1,1,2)
> convolve(x,h,type="o")
> [1]  4  6 10 11 14  7  4
> I think "open" is the right type. I want to get length(x)+length(h)-1
> terms in the convolution output, which is what "open" gives. I definitely
> don't want "circular"

  The difference between the definitions is a sign change in the index of
the filter: f[i-j] vs. f[j-i] (if the filter is symmetric then it wouldn't
matter).  If you do

convolve(x,rev(h),type="o")

you'll get the answer from Bracewell.  The definition in R makes more
sense to me, but that's personal I guess.

> 
> h<-dnorm(x, mean=50, sd=10)
> # the filter is a Gaussian
> 
> plot(convolve(y,h, type="open"))
> #weird
> points(pnorm(x,mean=50,sd=10))
> #what I expected, since convolution with a step is integration

  Here the problem is that the filter is "running off the end".
If you extend the filter with ones

y <- c(y,rep(1,100))

and redo the analysis you'll see that it levels off at 1 for a while and
then drops back down.

  Hope that helps ...  (as usual, someone else has probably answered this
one already).

  Ben

-- 
--------------------------------------------
Ben Bolker                                  bolker at zoo.ufl.edu
Zoology Department, University of Florida   http://www.zoo.ufl.edu/bolker
318 Carr Hall/Box 118525                    tel: (352) 392-5697
Gainesville, FL 32611-8525                  fax: (352) 392-3704

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._