# Confusion about ensembles and averages in autocorrelation matrices

I just noticed that until now I often don't cared about the scaling of the autocorrelation Matrix in Matlab. I then checked with the book Statistical Digital Signal Processing from M. H. Hayes. There it is defined like this:

First of all, why does the author note it like $$x(0)x^*(0)$$? Written like this wouldm't this mean first element of the vector x times first element of the vector x? And the expected value is then what $$\frac{x(0)x^*(0)}{2}$$?

$$r_x(\tau)$$ can be calculated as $$E[X(t)X(t+\tau)]$$ This is more what I would expect. A pointwhise vector multiplication of a $$\tau$$ shifted vector.

Here comes my actual question. For e.g. $$x = [1 ,2, 3]^T$$, $$r_x(0)$$ should become $$\frac{1^2+2^2+3^2}{3}=\frac{14}{3}$$ in my oppinion. Matlab shows for xcorr(x)

3.0000    8.0000   14.0000    8.0000    3.0000


$$r_x(0)$$ should be the value in the middle. But why don't they take the expected value? On this page they suggest that the autocorrelation matrix can be calculated with

z=autocorr(x)


and

Rxx=toeplitz(z)


Instead of autocorr, xcorr(x) should work too stated other posts. As expected this results in

Rx =

3.0000    8.0000   14.0000    8.0000    3.0000
8.0000    3.0000    8.0000   14.0000    8.0000
14.0000    8.0000    3.0000    8.0000   14.0000
8.0000   14.0000    8.0000    3.0000    8.0000
3.0000    8.0000   14.0000    8.0000    3.0000


Again, no expected values. And this solution is also wrong because it is shifted. $$r_x(0)$$ is not in the upper left corner. I could shift it by hand, but is that how autocorrelation matrixes are calculated in Matlab? I checked with some other sources in the internet. Often the algorithm to calculate an entry of an autocorrelation contains just the sum.

• one averages over an ensemble of samples, not lags
– user28715
Feb 3, 2019 at 22:21

First, the output of the xcorr() function returns lag-0 of the auto-correlation sequence (ACS) estimate at its middle sample, as you recognize. The function argument scaleopt provides normalization of the output: the default value none causes your confusion. If you select biased, instead, then Eq-7 (below) and xcorr() outputs will be identical, like this:

rx = xcorr([1,2,3],'biased'); % identical to (1/3)*xcorr([1,2,3])


Second, do not use the full output of xcorr as the input to the toeplitz function, since the xcorr() output will include $$r_x[k]$$ estimates for all lags $$r_x[-2],r_x[-1],r_x[0],r_x[1],r_x[2]$$. But you only need $$r_x[0],r_x[1],r_x[2]$$. Instead, call it like:

L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) );  % now you're ok!


Finally, distinguish (theoretical) expectations from their (practical) estimations. When algorithms work on data they don't calculate expectations, they estimate. Hence, all expectations are replaced with their estimates in practical applications.

NOTE: Below is some optional discussion.

• Theoretical definition of the ACS $$r_{XX}[k]$$ of an RP $$X[n,s]$$.
• Practical estimation of ACS, from a sample sequence $$x[n]$$ (one realization of the RP) via ergodicity.
• Utilization of a random-vector

Consider the interpretation of an RP $$X[n,s]$$ as an ordered collection (set) of indexed RVs $$X_n(s)$$ for each time-index $$n$$. Each $$X_n(s)$$, in principle, has its own probability distribution/density function $$F_{X_n}(x)$$.

Unfortunately, this yields enormously large number and type of possible joint density functions needed to characterize the RP, therefore, things are simplified by assuming that all those RVs have the same (identical) PDF $$F_{X_n}(x)$$ for all $$n$$, and that all are independent of each other, and further that first & second moments (mean,variance and aurocorrelations) are independent of partciular time-instants. This simplification leaves an i.i.d. (independent, identically distributed) WSS RP $$X[n,s)$$ at hand.

One more simplification (notationally) is to replace the notation $$X_n(s)$$ with $$x[n]$$ to represent the RP. It may cause some confusions, and I will continue with $$X_n(s)$$ below.

The (theoretical) ACS of an i.i.d. WSS RP is defined as:

$$r_x[k] = r_{x x}[k] = r_{X_n X_{n+k}}[n,n+k]= E\{X_n X^*_{n+k}\} = E\{x[n]x^*[n+k] \} \tag{1}$$

