27

I am seeking some simple (i.e. - no maths notation, long-form reproducible code) examples for the filter function in R I think I have my head around the convolution method, but am stuck at generalising the recursive option. I have read and battled with various documentation, but the help is just a bit opaque to me.

Here are the examples I have figured out so far:

# Set some values for filter components
f1 <- 1; f2 <- 1; f3 <- 1;

And on we go:

# basic convolution filter
filter(1:5,f1,method="convolution")
[1] 1 2 3 4 5

#equivalent to:
x[1] * f1 
x[2] * f1 
x[3] * f1 
x[4] * f1 
x[5] * f1 

# convolution with 2 coefficients in filter
filter(1:5,c(f1,f2),method="convolution")
[1]  3  5  7  9 NA

#equivalent to:
x[1] * f2 + x[2] * f1
x[2] * f2 + x[3] * f1
x[3] * f2 + x[4] * f1 
x[4] * f2 + x[5] * f1 
x[5] * f2 + x[6] * f1

# convolution with 3 coefficients in filter
filter(1:5,c(f1,f2,f3),method="convolution")
[1] NA  6  9 12 NA

#equivalent to:
 NA  * f3 + x[1] * f2 + x[2] * f1  #x[0] = doesn't exist/NA
x[1] * f3 + x[2] * f2 + x[3] * f1
x[2] * f3 + x[3] * f2 + x[4] * f1 
x[3] * f3 + x[4] * f2 + x[5] * f1 
x[4] * f3 + x[5] * f2 + x[6] * f1

Now's when I am hurting my poor little brain stem. I managed to figure out the most basic example using info at this post: https://stackoverflow.com/a/11552765/496803

filter(1:5, f1, method="recursive")
[1]  1  3  6 10 15

#equivalent to:

x[1]
x[2] + f1*x[1]
x[3] + f1*x[2] + f1^2*x[1]
x[4] + f1*x[3] + f1^2*x[2] + f1^3*x[1]
x[5] + f1*x[4] + f1^2*x[3] + f1^3*x[2] + f1^4*x[1]

Can someone provide similar code to what I have above for the convolution examples for the recursive version with filter = c(f1,f2) and filter = c(f1,f2,f3)?

Answers should match the results from the function:

filter(1:5, c(f1,f2), method="recursive")
[1]  1  3  7 14 26

filter(1:5, c(f1,f2,f3), method="recursive")
[1]  1  3  7 15 30

EDIT

To finalise using @agstudy's neat answer:

> filter(1:5, f1, method="recursive")
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]  1  3  6 10 15
> y1 <- x[1]                                            
> y2 <- x[2] + f1*y1      
> y3 <- x[3] + f1*y2 
> y4 <- x[4] + f1*y3 
> y5 <- x[5] + f1*y4 
> c(y1,y2,y3,y4,y5)
[1]  1  3  6 10 15

and...

> filter(1:5, c(f1,f2), method="recursive")
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]  1  3  7 14 26
> y1 <- x[1]                                            
> y2 <- x[2] + f1*y1      
> y3 <- x[3] + f1*y2 + f2*y1
> y4 <- x[4] + f1*y3 + f2*y2
> y5 <- x[5] + f1*y4 + f2*y3
> c(y1,y2,y3,y4,y5)
[1]  1  3  7 14 26

and...

> filter(1:5, c(f1,f2,f3), method="recursive")
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]  1  3  7 15 30
> y1 <- x[1]                                            
> y2 <- x[2] + f1*y1      
> y3 <- x[3] + f1*y2 + f2*y1
> y4 <- x[4] + f1*y3 + f2*y2 + f3*y1
> y5 <- x[5] + f1*y4 + f2*y3 + f3*y2
> c(y1,y2,y3,y4,y5)
[1]  1  3  7 15 30
2
  • 2
    Think of filter as stepping through your original vector, applying the weights and summing at each step. The recursive filter is just like the convolution filter, except the weights f1, ..., fn automatically become c(1, f1, ..., fn), and at each step 1 is applied to the current value, while f1, ..., fn are applied to the last n values from the new corrected vector being created, instead of the original values. With convolution, (with default sides = 2), the weights straddle the current value, with next n/2 original values on one side, and the previous n/2 original values on the other. Jan 17, 2013 at 6:14
  • 1
    The behaviour of filter is fundamentally different for the two methods. IMO this is abominable package / function design: they should have different names. Apr 29, 2020 at 14:38

