matchedfilter {bspec}R Documentation

Filter a noisy time series for a signal of given shape


Computes the maximized likelihood ratio (as a test- or detection-statistic) of "signal" vs. "noise only" hypotheses. The signal is modelled as a linear combination of orthogonal basis vectors of unknown amplitude and arrival time. The noise is modelled either as Gaussian or Student-t-distributed, and coloured.


matchedfilter(data, signal, noisePSD, timerange = NA,
    reconstruct = TRUE, two.sided = FALSE)

studenttfilter(data, signal, noisePSD, df = 10, timerange = NA,
    deltamax = 1e-06, itermax = 100,
    reconstruct = TRUE, two.sided = FALSE) 



the data to be filtered, a time series (ts) object.


the signal waveform to be filtered for. May be a vector, a matrix, a time series (ts) or a multivariate time series (mts) object.


the noise power spectral density. May be a vector of appropriate length (length(data)/2+1) or a function of frequency.


the number of degrees-of-freedom (nu[j]) for each frequency bin. May be a vector of appropriate length (length(data)/2+1) or a function of frequency.


the range of times (with respect to the data argument's time scale) to maximize the likelihood ratio over.


the minimal difference in logarithmic likelihoods to be aimed for in the EM-iterations (termination criterion).


the maximum number of EM-iterations to be performed.


a logical flag indicating whether the output is supposed to include the best-fitting signal waveform.


a logical flag indicating whether the noisePSD argument is to be interpreted as a one-sided or a two-sided spectrum.


The time series data d(t) is modelled as a superposition of signal s(beta,t0,t) and noise n(t):


where the signal is a linear combination of orthogonal (!) basis vectors s[i](t), and whose time-of-arrival is given by the parameter t0:

s(beta,t-t0) = beta[1] s[1](t-t_0) + ... + beta[k] s[k](t-t_0).

The noise is modelled as either Gaussian (matchedfilter()) or Student-t distributed (studenttfilter()) with given power spectral density and, for the latter model only, degrees-of-freedom parameters.

The filtering functions perform a likelihood maximization over the time-of-arrival t0 and coefficients beta. In the Gaussian model, the conditional likelihood, conditional on t0, can be maximized analytically, while the maximization over t0 is done numerically via a brute-force search. In the Student-t model, likelihood maximization is implemented using an EM-algorithm. The maximization over t0 is restricted to the range specified via the timerange argument.

What is returned is the maximized (logarithmic) likelihood ratio of "signal" versus "noise-only" hypotheses (the result's $maxLLR component), and the corresponding ML-estimates \code{tHat} and \code{betaHat}, as well as the ML-fitted signal ("$reconstruction").


A list containing the following elements:


the maximized likelihood ratio of signal vs. noise only hypotheses.


the range of times to maximize the likelihood ratio over (see the timerange input argument).


the ML-estimated vector of coefficients.


the ML-estimated signal arrival time.


the ML-fitted signal (a time series (ts) object).


an object of class call giving the function call that generated the result.

elements only for the matchedfilter() function:


the time series of (conditionally) maximized likelihood ratio for each given time point (the profile likelihood).

elements only for the studenttfilter() function:


a matrix indicating the progress of the EM-fitting.


Christian Roever,


Roever, C.. A Student-t based filter for robust signal detection. Physical Review D, 84(12):122004, 2011. See also arXiv preprint 1109.0442.

Roever, C., Meyer, R., Christensen, N. Modelling coloured residual noise in gravitational-wave signal processing. Classical and Quantum Gravity, 28(1):015010, 2011. See also arXiv preprint 0804.3853.


# sample size and sampling resolution:
deltaT  <- 0.001
N       <- 1000

# For the coloured noise, use some AR(1) process;
# AR noise process parameters:
sigmaAR <- 1.0
phiAR   <- 0.9

# generate non-white noise
# (autoregressive AR(1) low-frequency noise):
noiseSample <- rnorm(10*N)
for (i in 2:length(noiseSample))
  noiseSample[i] <- phiAR*noiseSample[i-1] + noiseSample[i]
noiseSample <- ts(noiseSample, deltat=deltaT)

# estimate the noise spectrum:
PSDestimate <- welchPSD(noiseSample, seglength=1,

# show noise and noise PSD:
  plot(noiseSample, main="noise sample")
  plot(PSDestimate$freq, PSDestimate$pow, log="y", type="l",
       main="noise PSD", xlab="frequency", ylab="power")

# generate actual data:
noise <- rnorm(N)
for (i in 2:length(noise))
  noise[i] <- phiAR*noise[i-1] + noise[i]
noise <- ts(noise, start=0, deltat=deltaT)

# the "sine-Gaussian" signal to be injected into the noise:
t0    <- 0.6
phase <- 1.0
signal <- exp(-(time(noise)-t0)^2/(2*0.01^2)) * sin(2*pi*150*(time(noise)-t0)+phase)

t <- seq(-0.1, 0.1, by=deltaT)
# the signal's orthogonal (sine- and cosine-) basis waveforms:
signalSine   <- exp(-t^2/(2*0.01^2)) * sin(2*pi*150*t)
signalCosine <- exp(-t^2/(2*0.01^2)) * sin(2*pi*150*t+pi/2)
signalBasis <- ts(cbind("sine"=signalSine, "cosine"=signalCosine),
                  start=-0.1, deltat=deltaT)
plot(signalBasis[,1], col="red", main="the signal basis")
lines(signalBasis[,2], col="green")
# (the sine- and cosine- components allow to span
#  signals of arbitrary phase)
# Note that "signalBasis" may be shorter than "data",
# but must be of the same time resolution.

# compute the signal's signal-to-noise ration (SNR):
signalSnr <- snr(signal, psd=PSDestimate$pow)

# scale signal to SNR = 6:
rho <- 6
data <- noise + signal * (rho/signalSnr)
data <- data * tukeywindow(length(data))
# Note that the data has (and should have!)
# the same resolution, size, and windowing applied
# as in the PSD estimation step.

# compute filters:
f1 <- matchedfilter(data, signalBasis, PSDestimate$power)
f2 <- studenttfilter(data, signalBasis, PSDestimate$power)

# illustrate the results:
  plot(data, ylab="", main="data")
  lines(signal* (rho/signalSnr), col="green")
  legend(0,max(data),c("noise + signal","signal only"),
         lty="solid", col=c("black","green"), bg="white")

  plot(signal * (rho/signalSnr), xlim=c(0.55, 0.65), ylab="",
       main="original & recovered signals")
  lines(f1$reconstruction, col="red")
  lines(f2$reconstruction, col="blue")
  abline(v=c(f1$tHat,f2$tHat), col=c("red", "blue"), lty="dashed")
  legend(0.55, max(signal*(rho/signalSnr)),
         c("injected signal","best-fitting signal (Gaussian model)",
           "best-fitting signal (Student-t model)"),
         lty="solid", col=c("black","red","blue"), bg="white")

  plot(f1$maxLLRseries, type="n", ylim=c(0, f1$maxLLR),
       main="profile likelihood (Gaussian model)",
       ylab="maximized (log-) likelihood ratio")
  lines(f1$maxLLRseries, col="grey")
  lines(window(f1$maxLLRseries, start=f1$timerange[1], end=f1$timerange[2]))
  abline(v=f1$timerange, lty="dotted")
  lines(c(f1$tHat,f1$tHat,-1), c(0,f1$maxLLR,f1$maxLLR), col="red", lty="dashed")

[Package bspec version 1.5 Index]