deconvolve: Deconvolve a known linear filter

In jrevenaugh/TSAUMN: Time Series Analysis Tools Used in ESCI5201 at the University of Minnesota

Description

Deconvolve a known linear filter from sampled univariate or multivariate timeseries using standard spectral division with or without some regularization.

Usage

deconvolve(s, f, type = c("wl", "dls", "noreg"), alpha = 0.1,
  metric = c("median", "max", "percentile"))

Arguments

s	real-valued univariate or multivariate (as matrix columns) timeseries.
f	real-valued filter to deconvolve.
type	one of water-level ("wl", default), damped "least squqres" ("dls") or no regularization ("noreg"). See detail for more information
alpha	multiplier of the max or median value of the denominator for the corresponding metric, If metric is "percentile", this specifies the percentile value of the denominator to use in regularization.
metric	one of "max" (default), "median" or "percentile." Together with alpha this determines the regularization of the deconvolution.

Details

If the convolved filter has spectral zeros (or near zeros), any noise in that band will be greatly exaggerated by deconvolution. To control this, it is common to regularize the denominator during spectral division (deconvolution in the frequency domain). Two methods are provided. (1) water-level ("wl") doesn't allow the value of the denominator to drop below the scaled metric. Visually, the denominator plotted against frequency is filled with water up to the scaled metric, thus filling in "holes." (2) damped "least-squares" ("dls") simply adds the scaled metric to the denominator at all frequencies. Visually, this shifts the baseline of the denominator up. While this reduces the noise-amplification of zeros and near-zeros, it biases results at all frequencies. Type "noreg" applies no regularization and will produce useless results for noisy inputs.

The scaled metric is one of: alpha * max( D ), alpha * median( D ) or quantile( D, alpha ) where D is the denominator (function of frequency).

Deconvolution is circular. Long filters will cycle around. Padding s prior to deconvolution is often desirable.

Value

Deconvolved timeseries

Examples

# Create spike sequence
r <- runif( 1024, min = 0, max = 1 )
c <- which( r > 0.95 )
r <- rep( 0, 1024 )
r[c] <- 5 * rnorm( length( c ) )

# Create a wavelet
t <- seq( 0, 1023, 1 )
w <- 10 * dnorm( t, sd = 10 ) * sin( 2.0 * pi * t / 20 )

# Convolve spikes with wavelet and add white noise.
z <- convolve( r, w, conj = FALSE, type = "circular" ) + rnorm( 1024, sd = 0.1 )

# Deconvolve
rd <- deconvolve( z, w, type = "wl", metric = "max", alpha = 0.005 )
plotl( r, main = "Water Level" )
lines( rd, col = "red" )
rd <- deconvolve( z, w, type = "dls", metric = "max", alpha = 0.005 )
plotl( r, main = "Damped LS" )
lines( rd, col = "red" )

jrevenaugh/TSAUMN documentation built on Nov. 8, 2019, 2:20 p.m.