Calculation of $$r_x[k]$$ requires evaluation of expectations: $$E\{X_n X_m\} = \int_{-\infty}^{\infty} x_1 x_2 f_{X_n,X_m}(x_1,x_2) dx_1 dx_2$$. This is ok with questions on paper, but for practical applications it's estimated from data, based on ergodicty assumption: theoretical ensemble-averages can be estimated from sample function time-averages. Then you have

$$\hat{r}_X[k] = \hat{r}_x[k] = \lim_{N \to \infty} \frac{1}{2N+1} \sum_{n=-N}^{N} x[n]x^*[n+k] \tag{2}$$

Note that $$x[n]$$ on RHS is a particular sample function (a realization) of the RP $$X[n,s)$$. Note also that the summation limits are determined by data window boundaries.

Now consider the random vector $$\bar{X}_N(n) = [X_n , X_{n-1} , ..., X_{n-N+1}]^T$$ which includes $$N$$ consecutive RVs of the RP $$X[n,s)$$, beginning at index $$n$$.

From this vector, generate a Random Matrix (RM) via the Hermitian outer product $$\bar{X}_N \cdot \bar{X}_N^H$$: $$RM = \begin{bmatrix} X_n X^*_n & X_n X^*_{n-1} & ... & X_n X^*_{n-N+1} \\ X_{n-1} X^*_n & X_{n-1} X^*_{n-1} & ... & X_{n-1} X^*_{n-N+1} \\ ... & ... & ... & ... \\ X_{n-N+1} X^*_n & X_{n-N+1} X^*_{n-1} & ... & X_{n-N+1} X^*_{n-N+1} \\ \end{bmatrix} \tag{3}$$

Taking the expectation of this RM, yields the theoretical auto-correlation matrix (ACM) $$R_x$$ of the WSS RP $$X[n,s)$$ :

$$R_x = R_{XX} = E\{RM\} = E\{ \bar{X_N} \cdot \bar{X_N}^H \} = E\{\begin{bmatrix} X_n X^*_n & X_n X^*_{n-1} & ... & X_n X^*_{n-N+1} \\ X_{n-1} X^*_n & X_{n-1} X^*_{n-1} & ... & X_{n-1} X^*_{n-N+1} \\ ... & ... & ... & ... \\ X_{n-N+1} X^*_n & X_{n-N+1} X^*_{n-1} & ... & X_{n-N+1} X^*_{n-N+1} \\ \end{bmatrix} \} \tag{4}$$

and this results in: $$R_x = E\{ \bar{X}_N \cdot \bar{X}_N^H \} = \begin{bmatrix} r_x[0] & r_x[-1] & ... & r_x[-N+1] \\ r_x[1] & r_x[0] & ... & r_x[-N+2] \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x[0] \\ \end{bmatrix} \tag{5}$$

where $$E\{X_n {X^*}_m\}= r_x[n-m]$$, and becuse of Hermitian symmetry of Auto-Correlations, it becomes: $$R_x = \begin{bmatrix} r_x[0] & r_x^*[1] & ... & r_x^*[N-1] \\ r_x[1] & r_x[0] & ... & r_x^*[N-2] \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x[0] \\ \end{bmatrix} \tag{6}$$

Again, this theoretical ACM $$R_x$$ is estimated from sample observation data.

Estimation of the ACM $$R_x$$ can be done in a number of related ways, but basically, it reduces into an estimation of the ACS $$r_x[k]$$ for the lags $$k=0,1,...,N-1$$. Then from the Toeplitz and Hermitian symmetric property of $$R_x$$, you can construct it.

One way to practically estimate ACS $$r_x[k]$$ from the data $$x[n]$$ of length $$L$$ is :

$$\hat{r}[k] = \frac{1}{L} \left( x[k] \star x^*[-k] \right) \tag{7}$$

This yields an estimate of ACS for the lags $$-L < k < L$$, (you only need $$N$$ of those for $$k=0,1,...,N-1$$)

Alternatively, ACM can be estimated similarly to ACS as:

$$\hat{R}_X = \frac{1}{L} \sum_{k=0}^{L-1} \bar{X}_N(k) \cdot \bar{X}_N(k)^H \tag{8}$$ where the data vector $$X_N(k)$$ includes $$N$$ consecutive samples of $$x[n]$$ beginning from $$x[k]$$, and $$L$$ is the length of the observation sequence $$x[n]$$.

