matchedfilter {bspec} | R Documentation |

## 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 ( |

`signal` |
the signal waveform to be filtered for. May be a vector,
a matrix, a time series ( |

`noisePSD` |
the noise power spectral density. May be a vector of
appropriate length ( |

`df` |
the number of degrees-of-freedom ( |

`timerange` |
the range of times (with respect to the |

`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 |

`two.sided` |
a |

### 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`

").

### Value

A list containing the following elements:

`maxLLR` |
the maximized likelihood ratio of |

`timerange` |
the range of times to maximize the likelihood ratio
over (see the |

`betaHat` |
the ML-estimated vector of coefficients. |

`tHat` |
the ML-estimated signal arrival time. |

`reconstruction` |
the ML-fitted signal (a time series
( |

`call` |
an object of class |

elements only for the `matchedfilter()`

function:

`maxLLRseries` |
the time series of (conditionally) maximized
likelihood ratio for each given time point (the |

elements only for the `studenttfilter()`

function:

`EMprogress` |
a |

### Author(s)

Christian Roever, christian.roever@med.uni-goettingen.de

### References

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

### Examples

```
# 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))
```

*bspec*version 1.6 Index]