Lowpass Filtering

Description

Creates lowpass filters and offers further functionalities around them. Lowpass filters are commonly used in ion channel recordings.

Details

The main function of this package is lowpassFilter which creates lowpass filters, currently only Bessel filters are supported. randomGeneration and randomGenerationMA allow to generate random numbers that are filtered, i.e. follow a model for ion channel recordings, see (Pein et al., 2018, 2020). getConvolution, getConvolutionJump, and getConvolutionPeak allow to compute the convolution of a signal with the kernel of a lowpass filter.

References

Pein, F., Bartsch, A., Steinem, C., and Munk, A. (2020) Heterogeneous idealization of ion channel recordings - Open channel noise. Submitted.

Pein, F., Tecuapetla-Gómez, I., Schütte, O., Steinem, C., Munk, A. (2018) Fully-automatic multiresolution idealization for filtered ion channel recordings: flickering event detection. IEEE Trans. Nanobioscience, 17(3):300-320.

Pein, F. (2017) Heterogeneous Multiscale Change-Point Inference and its Application to Ion Channel Recordings. PhD thesis, Georg-August-Universität Göttingen. http://hdl.handle.net/11858/00-1735-0000-002E-E34A-7.

Hotz, T., Schütte, O., Sieling, H., Polupanow, T., Diederichsen, U., Steinem, C., and Munk, A. (2013) Idealizing ion channel recordings by a jump segmentation multiresolution filter. IEEE Trans. Nanobioscience, 12(4):376-386.

See Also

lowpassFilter, randomGeneration, randomGenerationMA, getConvolution, getConvolutionJump, getConvolutionPeak

Examples

# creates a lowpass filter
filter <- lowpassFilter(type = "bessel", param = list(pole = 4, cutoff = 0.1), sr = 1e4)
time <- 1:4000 / filter$sr

# creates a piecewise constant signal with a single peak
stepfun <- getSignalPeak(time, cp1 = 0.2, cp2 = 0.2 + 3 / filter$sr, 
                         value = 20, leftValue = 40, rightValue = 40)

# computes the convolution of the signal with the kernel of the lowpass filter
signal <- getConvolutionPeak(time, cp1 = 0.2, cp2 = 0.2 + 3 / filter$sr, 
                             value = 20, leftValue = 40, rightValue = 40,
                             filter = filter)

# generates random numbers that are filtered                             
data <- randomGenerationMA(n = 4000, filter = filter, signal = signal, noise = 1.4)

# generated data
plot(time, data, pch = 16)

# zoom into the single peak
plot(time, data, pch = 16, xlim = c(0.199, 0.202), ylim = c(19, 45))
lines(time, stepfun, col = "blue", type = "s", lwd = 2)
lines(time, signal, col = "red", lwd = 2)

# use of data randomGeneration instead
data <- randomGeneration(n = 4000, filter = filter, signal = signal, noise = 1.4)

# similar result
plot(time, data, pch = 16, xlim = c(0.199, 0.202), ylim = c(19, 45))
lines(time, stepfun, col = "blue", type = "s", lwd = 2)
lines(time, signal, col = "red", lwd = 2)

Time discrete convolution

Description

For developers only; computes a time discrete convolution.

Usage

.convolve(val, kern)

Arguments

Value

A numeric vector giving the convolution.

See Also

lowpassFilter

Deconvolution of a single jump / isolated peak

Description

For developers only; computes the deconvolution of a single jump or an isolated peak assuming that the observations are lowpass filtered. More details are given in (Pein et al., 2018).

Usage

.deconvolveJump(grid, observations, time, leftValue, rightValue,
                typeFilter, inputFilter, covariances) 
.deconvolvePeak(gridLeft, gridRight, observations, time, leftValue, rightValue,
                typeFilter, inputFilter, covariances, tolerance)

Arguments

Value

For .deconvolveJump a single numeric giving the jump point. For .deconvolvePeak a list containing the entries left, right and value giving the left and right jump point and the value of the peak, respectively.

References

Pein, F., Tecuapetla-Gómez, I., Schütte, O., Steinem, C., Munk, A. (2018) Fully-automatic multiresolution idealization for filtered ion channel recordings: flickering event detection. IEEE Trans. Nanobioscience, 17(3):300-320.

See Also

lowpassFilter

Convolved piecewise constant signals

Description

Creates piecewise constant signals with a single jump / peak. Computes the convolution of piecewise constant signals with the kernel of a lowpass filter.

