Last chance! 50% off unlimited learning
Sale ends in
One-dimensional second-order (biquadratic) sections IIR digital filtering.
sosfilt(sos, x, zi = NULL)The filtered signal, of the same dimensions as the input signal. In
case the zi input argument was specified, a list with two elements
is returned containing the variables y, which represents the output
signal, and zf, which contains the final state vector or matrix.
Second-order section representation, specified as an nrow-by-6
matrix, whose rows contain the numerator and denominator coefficients of
the second-order sections:
sos <- rbind(cbind(B1, A1), cbind(...),
cbind(Bn, An)), where B1 <- c(b0, b1, b2), and A1 <- c(a0,
a1, a2) for section 1, etc. The b0 entry must be nonzero for each section.
the input signal to be filtered, specified as a numeric or complex
vector or matrix. If x is a matrix, each column is filtered.
If zi is provided, it is taken as the initial state of the
system and the final state is returned as zf. If x is a vector,
zi must be a matrix with nrow(sos) rows and 2 columns. If
x is a matrix, then zi must be a 3-dimensional array of size
(nrow(sos), 2, ncol(x)). Alternatively, zi may be the
character string "zf", which specifies to return the final state
vector even though the initial state vector is set to all zeros. Default:
NULL.
Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.
The filter function is implemented as a series of second-order filters with direct-form II transposed structure. It is designed to minimize numerical precision errors for high-order filters [1].
Smith III, J.O. (2012). Introduction to digital filters, with audio applications (3rd Ed.). W3K Publishing.
filter, filtfilt, Sos
fs <- 1000
t <- seq(0, 1, 1/fs)
s <- sin(2* pi * t * 6)
x <- s + rnorm(length(t))
plot(t, x, type = "l", col="light gray")
lines(t, s, col="black")
sosg <- butter(3, 0.02, output = "Sos")
sos <- sosg$sos
sos[1, 1:3] <- sos[1, 1:3] * sosg$g
y <- sosfilt(matrix(sos, ncol=6), x)
lines(t, y, col="red")
## using 'filter' will handle the gain for you
y2 <- filter(sosg, x)
all.equal(y, y2)
## The following example is from Python scipy.signal.sosfilt
## It shows the instability that results from trying to do a
## 13th-order filter in a single stage (the numerical error
## pushes some poles outside of the unit circle)
arma <- ellip(13, 0.009, 80, 0.05, output='Arma')
sos <- ellip(13, 0.009, 80, 0.05, output='Sos')
x <- rep(0, 700); x[1] <- 1
y_arma <- filter(arma, x)
y_sos <- filter(sos, x)
plot(y_arma, type ="l")
lines (y_sos, col = 2)
legend("topleft", legend = c("Arma", "Sos"), lty = 1, col = 1:2)
Run the code above in your browser using DataLab