convolution in R

I tried to do convolution in R directly and using FFTs then taking inverse. But it seems from simple observation it is not correct. Look at this example:

# DIRECTLY
> x2$xt
[1] 24.610 24.605 24.610 24.605 24.610
> h2$xt
[1] 0.003891051 0.003875910 0.003860829 0.003845806 0.003830842
> convolve(h2$xt,x2$xt)
[1] 0.4750436 0.4750438 0.4750435 0.4750437 0.4750435

# USING INVERSE FOURIER TRANSFORM
> f=fft(fft(h2$xt)*fft(x2$xt), inv=TRUE)
> Re(f)/length(f) 
[1] 0.4750438 0.4750435 0.4750437 0.4750435 0.4750436
>

Lets take the index 0. At 0, the convolution should simply be the last value of x2$xt (24.610) multiplied by first value of h开发者_高级运维2$xt (0.003891051) which should give convolution at index 0 = 24.610*0.003891051 = 0.09575877 which is way off from 0.4750436.

Am I doing something wrong? Why is the values so different from expected?

Both convolve and fft are circular. The first element of convolution must be the dot product of these two series. The results you obtain are correct in this sense.

To perform a linear convolution use:

convolve(h2$xt,x2$xt,type="open")

Circular convolution is also applied in this case but a required amount of zeros are padded to inputs to achieve linear convolution.

I believe there is not a direct way to achieve linear convolution with fft in R. However this doesn't really matter beacuse convolve itself uses the FFT approach you posted.

Circular Convolution

A discrete signal x is periodic if there is a period N such that x[n] = x[n+N] for all n. Such signals can be represented by N samples from x[0] to x[N-1].

... x[-2] x[-1] x[0] x[1] x[2] ... x[N-2] x[N-1] x[N] x[N+1] ...
                ^    this part is sufficient   ^

A regular definition of convolution between aperiodic x and y is defined as:

(x * y)[n] = sum{k in [-inf, inf]}(x[k]y[n-k])

However, for periodic signals, this formula does not produce finite results. To overcome this problem we define the circular convolution between periodic x and y.

(x * y)[n] = sum{k in [0, N-1]}(x[i]y[n-k])

When these two signals are represented with N values only, we can use y[n-k+N] in place of y[n-k] for negative values of n-k.

The cool thing with circular convolution is that it can calculate the linear convolution between box signals, which are discrete signals that have a finite number of non-zero elements.

Box signals of length N can be fed to circular convolution with 2N periodicity, N for original samples and N zeros padded at the end. The result will be a circular convolution with 2N samples with 2N-1 for linear convolution and an extra zero.

Circular convolution is generally faster than a direct linear convolution implementation, because it can utilize the Fast Fourier Transform, a fast algorithm to calculate the Discrete Fourier Transform, which is only defined for periodic discrete signals.

Please see:

• http://www.centerspace.net/blog/nmath/convolution-correlation-and-the-fft/

Also see:

• http://en.wikipedia.org/wiki/Circular_convolution
• http://en.wikipedia.org/wiki/Discrete_Fourier_transform