rollconv: Windowed Convolution Filtering

In AkselA/R-rollfun:

Description

A windowed rolling mean function based on optimized FFT convolution filtering

Usage

rollconv(x, w, B = 7, scale.window = TRUE, fill = NA, partial = FALSE)

Arguments

x	a numeric vector
w	integer or numeric vector; the convolution window. Can either be a single integer specifying the width of a rectangular window, or a vector describing a discrete window of any shape
B	prime; the smoothness of the convolution sequences. Default is 7, meaning that the sequences will be padded to a length that has prime factors 1, 2, 3, 5 and 7
scale.window	logical; should the window be scaled so that its values sum to 1?
fill	single character or numeric; used for filling out start/end. default is NA. set to NULL for no fill

Details

This convolution filtering relies on the convolution theorem and the
Cooley–Tukey FFT algorithm to ensure efficient computation.
The FFT is fastest when the length of the series being transformed is smooth
(i.e., has many factors). If this is not the case, the transform
may take a long time to compute and will use a large amount of memory.
There is a trade-off between smoothness and series length. A smoothness of 2 will
make calculations more efficient, but will also on average require the series
to be padded more, making it longer, and hence require more resources.
A smoothness of 5 or 7 is generally OK.

Examples

rollconv(rep(0:1, 3, each=4), 5, fill=NULL)


x <- rep(0:1, 2, each=6)
plot(x)
matlines(cbind(
  rollconv(x, w=c(1, 2, 3, 4, 3, 2, 1)),
  rollconv(x, w=7),
  rollconv(x, w=c(1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1)),
  rollconv(x, w=11)
  ), lty=1, lwd=2)


### filtering out harmonics
x <- rep(c(1, 2, 4, 2), 10, each=3)
x <- x + (seq_along(x)/50) * sin(seq_along(x)/20)

# a triangular filter will in general produce a result with less
# ringing artefacts than a simple box/rectangular filter.

# Here the 'steppyness' is filtered out using a triangular filter
# of about twice the width of the steps.
r1 <- rollconv(x, w=c(1, 2, 3, 2, 1), fill=NA)

# using a rectangular filter retains the triangular shape of the wave
r2 <- rollconv(x, w=3, fill=NA)

# A simple box filter will isolate the underlying smooth wave as long as
# the width is carefully tuned
r3 <- rollconv(x, w=12, fill=NA)
r4 <- rollconv(x, w=13, fill=NA)

plot(x, pch=16, cex=0.4)
lines(r2, col="skyblue", lwd=2)
lines(r1, col="blue", lwd=2)
lines(r4, col="green", lwd=2)
lines(r3, col="darkgreen", lwd=2)


opar <- par(no.readonly=TRUE)
par(mar=c(2, 2, 1, 1), xaxs="r")

set.seed(1)
x <- rollconv(rnorm(400), 5, fill=NULL)

w <- 60
win1 <- trapezwin(w, l.slopes=1)
win2 <- trapezwin(w, l.slopes=3)
win3 <- trapezwin(w, l.slopes=5)
win4 <- trapezwin(w, l.slopes=7)

y1 <- rollconv(x, win1)
y2 <- rollconv(x, win2)
y3 <- rollconv(x, win3)
y4 <- rollconv(x, win4)

plot(x, type="l", col="grey", ylim=c(
  min(c(y1, y2, y3, y4), na.rm=TRUE)*1.1, 
  max(c(y1, y2, y3, y4), na.rm=TRUE)*1.1),
  lwd=0.5)
lines(y1, col="#FF000088", lwd=2)
lines(y2, col="#FFAA0088", lwd=2)
lines(y3, col="#00AAFF88", lwd=2)
lines(y4, col="#0000FF88", lwd=2)

par(opar)


# calculating partial results gets slow with wide windows
set.seed(1)
l <- 2e4
r <- log1p(rgamma(l, 1.2+sin(1:l/l*pi*4)))
plot(r, type="l", col="#00000044")
rr1 <- rollconv(r, w=trapezwin(1e4, 0.8), partial=TRUE)
lines(rr1, col=5, lwd=1)
rr1 <- rollconv(r, w=trapezwin(1e4, 0.8), partial=FALSE)
lines(rr1, col=4, lwd=3)

AkselA/R-rollfun documentation built on May 31, 2019, 8:41 a.m.