# trying to reproduce the behavior of initial conditions for the filter delays of scipy lfilter with R

I'm trying to reproduce the behavior of lfilter and lfilter_zi (from scipy) in R (using signal package if possible).

This is the example provided by scipy.signal.lfilter (I just removed a random component to be able to compare it with R)

from scipy import signal

import matplotlib.pyplot as plt
import numpy as np

t = np.linspace(-1, 1, 201)

x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +

0.1*np.sin(2*np.pi*1.25*t + 1) +

0.18*np.cos(2*np.pi*3.85*t))

# Create an order 3 lowpass butterworth filter:

b, a = signal.butter(3, 0.05)

# Apply the filter to x. Use lfilter_zi to choose the initial condition of the filter:

zi = signal.lfilter_zi(b, a)

z, _ = signal.lfilter(b, a, x, zi=zi*x[0])


zi and z look like this:

zi
array([ 0.99958345, -1.68782358,  0.73058189])

z

array([-0.81143825, -0.81141815, -0.81127909, -0.81079006, -0.80958088,
-0.80715887, -0.80292746, -0.79620658, -0.78625483, -0.77229353,
-0.75353266, -0.72919883, -0.6985649 , -0.66098087, -0.61590557,
-0.5629379 , -0.50184682, -0.43259855, -0.35537974, -0.27061505,
-0.178978  , -0.08139367,  0.02096739,  0.12670479,  0.23421348,
0.34171818,  0.44731638,  0.54902848,  0.64485359,  0.7328288 ,
0.81108977,  0.87793009,  0.93185693,  0.97164066,  0.99635612,
1.00541372,  0.99857903,  0.97597976,  0.93809993,  0.88576138,
0.82009349,  0.74249248,  0.65457214,  0.55810822,  0.45497912,
0.34710537,  0.23639076,  0.12466757,  0.01364821, -0.09511478,
-0.20025897, -0.3006316 , -0.39529816, -0.48354019, -0.56484324,
-0.63887624, -0.70546404, -0.76455525, -0.8161876 , -0.86045305,
-0.89746513, -0.92733025, -0.95012489, -0.96587999, -0.97457311,
-0.97612896, -0.97042782, -0.95732115, -0.93665319, -0.90828687,
-0.87213219, -0.82817497, -0.77650384, -0.71733361, -0.65102302,
-0.57808558, -0.49919234, -0.41516614, -0.32696721, -0.23567056,
-0.14243637, -0.04847447,  0.04499495,  0.13678266,  0.2257712 ,
0.31095339,  0.39146601,  0.46661712,  0.53590537,  0.59903023,
0.65589258,  0.70658571,  0.75137705,  0.79068186,  0.82503022,
0.85502924,  0.88132248,  0.90454911,  0.92530473,  0.94410635,
0.96136324,  0.9773553 ,  0.99222013,  1.00594934,  1.01839434,
1.02928094,  1.03823213,  1.0447972 ,  1.04848578,  1.04880455,
1.04529431,  1.03756542,  1.02532925,  1.00842383,  0.98683223,
0.96069232,  0.9302976 ,  0.89608879,  0.85863664,  0.81861716,
0.77678038,  0.73391465,  0.69080854,  0.64821242,  0.60680207,
0.56714639,  0.5296811 ,  0.49468995,  0.46229459,  0.43245354,
0.40497058,  0.37951183,  0.35563067,  0.33279901,  0.31044314,
0.28798196,  0.26486556,  0.24061173,  0.21483836,  0.18728986,
0.15785602,  0.12658225,  0.0936706 ,  0.05947144,  0.02446644,
-0.01075634, -0.04553372, -0.07916315, -0.11094097, -0.14020089,
-0.16635045, -0.18890337, -0.20750588, -0.22195553, -0.23221132,
-0.23839469, -0.24078114, -0.23978314, -0.23592523, -0.22981275,
-0.22209607, -0.21343243, -0.20444753, -0.1956991 , -0.18764473,
-0.18061561, -0.17479785, -0.17022247, -0.16676453, -0.16415158,
-0.16198079, -0.15974396, -0.15685869, -0.15270414, -0.14665909,
-0.13814026, -0.12663849, -0.11175076, -0.09320618, -0.07088438,
-0.04482529, -0.01522967,  0.01754961,  0.05302499,  0.09059946,
0.12959802,  0.16930358,  0.20899523,  0.24798663,  0.28566232,
0.32150973,  0.35514518,  0.38633225,  0.41499146,  0.44120071,
0.46518647,  0.48730618,  0.50802294,  0.52787388,  0.54743416,
0.56727857])


This my attempt of R code for the same thing:

t <- seq(-1, 1, length.out =  201)

x <- (sin(2*pi*0.75*t*(1-t) + 2.1) +
0.1*sin(2*pi*1.25*t + 1) +
0.18*cos(2*pi*3.85*t))

# Create an order 3 lowpass butterworth filter:
h = signal::butter(3, 0.05)
b <- h$$b a <- h$$a

# does lfilter_zi:
n <- max(length(a), length(b))
B <- b[2:length(b)] - a[2:length(a)] * b[1]
A <- t(pracma::compan(a))
zi <- solve( diag(rep(1,n-1)) - A, B)

# doesn't do lfilter with initial conditions
z <- signal::filter(filt = b, a = a, x = x, init = zi * x[1])


Now all the intermediate steps give me the same results as python, except the last one that is off. But if I don't set initial conditions, scipy.signal.lfilter and filter give identical results. Any ideas?

(I've seen this https://stackoverflow.com/questions/27220671/matlab-filtfilt-function-implementation-in-java/52742099#52742099, but still doesn't help me to scale it to R) (And this Simulate butterworth filter with initial condition of state vector in R doesn't really solve my issue).

I think the problem might be the following.

The signal::filter() function is a direct-form I implementation of a digital filter. It first calls stats::filter() with the moving-average coefficients (b) and method = "convolution", and applies stats:filter() again on the output of the first run, with the autoregressive coefficients (a) and method = "recursive". But stats::filter() does not return the final state of the filter like the comparable functions in Matlab/Octave or Python. The initial conditions in signal::filter() are therefore used in a slightly different way, which can be useful in many cases, but does appear to work in your situation.

I developed a package called gsignal, which complements and extends signal. It can be installed in R using devtools::install_github("gjmvanboxtel/gsignal"). The gsignal::filter() function contained in that package is a direct-form II implementation of the filter, coded in C++, and handles initial conditions just as the Python and Matlab/Octave equivalents. It also contains a function filter_zi() as in Python scipy.signal.lfilter_zi (the same functionality could also be obtained by filtic(), by the way).

library(gsignal)
t <- seq(-1, 1, length.out =  201)
x <- (sin(2 * pi * 0.75 * t * (1 - t) + 2.1)
+ 0.1 * sin(2 * pi * 1.25 * t + 1)
+ 0.18 * cos(2 * pi * 3.85 * t))
h <- butter(3, 0.05)
zi <- filter_zi(h)
z1 <- filter(h, x)
z2 <- filter(h, x, zi * x[1])
plot(t, x, type = "l")
lines(t, z1, col = "red")
lines(t, z2\$y, col = "green")
legend("bottomright", legend = c("Original signal",
"Filtered without initial conditions",
"Filtered with initial conditions"),
lty = 1, col = c("black", "red", "green"))


Hope this helps,

Geert