Correlation Analysis
The correlation operator (defined in §7.2.5) plays a major role in statistical signal processing. For a proper development, see, e.g., [27,33,65]. This section introduces only some of the most basic elements of statistical signal processing in a simplified manner, with emphasis on illustrating applications of the DFT.
Cross-Correlation
Definition: The circular cross-correlation of two signals and
in
may be defined by
![$\displaystyle \zbox {{\hat r}_{xy}(l) \isdef \frac{1}{N}(x\star y)(l)
\isdef \frac{1}{N}\sum_{n=0}^{N-1}\overline{x(n)} y(n+l), \; l=0,1,2,\ldots,N-1.}
$](http://www.dsprelated.com/josimages_new/mdft/img1550.png)
![$ l$](http://www.dsprelated.com/josimages_new/mdft/img1446.png)
![$ 1$](http://www.dsprelated.com/josimages_new/mdft/img111.png)
![$ \star$](http://www.dsprelated.com/josimages_new/mdft/img1406.png)
The term ``cross-correlation'' comes from
statistics, and what we have defined here is more properly
called a ``sample cross-correlation.''
That is,
is an
estimator8.8 of the true
cross-correlation
which is an assumed statistical property
of the signal itself. This definition of a sample cross-correlation is only valid for
stationary stochastic processes, e.g., ``steady noises'' that
sound unchanged over time. The statistics of a stationary stochastic
process are by definition time invariant, thereby allowing
time-averages to be used for estimating statistics such
as cross-correlations. For brevity below, we will typically
not include ``sample'' qualifier, because all computational
methods discussed will be sample-based methods intended for use on
stationary data segments.
The DFT of the cross-correlation may be called the cross-spectral density, or ``cross-power spectrum,'' or even simply ``cross-spectrum'':
![$\displaystyle {\hat R}_{xy}(\omega_k) \isdef \hbox{\sc DFT}_k({\hat r}_{xy}) = \frac{\overline{X(\omega_k)}Y(\omega_k)}{N}
$](http://www.dsprelated.com/josimages_new/mdft/img1553.png)
Unbiased Cross-Correlation
Recall that the cross-correlation operator is cyclic (circular)
since is interpreted modulo
. In practice, we are normally
interested in estimating the acyclic cross-correlation
between two signals. For this (more realistic) case, we may define
instead the unbiased cross-correlation
![$\displaystyle \zbox {{\hat r}^u_{xy}(l) \isdef \frac{1}{N-l}\sum_{n=0}^{N-1-l} \overline{x(n)} y(n+l),\quad
l = 0,1,2,\ldots,L-1}
$](http://www.dsprelated.com/josimages_new/mdft/img1555.png)
![$ L\ll N$](http://www.dsprelated.com/josimages_new/mdft/img1556.png)
![$ L\approx\sqrt{N}$](http://www.dsprelated.com/josimages_new/mdft/img1557.png)
![$ \overline{x(n)} y(n+l)$](http://www.dsprelated.com/josimages_new/mdft/img1558.png)
![$ L-1$](http://www.dsprelated.com/josimages_new/mdft/img1218.png)
![$ n=N-l-1$](http://www.dsprelated.com/josimages_new/mdft/img1559.png)
![$ n$](http://www.dsprelated.com/josimages_new/mdft/img80.png)
![$ N$](http://www.dsprelated.com/josimages_new/mdft/img35.png)
![$ {\hat r}^u_{xy}(l)$](http://www.dsprelated.com/josimages_new/mdft/img1560.png)
![$ r_{xy}(l)$](http://www.dsprelated.com/josimages_new/mdft/img1552.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
![$ y$](http://www.dsprelated.com/josimages_new/mdft/img26.png)
An unbiased acyclic cross-correlation may be computed faster via DFT (FFT) methods using zero padding:
![$\displaystyle \zbox {{\hat r}^u_{xy}(l) = \frac{1}{N-l}\hbox{\sc IDFT}_l(\overline{X}\cdot Y), \quad
l = 0,1,2,\ldots,L-1}
$](http://www.dsprelated.com/josimages_new/mdft/img1561.png)
![\begin{eqnarray*}
X &=& \hbox{\sc DFT}[\hbox{\sc CausalZeroPad}_{N+L-1}(x)]\\
Y &=& \hbox{\sc DFT}[\hbox{\sc CausalZeroPad}_{N+L-1}(y)].
\end{eqnarray*}](http://www.dsprelated.com/josimages_new/mdft/img1562.png)
Note that and
belong to
while
and
belong to
. The zero-padding may be causal (as defined in
§7.2.8)
because the signals are assumed to be be stationary, in which case all
signal statistics are time-invariant. As usual when embedding acyclic
correlation (or convolution) within the cyclic variant given by the
DFT, sufficient zero-padding is provided so that only zeros are ``time
aliased'' (wrapped around in time) by modulo indexing.
Cross-correlation is used extensively in audio signal processing for applications such as time scale modification, pitch shifting, click removal, and many others.
Autocorrelation
The cross-correlation of a signal with itself gives its autocorrelation:
![$\displaystyle \zbox {{\hat r}_x(l) \isdef \frac{1}{N}(x\star x)(l)
\isdef \frac{1}{N}\sum_{n=0}^{N-1}\overline{x(n)} x(n+l)}
$](http://www.dsprelated.com/josimages_new/mdft/img1564.png)
![$\displaystyle {\hat r}_x(-l) = \overline{{\hat r}_x(l)}
$](http://www.dsprelated.com/josimages_new/mdft/img1565.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
The unbiased cross-correlation similarly reduces to an unbiased
autocorrelation when :
The DFT of the true autocorrelation function
is the (sampled)
power spectral density (PSD), or power spectrum, and may
be denoted
![$\displaystyle R_x(\omega_k) \isdef \hbox{\sc DFT}_k(r_x).
$](http://www.dsprelated.com/josimages_new/mdft/img1569.png)
![$ R_x(\omega) \isdef
\hbox{\sc DTFT}_k(r_x)$](http://www.dsprelated.com/josimages_new/mdft/img1570.png)
![$ {\hat r}_x$](http://www.dsprelated.com/josimages_new/mdft/img1571.png)
![$\displaystyle {\hat R}_x(\omega_k)=\hbox{\sc DFT}_k({\hat r}_x) = \frac{\left\vert X(\omega_k)\right\vert^2}{N}
$](http://www.dsprelated.com/josimages_new/mdft/img1572.png)
![$ {\hat R}_x(\omega_k)$](http://www.dsprelated.com/josimages_new/mdft/img1573.png)
At lag zero, the autocorrelation function reduces to the average power (mean square) which we defined in §5.8:
![$\displaystyle {\hat r}_x(0) \isdef \frac{1}{N}\sum_{m=0}^{N-1}\left\vert x(m)\right\vert^2 % \isdef \Pscr_x^2
$](http://www.dsprelated.com/josimages_new/mdft/img1574.png)
Replacing ``correlation'' with ``covariance'' in the above definitions gives corresponding zero-mean versions. For example, we may define the sample circular cross-covariance as
![$\displaystyle \zbox {{\hat c}_{xy}(n)
\isdef \frac{1}{N}\sum_{m=0}^{N-1}\overline{[x(m)-\mu_x]} [y(m+n)-\mu_y].}
$](http://www.dsprelated.com/josimages_new/mdft/img1575.png)
![$ \mu_x$](http://www.dsprelated.com/josimages_new/mdft/img743.png)
![$ \mu_y$](http://www.dsprelated.com/josimages_new/mdft/img1576.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
![$ y$](http://www.dsprelated.com/josimages_new/mdft/img26.png)
![$ {\hat c}_x(0)$](http://www.dsprelated.com/josimages_new/mdft/img1577.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
![$\displaystyle {\hat c}_x(0) \isdef \frac{1}{N}\sum_{m=0}^{N-1}\left\vert x(m)-\mu_x\right\vert^2 \isdef {\hat \sigma}_x^2
$](http://www.dsprelated.com/josimages_new/mdft/img1578.png)
Matched Filtering
The cross-correlation function is used extensively in pattern
recognition and signal detection. We know from Chapter 5
that projecting one signal onto another is a means of measuring how
much of the second signal is present in the first. This can be used
to ``detect'' the presence of known signals as components of more
complicated signals. As a simple example, suppose we record
which we think consists of a signal
that we are looking for
plus some additive measurement noise
. That is, we assume the
signal model
. Then the projection of
onto
is
(recalling §5.9.9)
![$\displaystyle {\bf P}_s(x) \isdef \frac{\left<x,s\right>}{\Vert s\Vert^2} s
= \...
...}{\Vert s\Vert^2} s
= s + \frac{N}{\Vert s\Vert^2} {\hat r}_{se}(0)s
\approx s
$](http://www.dsprelated.com/josimages_new/mdft/img1582.png)
![$ e$](http://www.dsprelated.com/josimages_new/mdft/img92.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ \hbox{\sc Flip}(s)$](http://www.dsprelated.com/josimages_new/mdft/img1583.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ \hbox{\sc Flip}(\overline{s})$](http://www.dsprelated.com/josimages_new/mdft/img1584.png)
![$ s$](http://www.dsprelated.com/josimages_new/mdft/img29.png)
![$ x$](http://www.dsprelated.com/josimages_new/mdft/img25.png)
![$ {\hat r}_{sx}(l)$](http://www.dsprelated.com/josimages_new/mdft/img1585.png)
In the same way that FFT convolution is faster than direct convolution (see Table 7.1), cross-correlation and matched filtering are generally carried out most efficiently using an FFT algorithm (Appendix A).
FIR System Identification
Estimating an impulse response from input-output measurements is called system identification, and a large literature exists on this topic (e.g., [39]).
Cross-correlation can be used to compute the impulse response
of a filter from the cross-correlation of its input and output signals
and
, respectively. To see this, note that, by
the correlation theorem,
![$\displaystyle x\star y \;\longleftrightarrow\;\overline{X}\cdot Y = \overline{X}\cdot (H\cdot X) =
H\cdot\left\vert X\right\vert^2.
$](http://www.dsprelated.com/josimages_new/mdft/img1587.png)
![$\displaystyle H = \frac{\overline{X}\cdot Y}{\left\vert X\right\vert^2} = \frac{{\hat R}_{xy}}{{\hat R}_{xx}}
$](http://www.dsprelated.com/josimages_new/mdft/img1588.png)
![$ H(\omega_k) = \overline{X(\omega_k)}\cdot Y(\omega_k)/\vert X(\omega_k)\vert^2$](http://www.dsprelated.com/josimages_new/mdft/img1589.png)
% sidex.m - Demonstration of the use of FFT cross- % correlation to compute the impulse response % of a filter given its input and output. % This is called "FIR system identification". Nx = 32; % input signal length Nh = 10; % filter length Ny = Nx+Nh-1; % max output signal length % FFT size to accommodate cross-correlation: Nfft = 2^nextpow2(Ny); % FFT wants power of 2 x = rand(1,Nx); % input signal = noise %x = 1:Nx; % input signal = ramp h = [1:Nh]; % the filter xzp = [x,zeros(1,Nfft-Nx)]; % zero-padded input yzp = filter(h,1,xzp); % apply the filter X = fft(xzp); % input spectrum Y = fft(yzp); % output spectrum Rxx = conj(X) .* X; % energy spectrum of x Rxy = conj(X) .* Y; % cross-energy spectrum Hxy = Rxy ./ Rxx; % should be the freq. response hxy = ifft(Hxy); % should be the imp. response hxy(1:Nh) % print estimated impulse response freqz(hxy,1,Nfft); % plot estimated freq response err = norm(hxy - [h,zeros(1,Nfft-Nh)])/norm(h); disp(sprintf(['Impulse Response Error = ',... '%0.14f%%'],100*err)); err = norm(Hxy-fft([h,zeros(1,Nfft-Nh)]))/norm(h); disp(sprintf(['Frequency Response Error = ',... '%0.14f%%'],100*err)); |
Next Section:
Power Spectral Density Estimation
Previous Section:
Filters and Convolution