# 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. ```r 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. ## Details The time series data $d(t)$ is modelled as a superposition of signal $s(beta,t0,t)$ and noise $n(t)$: $$ d(t)=s(\beta,t-t_0)+n(t),d(t)=s(beta,t-t0)+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-t_0)=\sum_{i=1}^k \beta_i s_i(t-t_0).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`"). ## Returns 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. ## References Roever, C. A Student-t based filter for robust signal detection. **Physical Review D**, 84 (12):122004, 2011. tools:::Rd_expr_doi("10.1103/PhysRevD.84.122004") . See also [arXiv preprint 1109.0442](https://arxiv.org/abs/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. tools:::Rd_expr_doi("10.1088/0264-9381/28/1/015010") . See also [arXiv preprint 0804.3853](https://arxiv.org/abs/0804.3853). ## Author(s) Christian Roever, christian.roever@med.uni-goettingen.de ## Examples ```r # 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)) ```

matchedfilter function