bspec (version 1.6)

matchedfilter: Filter a noisy time series for a signal of given shape

Description

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.

Usage

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)

Arguments

data

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

signal

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

noisePSD

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

df

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.

timerange

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

deltamax

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

itermax

the maximum number of EM-iterations to be performed.

reconstruct

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

two.sided

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

Value

A list containing the following elements:

maxLLR

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

timerange

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

betaHat

the ML-estimated vector of coefficients.

tHat

the ML-estimated signal arrival time.

reconstruction

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

call

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

elements only for the matchedfilter() function:

maxLLRseries

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

elements only for the studenttfilter() function:

EMprogress

a matrix indicating the progress of the EM-fitting.

Details

The time series data \(d(t)\) is modelled as a superposition of signal \(s(\beta,t_0,t)\) and noise \(n(t)\): $$d(t)=s(\beta,t-t_0)+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 \(t_0\): $$s(\beta,t-t_0)=\sum_{i=1}^k \beta_i s_i(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 \(t_0\) and coefficients \(\beta\). In the Gaussian model, the conditional likelihood, conditional on \(t_0\), can be maximized analytically, while the maximization over \(t_0\) is done numerically via a brute-force search. In the Student-t model, likelihood maximization is implemented using an EM-algorithm. The maximization over \(t_0\) 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 \(\hat{t}_0\) and \(\hat{\beta}\), as well as the ML-fitted signal ("$reconstruction").

References

Roever, C. A Student-t based filter for robust signal detection. Physical Review D, 84(12):122004, 2011. 10.1103/PhysRevD.84.122004. 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. 10.1088/0264-9381/28/1/015010. See also arXiv preprint 0804.3853.

Examples

# NOT RUN {
# 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,
                        windowingPsdCorrection=FALSE)

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

# 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)
plot(signal)

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:
par(mfrow=c(3,1))
  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")
par(mfrow=c(1,1))
# }