4 Answers 4

23

In the recursive case, I think no need to expand the expression in terms of xi. The key with "recursive" is to express the right hand expression in terms of previous y's.

I prefer thinking in terms of filter size.

filter size =1

y1 <- x1                                            
y2 <- x2 + f1*y1      
y3 <- x3 + f1*y2 
y4 <- x4 + f1*y3 
y5 <- x5 + f1*y4 

filter size = 2

y1 <- x1                                            
y2 <- x2 + f1*y1      
y3 <- x3 + f1*y2 + f2*y1    # apply the filter for the past value and add current input
y4 <- x4 + f1*y3 + f2*y2
y5 <- x5 + f1*y4 + f2*y3
2
  • Well put: the key with "recursive" (in contrast to "convolution") is to express the RHS in terms of the y's. Jan 17, 2013 at 7:05
  • Thanks for this answer. You don't know how simple you've made a seemingly complicated concept! Jan 17, 2013 at 23:56
3

Here's the example that I've found most helpful in visualizing what recursive filtering is really doing:

(x <- rep(1, 10))
# [1] 1 1 1 1 1 1 1 1 1 1

as.vector(filter(x, c(1), method="recursive"))  ## Equivalent to cumsum()
#  [1]  1  2  3  4  5  6  7  8  9 10
as.vector(filter(x, c(0,1), method="recursive"))
#  [1] 1 1 2 2 3 3 4 4 5 5
as.vector(filter(x, c(0,0,1), method="recursive"))
#  [1] 1 1 1 2 2 2 3 3 3 4
as.vector(filter(x, c(0,0,0,1), method="recursive"))
#  [1] 1 1 1 1 2 2 2 2 3 3
as.vector(filter(x, c(0,0,0,0,1), method="recursive"))
#  [1] 1 1 1 1 1 2 2 2 2 2
2

With recursive, the sequence of your "filters" is the additive coefficient for the previous sums or output values of the sequence. With filter=c(1,1) you're saying "take the i-th component in my sequence x and add to it 1 times the result from the previous step and 1 times the results from the step before that one". Here's a couple examples to illustrate

I think the lagged effect notation looks like this:

## only one filter, so autoregressive cumsum only looks "one sequence behind"
> filter(1:5, c(2), method='recursive')
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]  1  4 11 26 57

1 = 1
2*1 + 2 = 4
2*(2*1 + 2) + 3 = 11
...

## filter with lag in it, looks two sequences back
> filter(1:5, c(0, 2), method='recursive')
Time Series:
Start = 1 
End = 5 
Frequency = 1 
[1]  1  2  5  8 15

1= 1
0*1 + 2 = 2
2*1 + 0*(0*1 + 2) + 3 = 5
2*(0*1 + 2) + 0 * (2*1 + 0*(0*1 + 2) + 3) + 4 = 8
2*(2*1 + 0*(0*1 + 2) + 3) + 0*(2*(0*1 + 2) + 0 * (2*1 + 0*(0*1 + 2) + 3) + 4) + 5 = 15

Do you see the cumulative pattern there? Put differently.

1 = 1
0*1 + 2 = 2
2*1 + 0*2 + 3 = 5
2*2 + 0*5 + 4 = 8
2*5 + 0*8 + 5 = 15
0

I spent one hour in reading this, below is my summary, by comparison with Matlab

NOTATION: command in Matlab = command in R

filter([1,1,1], 1, data) = filter(data, [1,1,1], method = "convolution") ; but the difference is that the first 2 elements are NA 


filter(1, [1,-1,-1,-1], data) = filter(data, [1,1,1], method = "recursive")

If you know some from DSP, then recursive is for IIR, convolution is for FIR

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.