Usage

getConvolution(t, stepfun, filter, truncated = TRUE)
getSignalJump(t, cp, leftValue, rightValue)
getConvolutionJump(t, cp, leftValue, rightValue, filter, truncated = TRUE)
getSignalPeak(t, cp1, cp2, value, leftValue, rightValue) 
getConvolutionPeak(t, cp1, cp2, value, leftValue, rightValue, filter, truncated = TRUE) 

Arguments

Value

a numeric of length length(t) giving the signal / convolution at time points t

References

Pein, F., Bartsch, A., Steinem, C., and Munk, A. (2020) Heterogeneous idealization of ion channel recordings - Open channel noise. Submitted.

Pein, F., Tecuapetla-Gómez, I., Schütte, O., Steinem, C., Munk, A. (2018) Fully-automatic multiresolution idealization for filtered ion channel recordings: flickering event detection. IEEE Trans. Nanobioscience, 17(3):300-320.

Pein, F. (2017) Heterogeneous Multiscale Change-Point Inference and its Application to Ion Channel Recordings. PhD thesis, Georg-August-Universität Göttingen. http://hdl.handle.net/11858/00-1735-0000-002E-E34A-7.

See Also

lowpassFilter

Examples

# creating and plotting a signal with a single jump at 0 from 0 to 1
time <- seq(-2, 13, 0.01)
signal <- getSignalJump(time, 0, 0, 1)
plot(time, signal, type = "l")

# setting up the filter
filter <- lowpassFilter(param = list(pole = 4, cutoff = 0.1))

# convolution with the truncated filter
convolution <- getConvolutionJump(time, 0, 0, 1, filter)
lines(time, convolution, col = "red")

# without truncating the filter, looks almost equal
convolution <- getConvolutionJump(time, 0, 0, 1, filter, truncated = FALSE)
lines(time, convolution, col = "blue")


# creating and plotting a signal with a single peak with jumps
# at 0 and at 3 from 0 to 1 to 0
time <- seq(-2, 16, 0.01)
signal <- getSignalPeak(time, 0, 3, 1, 0, 0)
plot(time, signal, type = "l")

# convolution with the truncated filter
convolution <- getConvolutionPeak(time, 0, 3, 1, 0, 0, filter)
lines(time, convolution, col = "red")

# without truncating the filter, looks almost equal
convolution <- getConvolutionPeak(time, 0, 3, 1, 0, 0, filter, truncated = FALSE)
lines(time, convolution, col = "blue")


# doing the same with getConvolution
# signal can also be an object of class stepblock instead,
# e.g. constructed by stepR::stepblock
signal <- data.frame(value = c(0, 1, 0), leftEnd = c(-2, 0, 3), rightEnd = c(0, 3, 16))

convolution <- getConvolution(time, signal, filter)
lines(time, convolution, col = "red")

convolution <- getConvolution(time, signal, filter, truncated = FALSE)
lines(time, convolution, col = "blue")


# more complicated signal
time <- seq(-2, 21, 0.01)
signal <- data.frame(value = c(0, 10, 0, 50, 0), leftEnd = c(-2, 0, 3, 6, 8),
                     rightEnd = c(0, 3, 6, 8, 21))

convolution <- getConvolution(time, signal, filter)
plot(time, convolution, col = "red", type = "l")

convolution <- getConvolution(time, signal, filter, truncated = FALSE)
lines(time, convolution, col = "blue")

Lowpass filtering

Description

Creates a lowpass filter.

Usage

lowpassFilter(type = c("bessel"), param, sr = 1, len = NULL, shift = 0.5)
## S3 method for class 'lowpassFilter'
print(x, ...)

Arguments

Value

An object of class lowpassFilter, i.e. a list that contains

"type", "param", "sr", "len"

the corresponding arguments

"kernfun"

the kernel function of the filter, obtained as the Laplace transform of the corresponding transfer function

"stepfun"

the step-response of the filter, i.e. the antiderivative of the filter kernel

"acfun"

the autocorrelation function, i.e. the convolution of the filter kernel with itself

"acAntiderivative"

the antiderivative of the autocorrelation function

"truncatedKernfun"

the kernel function of the at len / sr truncated filter, i.e. kernfun truncated and rescaled such that the new kernel still integrates to 1

"truncatedStepfun"

the step-response of the at len / sr truncated filter, i.e. the antiderivative of the kernel of the truncated filter

