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 signal vs. noise only hypotheses. |
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 profile likelihood). |
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))