This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. Please direct questions and suggestions to Sascha.Spors@uni-rostock.de.
Correlation is a statistical measure for the dependencies between random processes or between the samples of one random process. The auto-correlation function (ACF) characterizes the temporal dependencies within one random signal $x[k]$. It is an important measure for the analysis of signals in communications engineering, source coding and system identification.
For a continuous-amplitude real-valued random signal $x[k]$ the ACF is defined by the second-order ensemble avarage of the signal at two different time-instants $k_1$ and $k_2$
\begin{equation} \varphi_{xx}[k_1, k_2] = E\{ x[k_1] \cdot x[k_2] \} \end{equation}Under the assumption of wide-sense stationarity (WSS) the ACF does only depend on the difference $\kappa = k_1 - k_2$ between the considered sample indexes
\begin{equation} \varphi_{xx}[\kappa] = E\{x[k] \cdot x[k-\kappa] \} \end{equation}where $\kappa$ is commonly chosen as sample index instead of $k$ in order to indicate that it denotes a shift/lag. The ACF quantifies the similarity of a signal with a shifted version of itself. It has high values for high similarity and low values for low similarity.
If the process is additionally wide-sense ergodic, the ACF can be computed by averaging along one sample function
\begin{equation} \varphi_{xx}[\kappa] = \lim_{K \to \infty} \frac{1}{2K + 1} \sum_{k=-K}^{K} x[k] \cdot x[k-\kappa] \end{equation}Note that the normalization on the left side of the sum is discarded in some definitions of the ACF. Above summation resembles strongly the definition of the discrete convolution. For a random signal $x_N[k] = \text{rect}_N[k] \cdot x[k]$ of finite length $N$ and by exploiting the properties of a wide-sense ergodic random process one yields
\begin{equation} \varphi_{xx}[\kappa] = \frac{1}{N} \sum_{k=0}^{N-1} x_N[k] \cdot x_N[k-\kappa] = \frac{1}{N} \; x_N[k] * x_N[-k] \end{equation}where the ACF $\varphi_{xx}[\kappa] = 0$ for $|\kappa| > N-1$. Hence, the ACF can be computed by (fast) convolution of the random signal with a time reversed version of itself.
Note in numerical implementations (e.g. MATLAB, Python), the computed ACF is stored in a vector of length $2 N - 1$. The positive indexes $0, 1, \dots, 2 N - 1$ of this vector cannot be directly interpreted as $\kappa$. The indexes of the vector have to be shifted by $N-1$ for a proper interpretation.
The following properties of the ACF can be deduced from its definition
The ACF $\varphi_{xx}[\kappa]$ has a maximum for $\kappa = 0$. It is given as
$$ \varphi_{xx}[0] = E\{x^2[k]\} = \sigma_x^2 + \mu_x^2 $$
This is due to the fact that the signal is equal to itself for $\kappa = 0$. Please note that for periodic random signals more than one maximum will be present.
The ACF is a function with even symmetry
$$ \varphi_{xx}[\kappa] = \varphi_{xx}[-\kappa] $$
For typical random signals, the ACF approaches the limiting value
$$ \lim_{|\kappa| \to \infty} \varphi_{xx}[\kappa] = \mu_x^2 $$
The similarity of a typical random signal is often low for large lags $\kappa$.
The ACF of a periodic signal $x[k] = x[k + P]$ is also periodic
$$ \varphi_{xx}[\kappa] = \varphi_{xx}[\kappa + P] $$
with the period $P \in \mathbb{N} \setminus 0$
A random signal $x[k]$ is said to be uncorrelated if
\begin{equation} \varphi_{xx}[\kappa] = \mu_x^2 + \sigma_x^2 \cdot \delta[\kappa] \end{equation}and as correlated if this condition is not met. The samples of a signal which is uncorrelated show no dependencies between each other in a statistical sense.
A periodic signal $x[k]$ can be expressed as
\begin{equation} x[k] = x_0[k] * \sum_{n = - \infty}^{\infty} \delta[k - nP] \end{equation}where $x_0[k]$ denotes the signal of one period. The DTFT $X(e^{j \Omega})$ of a periodic signal consists of a series of equidistant Dirac impulses (line spectrum). Expressing the ACF as convolution, transforming this into the frequency domain by the DTFT and having the multiplication property of the Dirac delta in mind it can be concluded that the ACF is periodic with period $P$. The ACF is then given by convolution over one period only and periodic continuation
\begin{equation} \varphi_{xx}[\kappa] = \frac{1}{N} \left( x_0[\kappa] * x_0[-\kappa] \right) * \sum_{n = - \infty}^{\infty} \delta[k - nP] \end{equation}The ACF of a periodic signal is computed and plotted in the following.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
N = 8192 # total length of signal
P = 128 # period
K = 275 # upper/lower limit for lag in ACF
# generate periodic random signal
np.random.seed(1)
x0 = np.random.normal(size=P)
x = np.tile(x0, N//P)
# compute and truncate ACF
acf = 1/len(x) * np.correlate(x, x, mode='full')
acf = acf[(len(x)-1)-(K-1):(len(x)-1)+K]
kappa = np.arange(-(K-1), K)
# plot signal and its ACF
plt.figure(figsize = (10, 4))
plt.stem(x[:512])
plt.xlabel(r'$k$')
plt.ylabel(r'$x[k]$')
plt.grid()
plt.figure(figsize = (10, 4))
plt.stem(kappa, acf)
plt.xlabel(r'$\kappa$')
plt.ylabel(r'$\hat{\varphi}_{xx}[\kappa]$')
plt.axis([-K, K, 1.1*min(acf), 1.1*max(acf)]);
plt.grid()
Exercise
Solution: The period of the random signal can be estimated from the ACF, as both have the same period. This is most easy by inspection of the regular peaks of the ACF.
The following example estimates and plots the ACF of a short recorded speech signal.
from scipy.io import wavfile
K = 30 # upper/lower limit for lag in ACF
# read audio file
fs, x = wavfile.read('../data/speech_8k.wav')
x = np.asarray(x, dtype=float)/2**15
# compute and truncate ACF
acf = 1/len(x) * np.correlate(x, x, mode='full')
acf = acf[(len(x)-1)-(K-1):(len(x)-1)+K]
kappa = np.arange(-(K-1), K)
# plot ACF
plt.figure(figsize = (10, 8))
plt.stem(kappa, acf)
plt.xlabel(r'$\kappa$')
plt.ylabel(r'$\hat{\varphi}_{xx}[\kappa]$')
plt.axis([-K, K, 1.1*min(acf), 1.1*max(acf)]);
plt.grid()
Exercise
K
for $-K \leq \kappa < K$ and check if the last property is fulfilled for the speech signal.Solution: It can be observed that the plotted ACF has its maximum value at $\kappa = 0$ and that it shows an even symmetry. The third property can be confirmed by increasing K
in above example. The speech signal is correlated since $\varphi_{xx}[\kappa] \neq 0$ for $\kappa \neq 0$.
The auto-covariance function is defined as the ACF of a zero-mean random signal. For a random signal $x[k]$ with linear mean $\mu_x \neq 0$ it is given as
\begin{equation} \psi_{xx}[\kappa] = \varphi_{xx}[\kappa] - \mu_x^2 \end{equation}The cross-correlation function (CCF) is a measure of similarity that two random signals $x[k]$ and $y[k]$ show in a statistical sense.
For a continuous-amplitude real-valued random signal $x[k]$ the CCF is defined by the second-order ensemble avarage of the two signals $x[k]$ and $y[k]$ at two different time-instants $k_x$ and $k_y$
\begin{equation} \varphi_{xy}[k_x, k_y] = E\{ x[k_x] \cdot y[k_y] \} \end{equation}Under the assumption of wide-sense stationarity (WSS) the CCF does only depend on the difference $\kappa = k_x - k_y$ between the considered sample indexes
\begin{equation} \varphi_{xy}[\kappa] = E\{x[k] \cdot y[k - \kappa] \} = E\{x[k + \kappa] \cdot y[k] \} \end{equation}The cross-correlation function (CCF) is a measure of similarity that two random signals $x[k]$ and $y[k - \kappa]$ have with respect to the shift $\kappa \in \mathbb{Z}$. If $x[k]$ and $y[k]$ are wide-sense ergodic processes, the CCF can be computed by averaging along one sample function
\begin{equation} \varphi_{xy}[\kappa] = \lim_{K \to \infty} \frac{1}{2K + 1} \sum_{k=-K}^{K} x[k] \cdot y[k-\kappa] \end{equation}For random signals $x_N[k] = \text{rect}_N[k] \cdot x[k]$ and $y_M[k] = \text{rect}_M[k] \cdot y[k]$ of finite lengths $N$ and $M$ one yields
\begin{equation} \varphi_{xy}[\kappa] = \frac{1}{N} \sum_{k=0}^{N-1} x_N[k] \cdot y_M[k-\kappa] = \frac{1}{N} \; x_N[k] * y_M[-k] \end{equation}where the CCF $\varphi_{xy}[\kappa] = 0$ for $\kappa < -(M-1)$ and $\kappa > N-1$. The CCF can be computed by (fast) convolution of one random signal with a time reversed version of the other random signal. Note in numerical implementations (e.g. MATLAB, Python), the computed CCF is stored in a vector of length $N + M - 1$. The positive indexes $0, 1, \dots, N + M - 1$ of this vector cannot be directly interpreted as $\kappa$. The indexes of the vector have to be shifted by $M-1$ for a proper interpretation.
Above results hold also for the CCF $\varphi_{yx}[\kappa]$ when exchanging $x[k]$ with $y[k]$ and $N$ with $M$.
When exchanging the two random signals, the CCF exhibits the following symmetry
$$ \varphi_{xy}[\kappa] = \varphi_{yx}[-\kappa] $$
For typical random processes, the CCF approaches the limiting value
$$ \lim_{|\kappa| \to \infty} \varphi_{xy}[\kappa] = \mu_x \cdot \mu_y$$
Two random signals are said to be uncorrelated if
\begin{equation} \varphi_{xy}[\kappa] = \mu_x \cdot \mu_y \end{equation}and as correlated if this condition is not met. The samples of two signals which are uncorrelated to each other show no dependencies between each other in a statistical sense. $\varphi_{xy}[\kappa] = 0$ for uncorrelated signals if one of the two random processes is zero-mean.
The CCF can be used to estimate the time-of-arrival (TOA) of a signal at transmitted from a sender to a receiver. The TOA is used for instance for radiolocation or sound source localization. For ease of illustration it is assumed that a random signal $x[k]$ is only delayed by transmission. The received signal $y[k]$ is then just a delayed version of transmitted signal
\begin{equation} y[k] = x[k - k_0] = x[k] * \delta[k - k_0] \end{equation}The aim is to estimate the delay $k_0$ by observation of the received signal under knowledge of the transmitted signal. Under the assumption of finite length signals, the CCF between the transmitted and the received signal is given as
\begin{equation} \begin{split} \varphi_{xy}[\kappa] &= \frac{1}{N} x_N[k] * y_M[-k] \\ &= \frac{1}{N} \left( x_N[k] * x_N[-k] \right) * \delta[-k+k_0] \\ &= \frac{1}{N} \varphi_{xx}[\kappa] * \delta[k - k_0] \end{split} \end{equation}where the associativity of the convolution and the even symmetry of the Dirac pulse has been used. The delay between the two signals results in a shift of the ACF of the transmitted signal. The delay can be estimated by finding the maximum value of the CCF $\varphi_{xy}[\kappa]$, if the ACF $\varphi_{xx}[\kappa]$ of the transmitted signal $x[k]$ exhibits a pronounced maximum at $\kappa = 0$
\begin{equation} \hat{k}_0 = \underset{\kappa}{\mathrm{argmax}} \{ \varphi_{xy}[\kappa] \} \end{equation}One can differentiate between to types of transmitted signals: (i) ones which have been designed explicitly for TOA estimation and (ii) application specific signals. The latter are for instance wireless or speech signals which are used in the applications where specific measurement signals may not be appropriate. For the former, there exist various choices for finite length random sequences with good auto-correlation properties.
The Gold code is a pseudo random noise (PRN) sequence used in the Global Positioning System (GPS) to estimate the TOA between satellites and the receiver. It constitutes a series of binary sequences with good auto- and cross-correlation properties. The latter is required to lower the interference between the sequences of different GPS satellites. The properties of these specific random signals are investigated and the application of these sequences to TOA estimation is illustrated. First the pre-computed sequences are loaded and a function for computing and plotting the CCF (and ACF) is defined for convenience
# load set of PRN sequences (Gold codes)
prn = np.load('../data/gold_sequences.npz')['prn']
def compute_plot_CCF(x, y, ylabel, K=30):
# compute and truncate CCF
ccf = 1/len(x) * np.correlate(x, y, mode='full')
ccf = ccf[(len(y)-1)-K:len(y)+K]
kappa = np.arange(-K, K+1)
# plot CCF
plt.stem(kappa, ccf)
plt.xlabel(r'$\kappa$')
plt.ylabel(ylabel)
plt.axis([-K, K, 1.1*min(ccf), 1.1*max(ccf)]);
plt.grid()
return ccf
The ACF for one particular sequence and the CCF between two different sequences is computed and plotted
plt.figure(figsize = (10, 4))
plt.subplot(121)
compute_plot_CCF(prn[10, :], prn[10, :], r'$\hat{\varphi}_{x_n x_n}[\kappa]$')
plt.title('ACF')
plt.subplot(122)
compute_plot_CCF(prn[10, :], prn[11, :], r'$\hat{\varphi}_{x_n x_m}[\kappa]$')
plt.title('CCF between two PRN sequences')
plt.ylim([-1, 1])
plt.tight_layout()
The ACF features a pronounced peak for $\kappa = 0$ and the two PRN sequences are reasonably uncorrelated. Now one particular sequence is delayed in order to model the received signal. The CCF is plotted and the TOA is estimated by finding the maximum value of the CCF between the transmitted and received sequence
k0 = 10 # true TOA
x = prn[10, :] # pick one PRN sequence
y = x[k0:] # delay transmitted signal by k0 samples
# compute and plot CCF
plt.figure(figsize = (10, 4))
ccf = compute_plot_CCF(x, y, r'$\hat{\varphi}_{xx}[\kappa]$')
plt.title('CCF between transmitted and received signal')
# estimate the TOA
print('Estimated TOA is {:2.0f} samples'.format(np.argmax(ccf) - K))
Estimated TOA is 10 samples
Exercise
Solution: Additive noise can be added for instance by extending above example with y = y + 0.1 * np.random.normal(len(y))
. Additive noise corrupts the CCF. For high levels, the maximum value of the CCF might not correspond to the delay $k_0$ anymore.
The following example estimates and plots the CCF of two random signals $x[k]$ and $y[k]$ of finite lengths $N$ and $M = 2 N$.
N = 1024 # length of random signals
# generate two uncorrelated random signals
np.random.seed(2)
x = 2 + np.random.normal(size=N)
y = 1 + np.random.normal(size=2*N)
# compute CCF
ccf = 1/len(x) * np.correlate(x, y, mode='full')
kappa = np.arange(-(N-1), 2*N)
# print mean values of signals
print('Mean of signal x[k]: %f' %np.mean(x))
print('Mean of signal y[k]: %f' %np.mean(y))
# plot CCF
plt.figure(figsize = (10, 8))
plt.stem(kappa, ccf)
plt.title('Estimated cross-correlation function')
plt.ylabel(r'$\hat{\varphi}_{xy}[\kappa]$')
plt.xlabel(r'$\kappa$')
plt.axis([kappa[0], kappa[-1], 0, 1.1*max(ccf)]);
plt.grid()
Mean of signal x[k]: 1.954480 Mean of signal y[k]: 1.000330
Exercise
Solution: The CCF is approximately constant for $0 < \kappa < N$. Consequently, the random signals $x[k]$ and $y[k]$ can be assumed to be uncorrelated. The trapezoidal shape of the CCF results from the truncation of the random signals to a finite number of samples. This truncation can be modeled as a multiplication of the infinite length signals by a rectangular signal $\text{rect}_N[k]$ and $\text{rect}_M[k]$, respectively. Interpreting the CCF as convolution and having in mind that the convolution of two rectangular signals results in a signal of trapezoidal shape explains the shape of the CCF shown above. Its theoretic value is constant with a value given by the multiplication of the two linear means of the random signals $\varphi_{xy}[\kappa] = \mu_x \cdot \mu_y$.
The output of the CCF were both correlated signals overlap completely, here $0 \leq \kappa < N$, is also termed as valid part of the CCF, the remaining contributions as boundary effects. The computation of the CCF can be limted to the valid part only with the parameter mode=valid
, for instance np.correlate(x, y, mode='valid')
.
The cross-covariance function is defined as the CCF of two zero-mean random signals. For two random signals $x[k]$ and $y[k]$ with linear means $\mu_x \neq 0$ and $\mu_y \neq 0$ it is given as
\begin{equation} \psi_{xy}[\kappa] = \varphi_{xy}[\kappa] - \mu_x \mu_y \end{equation}It follows from the properties of the CCF that the cross-covariance function $\psi_{xy}[\kappa]$ of two uncorrelated signals $x[k]$ and $y[k]$ is zero, $\psi_{xy}[\kappa] = 0$.
Exercise
Copyright
This notebook is provided as Open Educational Resource. Feel free to use the notebook for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2018.