sosfilt {gsignal}

Second-order sections filtering

Description

One-dimensional second-order (biquadratic) sections IIR digital filtering.

Usage

sosfilt(sos, x, zi = NULL)

Arguments

sos

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.

x

the input signal to be filtered, specified as a numeric or complex vector or matrix. If x is a matrix, each column is filtered.

zi

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.

Details

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

Value

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.

Author(s)

Geert van Boxtel, [e-mail].

References

Smith III, J.O. (2012). Introduction to digital filters, with audio applications (3rd Ed.). W3K Publishing.

See Also

filter, filtfilt, Sos

Examples

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)