"truncatedAcfun"

the autocorrelation function of the at len / sr truncated filter, i.e. the convolution of the kernel of the truncated filter with itself

"truncatedAcAntiderivative"

the antiderivative of the autocorrelation function of the at len / sr truncated filter

"kern"

the digitised filter kernel normalised to one, i.e. kernfun((0:len + shift) / sr) / sum(kernfun((0:len + shift) / sr))

"step"

the digitised step-response of the filter, i.e. stepfun((0:len + shift) / sr)

"acf"

the discrete autocorrelation, i.e. acfun(0:len / sr)

"jump"

the last index of the left half of the filter, i.e. min(which(ret$step >= 0.5)) - 1L, it indicates how much a jump is shifted in time by a convolution of the signal with the digitised kernel of the lowpassfilter; if all values are below 0.5, len is returned with a warning

"number"

for developers; an integer indicating the type of the filter

"list"

for developers; a list containing precomputed quantities to recreate the filter in C++

References

Pein, F., Bartsch, A., Steinem, C., and Munk, A. (2020) Heterogeneous idealization of ion channel recordings - Open channel noise. Submitted.

Pein, F., Tecuapetla-Gómez, I., Schütte, O., Steinem, C., Munk, A. (2018) Fully-automatic multiresolution idealization for filtered ion channel recordings: flickering event detection. IEEE Trans. Nanobioscience, 17(3):300-320.

Pein, F. (2017) Heterogeneous Multiscale Change-Point Inference and its Application to Ion Channel Recordings. PhD thesis, Georg-August-Universität Göttingen. http://hdl.handle.net/11858/00-1735-0000-002E-E34A-7.

Hotz, T., Schütte, O., Sieling, H., Polupanow, T., Diederichsen, U., Steinem, C., and Munk, A. (2013) Idealizing ion channel recordings by a jump segmentation multiresolution filter. IEEE Trans. Nanobioscience, 12(4):376-386.

See Also

filter

Examples

filter <- lowpassFilter(type = "bessel", param = list(pole = 4L, cutoff = 1e3 / 1e4),
                        sr = 1e4)

# filter kernel, truncated version
plot(filter$kernfun, xlim = c(0, 20 / filter$sr))
t <- seq(0, 20 / filter$sr, 0.01 / filter$sr)
# truncated version looks very similar
lines(t, filter$truncatedKernfun(t), col = "red")

# filter$len (== 11) is chosen automatically
# this ensures that filter$acf < 1e-3 for this lag and at all larger lags
plot(filter$acfun, xlim = c(0, 20 / filter$sr), ylim = c(-0.003, 0.003))
abline(h = 0.001, lty = "22")
abline(h = -0.001, lty = "22")

abline(v = (filter$len - 1L) / filter$sr, col = "grey")
abline(v = filter$len / filter$sr, col = "red")

# filter with sr == 1
filter <- lowpassFilter(type = "bessel", param = list(pole = 4L, cutoff = 1e3 / 1e4))

# filter kernel and its truncated version
plot(filter$kernfun, xlim = c(0, 20 / filter$sr))
t <- seq(0, 20 / filter$sr, 0.01 / filter$sr)
# truncated version looks very similar
lines(t, filter$truncatedKernfun(t), col = "red")
# digitised filter
points((0:filter$len + 0.5) / filter$sr, filter$kern, col = "red", pch = 16)

# without a shift
filter <- lowpassFilter(type = "bessel", param = list(pole = 4L, cutoff = 1e3 / 1e4),
                        shift = 0)
# filter$kern starts with zero
points(0:filter$len / filter$sr, filter$kern, col = "blue", pch = 16)

# much shorter filter
filter <- lowpassFilter(type = "bessel", param = list(pole = 4L, cutoff = 1e3 / 1e4),
                        len = 4L)
points((0:filter$len + 0.5) / filter$sr, filter$kern, col = "darkgreen", pch = 16)

Random number generation

Description

Generate random numbers that are filtered. Both, signal and noise, are convolved with the given lowpass filter, see details. Can be used to generate synthetic data resembling ion channel recordings, please see (Pein et al., 2018, 2020) for the exact models.

Usage

randomGeneration(n, filter, signal = 0, noise = 1, oversampling = 100L, seed = n,
                 startTime = 0, truncated = TRUE)
randomGenerationMA(n, filter, signal = 0, noise = 1, seed = n,
                   startTime = 0, truncated = TRUE)                 