• That was a very good answer. However, I'm still confused about Matlabs output. If I use your second equation for $\hat{r}_X[k]$ for lag 0 (the upper left most entry) then I would not get $1+2^2+3^2=14$. From intuition and my daily use of convolution I know thats true for lag 0. But I can't connect it with the newest theory I'm learning about autocorrelation matrices. All those formulas do not result in what I used to know. I still make mistakes regarding ensembles and means. But again, with all I've read above, this doesn't fit. There should at least be a 1/N or (1/K) factor for the estimation Feb 4, 2019 at 23:39
• Regarding the scaling, all estimation equations do have a scale, as you can see too. You are somewhat right that statistical signal processing can provide puzzling equations at first but they all turn to be same later. May be you should try to focus on the very step that causes the first confusion, rather than bundling them together. Feb 5, 2019 at 1:06
• Yes, thats perfectly what I was confused about. I'm sorry I didn't explained it well the first time. It's actually pretty hard to say what you don't understand yet. On page 393 spectrum estimation Hayes states the same formula to estimate $\hat{r}_x(k)$ as you. However, in our exam cheat sheet ACS is calculated without 1/N. This is either a mistake or one of those damn conventions. I hope this gets better with time. At the moment I'm freaking out every time a formula doesn't look exactly like I'm used to. Feb 5, 2019 at 14:56
• Well... I hate it. Guess I just got it. In the All-Pole Modeling were the autocorrelation matrix is used ($R_x*a_p=-r_x$) you can just factor out 1/N on both sides. That's probably why they don't bother calculating the autocorrelation correctly. I'm freaking out. Why can't people just write it down like it is :) Feb 5, 2019 at 15:00
• @Mr.Sh4nnon yes! that's the reason. For computational efficiency purposes, if both sides of the linear equation are weighted by $1/N$ (which is the case for $R_x a_p = r_x$ all-pole ,Yule-Walker eqs) then yo ucan ignore them. However, be very cautious that ACM $R_x$ or ACS $r_x$ when computed alone, do need the weight in for their ergodic estimates from data. Feb 5, 2019 at 17:19

Thanks for the answer @Fat32! Please allow me to add my notes as an answer hopping to get any comments for corrections and further understanding.

Give, $$\tilde{X}(n) = [X(n),X(n-1),..,X(n-N+1)]^T \in \mathbb{R}^N$$, $$n=N,\dots, L$$, where $$N$$ is the lag size and $$L$$ the signal size, I was trying to understand how the sample covariance matrix

$$\tilde{R}_{X} = \frac{1}{L} \sum_{t=N}^L \tilde{Χ}(t)\tilde{Χ}^*(t) = \frac{1}{L}\begin{bmatrix} \sum_{t=N}^L X^2(t) & \sum_{t=N}^L X(t-1) X(t) & \cdots &\sum_{t=N}^L X(t-N+1)X(t)\\ \sum_{t=N}^L X(t)X(t-1) & \sum_{t=N}^L X^2(t-1) & \cdots &\sum_{t=N}^L X(t-N+1)X(t-1) \\ \vdots & \vdots & \cdots & \vdots \\ \sum_{t=N}^L X(t)X(t-N+1) & \sum_{t=N}^L X(t-1)X(t-N+1) & \cdots &\sum_{t=N}^L X^2(t-N+1) \end{bmatrix}\tag{A}$$

can be constructed using the auto-correlation matrix. By definition the auto-correlation matrix is given by

$$R_X = E\{ X \cdot X^H \}\tag{B}$$

Next, we show $$\tilde{R}_X \stackrel{\Delta}{=} \frac{1}{L}R_X \tag{C}$$

