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

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!

[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.

• What can you tell us about the statistical properties of these measurement systems? What we most need to know are their (a) biases and (b) second moments. The latter refers to possible correlations as well as different variances.
– whuber
Commented Dec 29, 2021 at 21:20
• @whuber I am unsure how to compute the bias, but have added mean, standard deviation and correlations. Commented Dec 29, 2021 at 22:03
• @litmus That's fine! You just have to use a lop-sided or even one-sided set of points of which to take the median when you are close to the edges of your overall time window. Commented Jan 6, 2022 at 14:31
• A lowpass filter is certainly going to do something similar. However, in my experience, there are several drawbacks: 1. You must know a fair amount about lowpass filters in order to get good performance. 2. The result of a lowpass filter does not contain actually measured values anymore. 3. In my experience, the signal-to-noise ratio gain from a lowpass filter isn't even close to a median filter. So what are the disadvantages of the median filter? The primary one is that the resulting data will not be as smooth as what you typically get from a lowpass filter. Commented Jan 6, 2022 at 19:09
• About how many up-periods and down-periods do you expect? As extremes, there might be three periods (if the event is a sine wave whose period is the length of the series) or hundreds (if the event is as spiky as the unsmoothed measurements), and those would suggest different models and techniques.
– user225256
Commented Jan 7, 2022 at 13:57

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)

• I have tried to use KF but it does not seem to work on the data I have provided (found in question, edit 2). I think it becomes much easier with artificial data, such as the type you created. Unfortunately it is my experiance so far that it is different from actual measured data. Would you like to try yourself with the data I provided in my question? Commented Jan 8, 2022 at 13:35

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

• Could you please elaborate? I looked up PCA on wikipedia but not sure what you mean. Also, I have attached links to the data if it is easy for you to try/show? Commented Jan 7, 2022 at 17:05
• take a look at example in a. different domain: clarusft.com/… - interest rates modeling in finance Commented Jan 7, 2022 at 17:08

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.

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