Consider a wide-sense stationary stochastic process $X(t)$ with constant mean $\mu$ and standard deviation $\sigma$. The (auto)correlation function $\rho(\tau)$ of $X(t)$ is defined as:

$$\rho(\tau) = \frac{\mathbb{E}[(X(t)-\mu)(X(t+\tau)-\mu)]}{\sigma^2}$$

My question relates to the estimation of the correlation function from a discretised realisation $\boldsymbol{X} = \{X_1,X_2,\dots,X_n\}$ of $X(t)$. A simple estimator $\hat{\rho}(\tau)$ is to compute Pearson's correlation coefficient for all pairs of lags separated by $k$:

\begin{equation} \hat{\rho}(k)=\frac{1}{(n-k) \sigma^2} \sum_{t=1}^{n-k} (X_t-\mu)(X_{t+k}-\mu) \end{equation}

Now suppose that instead of a single realisation of $X(t)$ being available there are $m$ realisations of $X(t)$ available. The estimator is easily extended to this case by considering all discretised points seperated by lag $k$ for all realisations:

\begin{equation} \hat{\rho}(k)=\frac{1}{j(n-k) \sigma^2} \sum_{j=1}^m \sum_{t=1}^{n-k} (X_{t}^j-\mu)(X_{t+k}^j-\mu) \end{equation}

However these estimators are computationally inefficient in comparison to an approach depending on the signal processing definition of autocorrelation as the convolution of the signal with itself:

$$\rho(\tau) = X(t) * \bar{X}(t) = \int_{-\infty}^\infty X(t)X(t-\tau)\, {\rm d}t$$

It is then possible to use the convolution theorem to write:

$$X(t) * \overline{{X}(t)} = \mathcal{F}^{-1}\big\{\mathcal{F}\{X(t)\}\cdot\overline{\mathcal{F}\{X(t)\}}\big\}$$

Where $\mathcal{F}\{X(t)\}$ is the Fourier transform of $X(t)$ and $\overline{\mathcal{F}\{X(t)\}}$ is its complex conjugate. This leads to a computationally efficient estimator for $\rho(\tau)$ that makes use of optimised fast Fourier transform algorithms:

\begin{equation} \hat{\rho}(k) = \mathcal{F}^{-1}\big\{\mathcal{F}\{\boldsymbol{X}\}\cdot\overline{\mathcal{F}\{\boldsymbol{X}\}}\big\} \end{equation}

Where $\boldsymbol{X}$ is a discretised single realisation of $X(t)$ (as above). My question is whether it possible to extend this approach to the case where there are $m$ realisations of $X(t)$? Say that $\boldsymbol{X}$ is an $n$ by $m$ matrix instead of a vector.

  • $\begingroup$ If $ \mathbb{E}[X_t]$ and $R_k = \mathbb{E}[(X_t-\mu)( X_{t+k}-\mu)]$ is constant then $R_k =\lim_{A \to \infty} \frac{1}{2A} \mathbb{E}[\int_{-A}^A (X_t-\mu)( X_{t+k}-\mu)dt] = \lim_{N \to \infty} \frac{1}{2N} \mathbb{E}[\sum_{n=-N}^N (X_n-\mu)(X_{n+k}- \mu)]=\frac{1}{2N} \sum_{n=-N}^N \mathbb{E}[(X_n-\mu)(X_{n+k}- \mu)]$. $\endgroup$
    – reuns
    Aug 29 '17 at 16:12
  • $\begingroup$ For this last one, a good estimator is $\widehat{R}_k=\frac{1}{2N m}\sum_{j=1}^m \sum_{n=-N}^N (x_n^{(j)}-\mu)(x_{n+k}^{(j)}- \mu)$ where $x^{(j)},j=1 \ldots m$ are realizations of your process. You can compute $\widehat{R}_k$ (for $k = 1 \ldots K$) efficiently from the FFT of $\frac{1}{m} \sum_{j=1}^m x^{(j)}_n-\mu, n = -N \ldots N-K$ $\endgroup$
    – reuns
    Aug 29 '17 at 16:13
  • $\begingroup$ Did you get what I wrote ? $\endgroup$
    – reuns
    Aug 31 '17 at 22:33

Since the DFT is a Linear Operator you can apply the averaging in the Fourier Domain and have the same result as in the time domain.

In practice, just pay attention that DFT which is what used applies Circular Convolution.
You either will be OK with that or pad your signals accordingly.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.