Arguments

Details

As discussed in (Pein et al., 2018) and (Pein et al., 2020), in ion channel recordings the recorded data points can be modelled as equidistant sampled at rate filter$sr from the convolution of a piecewise constant signal perturbed by Gaussian white noise scaled by the noise level with the kernel of an analogue lowpass filter. The noise level is either constant (homogeneous noise, see (Pein et al., 2018)) or itself varying (heterogeneous noise, see (Pein et al., 2020)). randomGeneration and randomGenerationMA generate synthetic data from such models. randomGeneration allows homogeneous and heterogeneous noise, while randomGenerationMA only allows homogeneous noise, i.e. noise has to be a single numeric giving the constant noise level. The resulting observations represent the conductance at time points startTime + 1:n / filter$sr.

The generated observations are the sum of a convolved signal evaluated at those time points plus centred Gaussian errors that are correlated (coloured noise), because of the filtering, and scaled by the noise level. The convolved signal evaluated at those time points can either by specified in signal directly or signal can specify a piecewise constant signal that will be convolved with the filter using getConvolution and evaluated at those time points. randomGenerationMA computes a moving average process with the desired autocorrelation to generate random errors. randomGeneration oversamples the error, i.e. generates errors at time points startTime + (seq(1 - filter$len + 1 / oversampling, n, 1 / oversampling) - 1 / 2 / oversampling) / filter$sr, which will then be convolved with the filter. For this function noise can either give the noise levels at those oversampled time points or specify a piecewise constant function that will be automatically evaluated at those time points.

Value

a numeric vector of length n giving the generated random numbers

References

Pein, F., Bartsch, A., Steinem, C., and Munk, A. (2020) Heterogeneous idealization of ion channel recordings - Open channel noise. Submitted.

Pein, F., Tecuapetla-Gómez, I., Schütte, O., Steinem, C., Munk, A. (2018) Fully-automatic multiresolution idealization for filtered ion channel recordings: flickering event detection. IEEE Trans. Nanobioscience, 17(3):300-320.

Pein, F. (2017) Heterogeneous Multiscale Change-Point Inference and its Application to Ion Channel Recordings. PhD thesis, Georg-August-Universität Göttingen. http://hdl.handle.net/11858/00-1735-0000-002E-E34A-7.

See Also

lowpassFilter, getConvolution

Examples

filter <- lowpassFilter(type = "bessel", param = list(pole = 4, cutoff = 0.1), sr = 1e4)
time <- 1:4000 / filter$sr
stepfun <- getSignalPeak(time, cp1 = 0.2, cp2 = 0.2 + 3 / filter$sr, 
                         value = 20, leftValue = 40, rightValue = 40)
signal <- getConvolutionPeak(time, cp1 = 0.2, cp2 = 0.2 + 3 / filter$sr, 
                             value = 20, leftValue = 40, rightValue = 40, filter = filter)
data <- randomGenerationMA(n = 4000, filter = filter, signal = signal, noise = 1.4)

# generated data
plot(time, data, pch = 16)

# zoom into the single peak
plot(time, data, pch = 16, xlim = c(0.199, 0.202), ylim = c(19, 45))
lines(time, stepfun, col = "blue", type = "s", lwd = 2)
lines(time, signal, col = "red", lwd = 2)

# use of randomGeneration instead
data <- randomGeneration(n = 4000, filter = filter, signal = signal, noise = 1.4)

# similar result
plot(time, data, pch = 16, xlim = c(0.199, 0.202), ylim = c(19, 45))
lines(time, stepfun, col = "blue", type = "s", lwd = 2)
lines(time, signal, col = "red", lwd = 2)

## heterogeneous noise
# manual creation of an object of class 'stepblock'
# instead the function stepblock in the package stepR can be used
noise <- data.frame(leftEnd = c(0, 0.2, 0.2 + 3 / filter$sr),
                    rightEnd = c(0.2, 0.2 + 3 / filter$sr, 0.4),
                    value = c(1, 30, 1))
attr(noise, "x0") <- 0
class(noise) <- c("stepblock", class(noise))

data <- randomGeneration(n = 4000, filter = filter, signal = signal, noise = noise)

plot(time, data, pch = 16, xlim = c(0.199, 0.202), ylim = c(19, 45))
lines(time, stepfun, col = "blue", type = "s", lwd = 2)
lines(time, signal, col = "red", lwd = 2)