Methods to extract signal from three very noisy time series of same event

Question

I have three time series of same length, all containing magnitude measurements of the same event "A". But each time series is using a different method of measurement.

My goal is to merge the three time series into a single one, so that I can more easily find the "upward" and "downward" phases. (Like a smooth sinus-curve or binary data-set 1=up, 0=down.)

For instance, as you can see the "raw" measurement data is very noisy and looks like this for measurement system number 3:

To get a better understanding of the data I have used a Moving Average (window length 2000) to smooth the time series data for each measurement system 1 - 3, which yields the following figure:

Since they are measurements of the same event, and often have their peaks and valleys at similar times, I would like to merge the three time series into a single time series, with as little noise as possible. What methods should I try? I'm open to try anything!

I have tried using Fourier analysis (FFT in the SciPy package), but I cannot find any significant frequencies in any part of the data.

[edit for Whuber's comment] I unfortunately only have access to the MA smoothed data for all three, which has the following statistics (using DataFrame.describe() function):

• MA_measurement 1: mean 0.991957, std 0.156941
• MA_measurement 2: mean -0.000003, std 0.000016
• MA_measurement 3: mean -0.000800, std 0.000856

And using DataFrame.corr() to get the correlations between the (smoothed) measurement systems:

$$\begin{array}{c|c|c|c|} & \text{MA_meas3} & \text{MA_meas2} & \text{MA_meas1} \\ \hline \text{MA_meas3} & 1.000000 & 0.337050 & 0.297922\\ \hline \text{MA_meas2} & 0.337050 & 1.000000 & 0.807282 \\ \hline \text{MA_meas1} & 0.297922 & 0.807282 & 1.000000 \\ \hline \end{array}$$

For the unsmoothed (i.e. "raw") data I only have data for time series number 2 and 3:

• measurement 2: mean -0.000003, std 0.000747
• measurement 3: mean -0.000812, std 0.022399

And they have a correlation of Corr(meas2,meas3)=0.027199.

[edit 2] I have been able to get hold of the the MA data shown in the second graph (i.e. 3 subplots), hope this can be of use!

MA data measurement 1 MA data measurement 2 MA data measurement 3

[edit 3] To elaborate on Matt F.'s comment: I hope to find around 33 upward and downward phases in total (peaks + troughs) in each MA measurement series (see edit 2) of event A. In theory it should be cyclical, i.e. up -> down -> up -> down... etc.

Answer 1

We have only few details given here, but in general, the problem of estimating the underlying signal from noisy measurements seems to me like a Kalman filtering problem (or structural time series more generally).

In the following, I follow very closely this textbook: https://www.stat.pitt.edu/stoffer/tsa4/ , Chapter 6 and in the code I'm using the R package accompanying the textbook.

(Here is how Kalman filter works in picutres: https://www.bzarg.com/p/how-a-kalman-filter-works-in-pictures/ )

To illustrate how one would approach such a problem, conisder $p=3$ noisy measurements $y_t$ of the same underlying phenomenon $x_t$. We can write this as

$y_t = A x_t + \epsilon_t$

$x_{t+1} = \delta + \Phi x_t + \eta_t$

$\epsilon_t \sim N(0,R)$, $\eta_t \sim N(0,Q)$, $x_1 \sim N(\mu_0,\Sigma_1)$

For simplicity, let's say that in our case $\Phi = 1$ and $A =\begin{bmatrix} 1 \\ 1 \\ 1\end{bmatrix}$, $\delta=0$, and that $R = \begin{bmatrix} r_{11} & 0 & 0 \\ 0 & r_{22} & 0 \\ 0 & 0 & r_{33} \\\end{bmatrix}$ is diagonal matrix (meaning that the three measurements $y_{t,1}, y_{t,2}, y_{t,3}$ are uncorrelated). The four unknown parameters that we want to estimate are variance of $\eta$ ,$Q$, and variances of $\epsilon$, $r_{ii}$. The following code shows how to do setup the above model and estimate it.

We can estimate this parameters by maximum likelihood, meaning the we will look for the values of these parameters for which the model-implied series of observations $\hat{y}_t$ resembles the actual one as much as possible.

---

---

We start by creating some artificial data

# library accompanying the book

library("astsa")

# for reproducibility
set.seed(1);

# dimensions of the problem in our case
N <- 400; # number of periods
p <- 3; # no. of observation series

  

# we assumed that this holds
A <- matrix(c(1, 1, 1), p, 1);
Phi <- 1;drift<-0;

# generate shocks
shocks_eps1 <- rnorm(N,mean=0,sd=1);
shocks_eps2 <- rnorm(N,mean=0,sd=2);
shocks_eps3 <- rnorm(N,mean=0,sd=3);
shocks_eta <- rnorm(N,mean=0,sd=.5);

shocks_eps=cbind(shocks_eps1,shocks_eps2,shocks_eps3)

# create artificial data
y_t <- matrix(0, N, p) # Stores observations
x_t <- matrix(0, N, 1) # stores unobserved true state
x_t[1, 1] <- 0;
y_t[1, ] <- A * x_t[1, 1] + shocks_eps[1, ];

for (i in 2:N) {

x_t[i, 1] <- drift + Phi*x_t[i - 1, 1] + shocks_eta[i - 1];
y_t[i, ] <- A * x_t[i, 1] + shocks_eps[i, ];

}

Next, let's estimate the four parameters as described above

# now setup the estimation procedure
# see the textbook Ch. 6 for details

input <- rep(1,N)
mu0 <- 0; Sigma0 <- 1;

A <- array(A, dim=c(p,1,N))

# Function to Calculate Likelihood
Linn = function(para){
cQ <- para[1] # Q
cR1 <- para[2] # 11 element of R
cR2 <- para[3] # 22 element of R
cR3 <- para[4] # 22 element of R
cR <- diag(c(cR1,cR2, cR3)) # put the matrix together
kf <- Kfilter1(N,y_t,A,mu0,Sigma0,Phi,0,drift,cQ,cR,input)
return(kf$like)
}

# Estimation
init.par <- c(.5,1,2,3) # initial values of parameters
(est <- optim(init.par, Linn, NULL, method='Nelder-Mead', hessian=TRUE,
    control=list(trace=1,REPORT=1))) # output not shown
(est <- optim(est$par, Linn, NULL, method='BFGS', hessian=TRUE,
    control=list(trace=1,REPORT=1))) # output not shown

cQ <- est$par[1]
cR1 <- est$par[2]
cR2 <- est$par[3]
cR3 <- est$par[4]
cR <- diag(c(cR1,cR2, cR3)) # put the matrix together

ks <- Ksmooth1(N,y_t,A,mu0,Sigma0,Phi,0,drift,cQ,cR,input)

# extract estimated signal
xsm <- ts(as.vector(ks$xs), start=1)

# and correspodning uncertainty
rmse <- ts(sqrt(as.vector(ks$Ps)), start=1)

Answer 2

I would apply PCA to three series, then assess how much variance is explained by the first principal component. If it is very high, then use the PC1 as your signal. The PCA transforms your three signals into three orthogonal linear combinations of signals in decreasing order of variance. So, your first principal component can be seen as an average level of three signals. If PCA works very well, then PC1 will explain the most variance in the data

You may need to mess with phase differences, depending on how the measurements were done. For instance, if it's three microphones picking up the source, then sound waves reach with a lag. This can easily be taken care of with coherence/phase shift analysis. Once you identify the phase shifts, it's still the same PCA but with appropriate lags

Answer 3

As a base reference you could try to directly merge the three data sets by interlacing them and then analyzing this new data set.

First you need to do a linear transformation on each data set individually so they have the same mean and variance (mean=0 and variance=1 is the simplest choice). Then you transform the x-scales of the three data sets by $x \mapsto 3x+1, 2, 3$ respectively. If you now look at the union of the three data sets you have created one new data set merged out of the three individual ones.

Next I would apply the moving average again (with a window length of now 6000) and compare the numbers and locations of peaks and throughs with the results on the three individual data sets and see whether this gives anything useful.

Answer 4

I cannot find any significant frequencies in any part of the data.

• Are your data stationary? - make tests (e.g. Dicky-Fuller) - if negative result transform your TSdata to statonary (e.g. taking first differencies)
• decompose for trend, cycles, noise - e.g. with python's statsmodels.tsa - if cycles exists, if residuals have some patterns
• in addition to seasonality (frequency) main feature of time-series is Volatility (for heteroscedastic data) - try use GARCH
• can estimate changing Volatility regime by Markov switching dynamic regression models or python's hmmlearn library -- for estimation of the probability of the phase
• if the final goal is forecast - you can benefit from quantile regression to estimate prediction uncertainty

My goal is to merge the three time series into a single one, so that I can more easily find the "upward" and "downward" phases. (Like a smooth sinus-curve or binary data-set 1=up, 0=down.)

it is difficalt to understand whether your 3 timeseries are 3 different features & if it is necessary to merge them in order to advise whether you need any sin/cos transformations - especially if you mentioned that you are not having any frequency. You can just investigate for the correlations of these features, as whuber♦ suggested in comments - e.g. with partial dependance for further creation of Multiple Linear Regression model (with possible interactions if exists) - can fit with any ML-tools (e.g. GradienBoosting)

if the goal is to identify "upward" and "downward" preassure, than having already differenced data can estimate Extremes (e.g. with python's pyextreme lib) relative any threshold - to mark them on unprocessed TS further - perhaps will see phase-changing points.

all possibilities are available for multi-dimentional data usually in available stat. soft & packages - if your data are such - just merge them to 3 columns and use such DataFrame in any stat soft/packages -- soft/package will do itself appropriate linear or non-linear transformations to smooth your data for signal detection in noice -- often according to some Signal/Noise Ratio estimation algorithms - either Kalman-Filtering or PCA or SpectrumEmbedding (last for AR also used)

Interactions among detrended (differenced) timeseries can also be estimated with ANOVA