Auto-Correlation Matrix in Matlab & Python: Toeplitz Form

by

An auto-correlation matrix collects the autocorrelation sequence of a signal into a Toeplitz matrix.

Learning path: from convolution to correlation matrices

Work through these in order, then come back here:

  1. Polynomials, convolution and Toeplitz matrices: the Toeplitz structure that this post relies on.
  2. Methods to compute linear convolution: the convolution used to get the autocorrelation sequence.
  3. AutoCorrelation (Correlogram) and persistence: what the autocorrelation sequence measures.
  4. Auto-correlation matrix (this post): building the matrix from the sequence.
  5. Yule Walker Estimation: an application that solves a system built from the same Toeplitz matrix.
  6. Generate correlated Gaussian sequence (colored noise): using correlated data in simulations.

Autocorrelation sequence and convention

For a stationary complex sequence $x[n]$, the autocorrelation at lag $k$ is

$$r(k) = E{x[n+k]\,\overline{x[n]}}$$

where the bar denotes complex conjugation. For real signals the conjugate has no effect. In practice the sequence is estimated from $N$ samples with the biased estimator, which divides by $N$ at every lag. Here $\overline{r(k)} = r(-k)$.

The cross-correlation of two sequences uses the same convention, with the conjugate on the first sequence. This is the Wikipedia convention. The original version of this post put the conjugate on the second sequence, which gives the complex conjugate of each lag. The two conventions differ only by $k \to -k$, but mixing them changes the sign of the imaginary parts of the matrix.

Autocorrelation matrix

Take $N$ consecutive samples as a vector $\mathbf{x}_n = [x[n], x[n+1], \dots, x[n+N-1]]^T$. The autocorrelation matrix is

$$\mathbf{R} = E{\mathbf{x}_n \mathbf{x}_n^H}$$

Its entries depend only on the lag between the row and column indices, so $\mathbf{R}$ is a Toeplitz matrix:

$$R_{ij} = r(i-j) \quad \text{for } i \ge j, \qquad R_{ij} = \overline{r(j-i)} \quad \text{for } i < j$$

The matrix is also Hermitian, $\mathbf{R} = \mathbf{R}^H$. These two properties let you build $\mathbf{R}$ from $N$ numbers, $r(0), \dots, r(N-1)$. For snapshot data, $\hat{\mathbf{R}} = \frac{1}{M}\sum_{m=1}^{M} \mathbf{x}_m \mathbf{x}_m^H$ gives the sample version.

Computing the autocorrelation sequence

There are several ways to compute $r(k)$. All four methods below give the same answer for $x = [1, 2, 3, 4]$, which is $[4, 11, 20, 30, 20, 11, 4]$. The center element is lag 0.

Method 1: xcorr (Matlab) and numpy.correlate (Python)

Matlab:

x = [1 2 3 4];
acf = xcorr(x)          % 4 11 20 30 20 11 4

Python:

import numpy as np
x = np.asarray([1, 2, 3, 4])
np.correlate(x, x, mode='full')   # array([ 4, 11, 20, 30, 20, 11,  4])

Method 2: convolution with the flipped conjugate

The autocorrelation is the convolution of $x$ with the flipped conjugate of $x$. The conjugate is not needed for real input.

Matlab:

x = [1 2 3 4];
acf = conv(x, fliplr(conj(x)))   % 4 11 20 30 20 11 4

Python:

import numpy as np
x = np.asarray([1, 2, 3, 4])
np.convolve(x, np.conj(x)[::-1])   # array([ 4, 11, 20, 30, 20, 11,  4])

Method 3: Toeplitz matrix times the flipped input

Matlab:

x = [1 2 3 4];
acf = toeplitz([x zeros(1, length(x)-1)], ...
      [x(1) zeros(1, length(x)-1)]) * fliplr(x).'   % 4 11 20 30 20 11 4

Python:

import numpy as np
from scipy.linalg import toeplitz
x = np.asarray([1, 2, 3, 4])
col = np.pad(x, (0, len(x) - 1))
row = np.pad([x[0]], (0, len(x) - 1))
toeplitz(col, row) @ x[::-1]   # array([ 4, 11, 20, 30, 20, 11,  4])

Method 4: FFT and IFFT

Zero-pad both sequences to length $2N-1$ so the circular convolution equals the linear one.

