[R] convolve bug?

[Bill Simpson]

Hi Martin,

>     Yudi> I agree with Bill that convolve() should mean convolution
>     Yudi> operation. So the default 'conj' should be set to F (Fourier
>     Yudi> transform of a convolution is a product of the transfroms), and
>     Yudi> the zero padding is on the right.
> 
> I tend to agree....  see below
> 
>     Yudi> E.g., target:  (2,2,3,3,4)*(1,1,2)  =  (2  4  9 10 13 10  8 )
> 
>     Yudi> convolve(c(2,2,3,3,4),c(2,1,1),type='o')
>     Yudi> #[1]  2  4  9 10 13 10  8  
>     Yudi> convolve(c(2,2,3,3,4),c(1,1,2),conj=F,type='o')
>     Yudi> #[1] 10  8  2  4  9 10 13      # almost 
>     Yudi> x _ c(2,2,3,3,4,0,0)   # zero padding
>     Yudi> y _ c(1,1,2,0,0,0,0)
>     Yudi> aa_ fft(fft(x) * fft(y), inv = TRUE)
>     Yudi> Re(aa)/length(aa)
>     Yudi> #[1]  2  4  9 10 13 10  8   # wanted
> (there are indeed several ways to get this result,  using
>  combinations of  Conj(), rev() and/or  0- padding...      
I would strongly suggest that R uses the "standard" way of doing it--as
described in the linear systems and digital signal processing books.
The standard way of doing linear (aperiodic) convolution using fft:
Bracewell(1978, p.32) example:
x={2 2 3 3 4}, h={1 1 2}, x*h= {2  4  9 10 13 10 8} [* is convolution]
length:  m         n                m+n-1
zero padded versions for linear convolution:
x={ 2 2 3 3 4 0 0},    h={1 1 2 0 0 0 0}
add n-1 zeros,         add m-1 zeros
Then linear convolution is: ifft(fft(x)*fft(h))

>     Yudi> I would suggest that the 'conj' option not be provided; when
>     Yudi> conj=F, the type='f' option gives meaningless answer:
> ?? Contradiction:
> Above you say that  conj=F should be the *default* ...
Yudi is right the second time: the conj option should NOT be there.
Instead, there should be TWO functions:
convolve does: fft(fft(a)*fft(b), inv=TRUE)
crosscorr does: fft(fft(a)*fft(Conj(b)), inv=TRUE)
This behaviour is fixed, not alterable.

> I agree that at least    convolve(x,y) 
> should really give what everyone expects..

Here is what I expect linear convolution to do:
myconvolve<-function (x,h) 
{
    nx <- length(x)
    nh <- length(h)
    #zero pad
    if(nx>nh)
       {
       x <- c(x,rep(0, nh-1))
       h<-c(h,rep(0,nx-1))
       }
    else
       {
       h <- c(h,rep(0, nx-1))
       x<-c(x,rep(0,nh-1))
       }
    x <- fft(fft(x) * fft(h), inv = TRUE)
    Re(x)/length(x)
#I am not sure about this, the R IFFT is weird
}

I am sure you R gurus can make this better.
I don't know how this works if x and h are 2D, say.

> One way to achieve this would be a new type, 
> e.g. "regular" , "default" or "polynomial".
> 
> And we would change the *default* 
>    -   from type = "circular" to  type = "regular" (or whatever we call it).
I think a nice pairing of terms is "circular" and "linear". I don't know
what the present "filter" is at all--I thought it was "linear", since I
think of linear convolution as filtering.


>    -   from conj = TRUE  to   conj = type != "circular"
> 	    (conj defaulting to FALSE when type="open", ..)
> 
> or do you think conj should default to FALSE even for circular?
NO, the structure of the code should be changed as I suggested above, so
that conj is gone as an argument. What "circular" should do is just
eliminate the zero-padding:

myconvolve2<-function (x,h)
{
#no padding, circular convolution
    nx <- length(x)
    nh <- length(h)
    x <- fft(fft(x) * fft(h), inv = TRUE)
    Re(x)/(nx)
}

> (I must admit not having seen a definition of a circular convolution in a
>  text book..)
See Oppenheim &Schafer Discrete-time signal processing, 2nd ed p. 548 sec
8.2.5

> Yes, I'm waiting for further propositions, and I agree the defaults should
> be changed.
Great, I have presented my arguments above.

Bill

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._