1
$\begingroup$

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).

$\endgroup$
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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