public
Created

An example showing the benefits of FFT convolution vs. manual calculation

  • Download Gist
fft_convolution.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
require(plyr)
set.seed(12345)
 
n = 10
n.sum = 2
 
# Simulate data
a = sample.int(10, n, replace=T)
df = data.frame(n=1:n, a)
 
# Precompute n-a
n.minus.a = with(df, n - a)
 
# Define a kernel for our sliding window
k = rep(0, n)
k[1:n.sum] = 1
 
# Instead of calculating the result of the sliding window manually, we can do
# the convolution quickly with fft
myConv <- function(x, k){
Fx = fft(x)
Fk = fft(k)
Fxk = Fx * Fk
xk = fft(Fxk, inverse=T)
(Re(xk) / n)[-(1:(n.sum-1))]
}
 
myConv(n.minus.a, k)
 
# Which is what happens under the hood when we do:
convolve(n.minus.a, k)[1:(length(n.minus.a)-n.sum+1)]
 
# Manual method
sliding = function(df, n, f)
ldply(1:(nrow(df) - n + 1), function(k)
f(df[k:(k + n - 1), ])
)
 
sliding(df, 2, function(df) with(df, data.frame(n = n[1], a = a[1], b = sum(n - a))))
 
# Reset n to something bigger like 10^4
system.time(myConv(n.minus.a, k))
system.time(convolve(n.minus.a, k, type='circ')[1:(length(n.minus.a)-n.sum+1)])
system.time(sliding(df, 2, function(df) with(df, data.frame(n = n[1], a = a[1], b = sum(n - a)))))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.