Matlab:

x = [1 2 3 4];
L = 2*length(x) - 1;
acf = ifft(fft(x, L) .* fft(fliplr(conj(x)), L))   % 4 11 20 30 20 11 4

Python:

import numpy as np
x = np.asarray([1, 2, 3, 4])
L = 2 * len(x) - 1
np.fft.ifft(np.fft.fft(x, L) * np.fft.fft(np.conj(x)[::-1], L))   # 4, 11, 20, 30, 20, 11, 4 (complex, zero imaginary part)

Building the matrix

Once $r(0), \dots, r(N-1)$ are known, the matrix is one call to toeplitz. Matlab’s toeplitz(c, r) takes the first column $c$ and the first row $r$. Use the lag values for $c$ and their conjugates for the first row.

Matlab, for a complex input:

x = [1+1j, 2+1j, 3-1j];
acf = conv(x, fliplr(conj(x)));   % Method 2
Rxx = acf(3:end);                 % Rxx(1) is lag 0, the center element
Rx = toeplitz(Rxx, [Rxx(1), conj(Rxx(2:end))])

Python, for the same input:

import numpy as np
from scipy.linalg import toeplitz
x = np.asarray([1+1j, 2+1j, 3-1j])
acf = np.convolve(x, np.conj(x)[::-1])   # Method 2
Rxx = acf[2:]                            # Rxx[0] is lag 0
Rx = toeplitz(Rxx, np.hstack((Rxx[0], np.conj(Rxx[1:]))))

Both give the same $3 \times 3$ matrix, which is Hermitian:

17 + 0i    8 + 6i    2 + 4i
 8 - 6i   17 + 0i    8 + 6i
 2 - 4i    8 - 6i   17 + 0i

The first row is $r(0)$, then the conjugates $\overline{r(1)}$ and $\overline{r(2)}$. The first column is $r(0)$, $r(1)$, $r(2)$. Each diagonal is constant, which is the Toeplitz property.

Simulation: AR(1) data

To see the structure on real data, generate an AR(1) process $x[n] = 0.8\,x[n-1] + w[n]$ and estimate its autocorrelation. The theoretical autocorrelation is $r(k)/r(0) = 0.8^{k}$. The matrix below is the $8 \times 8$ estimate built from the same sequence.

Estimated autocorrelation of AR(1) data against 0.8 to the power k, and the 8 by 8 Toeplitz autocorrelation matrix
Figure 1: Left, the estimated autocorrelation of 4000 AR(1) samples matches $0.8^{k}$. Right, the $8 \times 8$ matrix is constant along each diagonal.

The matrix is what the Wiener filter need. The Wiener filter solves $\mathbf{R}\mathbf{w} = \mathbf{p}$ for the filter weights.

Summary

The auto-correlation matrix is a Toeplitz, Hermitian matrix built from the autocorrelation sequence. Compute the sequence with any of the four methods, keep the conjugate convention consistent, and build $\mathbf{R}$ with toeplitz.

References

[1] S. S. Reddi, “Eigenvector properties of Toeplitz matrices and their application to spectral analysis of time series”, Signal Processing, vol. 7, North-Holland, 1984, pp. 46-56. https://doi.org/10.1016/0165-1684(84)90023-9

[2] Robert M. Gray, “Toeplitz and Circulant Matrices: A Review”, Stanford University. https://ee.stanford.edu/~gray/toeplitz.pdf

Further reading: Autocorrelation (Wikipedia), Toeplitz matrix (Wikipedia), Hermitian matrix (Wikipedia), Wiener filter (Wikipedia), Beamforming (Wikipedia), numpy.correlate, and scipy.linalg.toeplitz.

Build the simulations yourself. The GaussianWaves eBooks include complete, runnable simulation code, beyond what fits in a single blog post (check table of contents)

  • Digital Modulations using Python: modem and BER simulations built from scratch with NumPy, SciPy, and Matplotlib. Paperback on Amazon, or PDF with instant download.
  • Wireless Communication Systems in Matlab (Second Edition): end-to-end wireless simulations, from modulation and coding to OFDM and fading channels, with Matlab scripts.

Browse all GaussianWaves eBooks and products

Leave a Comment

3 thoughts on “Auto-Correlation Matrix in Matlab & Python: Toeplitz Form”