Instantly share code, notes, and snippets.

johncolby/fft_convolution.R Created Oct 27, 2011

What would you like to do?
An example showing the benefits of FFT convolution vs. manual calculation
 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)))))