We proceed as, $$\tilde{R}_X = \frac{1}{L} R_X =\frac{1}{L}E\{ X \cdot X^H \} = \\ \frac{1}{L} \begin{bmatrix} E\{X_1 X^*_1\} & E\{X_1 X^*_2\} & ... & E\{X_1 X^*_N\} \\ E\{X_1 X^*_2\} & E\{X_2 X^*_2\} & ... & E\{X_2 X^*_N\} \\ ... & ... & ... & ... \\ E\{X_1 X^*_N\} & E\{X_1 X^*_N\} & ... & E\{X_N X^*_N\} \\ \end{bmatrix} = \\ \frac{1}{L}\begin{bmatrix} r_x[0] & r_x[-1] & ... & r_x[-N+1] \\ r_x[1] & r_x[0] & ... & r_x[-N+2] \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x[0] \\ \end{bmatrix} = \\\frac{1}{L} \begin{bmatrix} r_x[0] & r_x[1]^* & ... & r_x[N-1]^* \\ r_x[1] & r_x[0] & ... & r_x[N-2]^* \\ ... & ... & ... & ... \\ r_x[N-1] & r_x[N-2] & ... & r_x[0] \\ \end{bmatrix} \tag{D}$$ where

$$X = \begin{bmatrix} {X}_1 \\ {X}_2 \\ {X}_{N} \end{bmatrix} = \begin{bmatrix} \tilde{X}(1) \\ \tilde{X}(2) \\ \tilde{X}(L-N+1) \end{bmatrix} = \\ \begin{bmatrix} X(1) & X(2) & \dots & X(N) \\ X(2) & X(3) & \dots & X(N+1) \\ X(L-N+1) & X(L-N+2) & \dots & X(L) \end{bmatrix}, \tag{E}$$ $$E\{X_n X^*_m\} = r_x[n-m]$$ and

$$r_x [k] = \sum_{t=N}^L X[t] X^*[t+k] \tag{F}$$

the auto-correlation coefficient [Section 3.3.4, 1]. Using, (F) in (D) we get (A) and the proof of (C) is complete.

Here is a Python code

import numpy as np
from scipy import signal
from scipy.linalg import toeplitz

if __name__ == '__main__':
rng = np.random.default_rng()

sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 12)
sig_noise = sig + rng.standard_normal(len(sig))

L = len(sig_noise)
mid_point = int(np.floor(L / 2))  # because of symmetry in auto-correlation, we keep the half of coefficients
corr = signal.correlate(sig_noise, sig_noise, mode='same')[mid_point:] / L
tildeRx = toeplitz(corr, np.conjugate(corr))


Here is also a Python brute force computation of

$$\tilde{R}_{X} = \frac{1}{L} \sum_{t=N}^L \tilde{Χ}(t)\tilde{Χ}^*(t)\tag{E}$$

# Brute force computation
_tildeRx = np.zeros((mid_point, mid_point))
t_range = np.arange(mid_point, len(sig_noise))
for t in t_range:
y_tilde = np.matrix(sig_noise[t - mid_point: t][::-1]).T
_tildeRx += y_tilde * y_tilde.H
_tildeRx /= L


[1] Monson H. Hayes, Statistical Digital Signal Processing and Modeling 1st edition, 1996

Note: I am not getting the same output with the two methods I am still working on it.

• I believe this is not meant to be an answer. It seems like you are trying to initiate a conversation here instead of actually providing a (possible) solution to the question. This is not how things work here though. This is a Q&A site and not a forum. You seem to already have enough reputation to know that, so if I am mistaken please let me know. Otherwise, I urge you to delete the post and when you finally come up with an answer then post it here. If you would like to chat with people about an answer or a problem (which is quite normal) please use the chat. Aug 11, 2023 at 15:59
• Hi! I made an edit. May be it's clearer now. Aug 12, 2023 at 0:14
• @Fat32 thanks for the response! In your post, in $(8)$, did you mean $\tilde{X}_N$ instead of $X_N$ ? In addition, in the same equation, the range of sum shouldn't be $N$ to $L$ instead of $0$ to $L-1$? Moreover, in the same equation, I am not sure about the scale fraction. Should it be $1/N$ or $1/L$. Can you take a closer look? Aug 12, 2023 at 10:57
• @Thoth No there's no tilde in Eq.8. But my definition of rx[k] and X are different than standard convention. Sepcifically, my random vector X uses causal samples, whereas those in some books use future samples, etc. Regarding the scale: you should consider data window boundary (edge) conditions specifically while performing estimations. There are different conventions which assume "0" or "nothing" out of the boundaries. Eq.8 is generic, you should implement it with specific assumption of yours. I will later check it again. Aug 12, 2023 at 11:17
• there's no tilde, but a bar . I ignored it assuming it was apparent. Again note that I used causal random vector, and most books use noncausal (offline) form. Yet both yield similar results. Aug 12, 2023 at 11:35