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
Community
  • 1
  • 1
thelatemail
  • 91,185
  • 12
  • 128
  • 188
  • 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. – Matthew Plourde Jan 17 '13 at 06:14
  • 1
    The behaviour of `filter` is _fundamentally different_ for the two `method`s. IMO this is abominable package / function design: they should have different names. – jessexknight Apr 29 '20 at 14:38

4 Answers4

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
agstudy
  • 119,832
  • 17
  • 199
  • 261
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
Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
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
AdamO
  • 4,283
  • 1
  • 27
  • 39
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

Michael Dorner
  • 17,587
  • 13
  • 87
  • 117
user2006050
  • 151
  • 2
  • 4