Accelerating the pace of engineering and science

# xcorr

Cross-correlation

## Description

example

r = xcorr(x,y) returns the cross-correlation of two discrete-time sequences, x and y. Cross-correlation measures the similarity between x and shifted (lagged) copies of y as a function of the lag. If x and y have different lengths, the function appends zeros at the end of the shorter vector so it has the same length as the other.

example

r = xcorr(x) returns the autocorrelation sequence of x. If x is a matrix, then r is a matrix whose columns contain the autocorrelation and cross-correlation sequences for all combinations of the columns of x.

example

r = xcorr(___,maxlag) limits the lag range from –maxlag to maxlag. This syntax accepts one or two input sequences. maxlag defaults to N – 1.

example

r = xcorr(___,scaleopt) additionally specifies a normalization option for the cross-correlation or autocorrelation. Any option other than 'none' (the default) requires x and y to have the same length.

example

[r,lags] = xcorr(___) also returns a vector with the lags at which the correlations are computed.

## Examples

expand all

### Delay Between Two Correlated Signals

Two sensors at different locations measure vibrations caused by a car as it crosses a bridge.

Load the signals and the sample rate, . Create time vectors and plot the signals. The signal from Sensor 2 arrives at an earlier time than the signal from Sensor 1.

load sensorData

t1 = (0:length(s1)-1)/Fs;
t2 = (0:length(s2)-1)/Fs;

subplot(2,1,1)
plot(t1,s1)
title('s_1')

subplot(2,1,2)
plot(t2,s2)
title('s_2')
xlabel('Time (s)')


The cross-correlation of the two measurements is maximum at a lag equal to the delay.

Plot the cross-correlation. Express the delay as a number of samples and in seconds.

[acor,lag] = xcorr(s2,s1);

[~,I] = max(abs(acor));
lagDiff = lag(I)
timeDiff = lagDiff/Fs

figure
plot(lag,acor)
a3 = gca;
a3.XTick = sort([-3000:1000:3000 lagDiff]);

lagDiff =

-350

timeDiff =

-0.0317



Align the two signals and replot them.

s1al = s1(-lagDiff:end);
t1al = (0:length(s1al)-1)/Fs;

subplot(2,1,1)
plot(t1al,s1al)
title('s_1, aligned')

subplot(2,1,2)
plot(t2,s2)
title('s_2')
xlabel('Time (s)')


### Echo Cancelation

A speech recording includes an echo caused by reflection off a wall. Use autocorrelation to filter it out.

In the recording, a person says the word MATLAB®. Load the data and the sample rate, .

load mtlb

% To hear, type soundsc(mtlb,Fs)


Model the echo by adding to the recording a copy of the signal delayed by samples and attenuated by a known factor : . Specify a time lag of 0.23 s and an attenuation factor of 0.5.

timelag = 0.23;
delta = round(Fs*timelag);
alpha = 0.5;

orig = [mtlb;zeros(delta,1)];
echo = [zeros(delta,1);mtlb]*alpha;

mtEcho = orig + echo;


Plot the original, the echo, and the resulting signal.

t = (0:length(mtEcho)-1)/Fs;

subplot(2,1,1)
plot(t,[orig echo])
legend('Original','Echo')

subplot(2,1,2)
plot(t,mtEcho)
legend('Total')
xlabel('Time (s)')

% To hear, type soundsc(mtEcho,Fs)


Compute an unbiased estimate of the signal autocorrelation. Select and plot the section that corresponds to lags greater than zero.

[Rmm,lags] = xcorr(mtEcho,'unbiased');

Rmm = Rmm(lags>0);
lags = lags(lags>0);

figure
plot(lags/Fs,Rmm)
xlabel('Lag (s)')


The autocorrelation has a sharp peak at the lag at which the echo arrives. Cancel the echo by filtering the signal through an IIR system whose output, , obeys .

[~,dl] = findpeaks(Rmm,lags,'MinPeakHeight',0.22);

mtNew = filter(1,[1 zeros(1,dl-1) alpha],mtEcho);


Plot the filtered signal and compare to the original.

subplot(2,1,1)
plot(t,orig)
legend('Original')

subplot(2,1,2)
plot(t,mtNew)
legend('Filtered')
xlabel('Time (s)')

% To hear, type soundsc(mtNew,Fs)


### Cross-Correlation with Multichannel Input

Generate three 11-sample exponential sequences given by , , and , with . Use stem3 to plot the sequences side by side.

N = 11;
n = (0:N-1)';

a = 0.4;
b = 0.7;
c = 0.999;

xabc = [a.^n b.^n c.^n];

stem3(n,1:3,xabc','filled')
ax = gca;
ax.YTick = 1:3;
view(37.5,30)


Compute the autocorrelations and mutual cross-correlations of the sequences. Output the lags so you do not have to keep track of them. Normalize the result so the autocorrelations have unit value at zero lag.

[cr,lgs] = xcorr(xabc,'coeff');

for row = 1:3
for col = 1:3
nm = 3*(row-1)+col;
subplot(3,3,nm)
stem(lgs,cr(:,nm),'.')
title(sprintf('c_{%d%d}',row,col))
ylim([0 1])
end
end


Restrict the calculation to lags between and .

[cr,lgs] = xcorr(xabc,5,'coeff');

for row = 1:3
for col = 1:3
nm = 3*(row-1)+col;
subplot(3,3,nm)
stem(lgs,cr(:,nm),'.')
title(sprintf('c_{%d%d}',row,col))
ylim([0 1])
end
end


Compute unbiased estimates of the autocorrelations and mutual cross-correlations. By default, the lags run between and .

cu = xcorr(xabc,'unbiased');

for row = 1:3
for col = 1:3
nm = 3*(row-1)+col;
subplot(3,3,nm)
stem(-(N-1):(N-1),cu(:,nm),'.')
title(sprintf('c_{%d%d}',row,col))
end
end


### Autocorrelation Function of Exponential Sequence

Compute the autocorrelation function of a 28-sample exponential sequence, for .

a = 0.95;

N = 28;
n = 0:N-1;
lags = -(N-1):(N-1);

x = a.^n;
c = xcorr(x);


Determine analytically to check the correctness of the result. Use a larger sample rate to simulate a continuous situation. The autocorrelation function of the sequence for , with , is

fs = 10;
nn = -(N-1):1/fs:(N-1);
cc = a.^abs(nn)/(1-a^2);

dd = (1-a.^(2*(N-abs(nn))))/(1-a^2).*a.^abs(nn);


Plot the sequences on the same figure.

stem(lags,c);
hold on
plot(nn,dd)
xlabel('Lag')
legend('xcorr','Analytic')
hold off


Repeat the calculation, but now find an unbiased estimate of the autocorrelation. Verify that the unbiased estimate is given by .

cu = xcorr(x,'unbiased');

du = dd./(N-abs(nn));

stem(lags,cu);
hold on
plot(nn,du)
xlabel('Lag')
legend('xcorr','Analytic')
hold off


Repeat the calculation, but now find a biased estimate of the autocorrelation. Verify that the biased estimate is given by .

cb = xcorr(x,'biased');

db = dd/N;

stem(lags,cb);
hold on
plot(nn,db)
xlabel('Lag')
legend('xcorr','Analytic')
hold off


Find an estimate of the autocorrelation whose value at zero lag is unity.

cz = xcorr(x,'coeff');

dz = dd/max(dd);

stem(lags,cz);
hold on
plot(nn,dz)
xlabel('Lag')
legend('xcorr','Analytic')
hold off


### Cross-Correlation of Two Exponential Sequences

Compute and plot the cross-correlation of two 16-sample exponential sequences, and , with .

N = 16;
n = 0:N-1;

a = 0.84;
b = 0.92;

xa = a.^n;
xb = b.^n;

r = xcorr(xa,xb);

stem(-(N-1):(N-1),r)


Determine analytically to check the correctness of the result. Use a larger sample rate to simulate a continuous situation. The cross-correlation function of the sequences and for , with , is

fs = 10;
nn = -(N-1):1/fs:(N-1);

cn = (1 - (a*b).^(N-abs(nn)))/(1 - a*b) .* ...
(a.^nn.*(nn>0) + (nn==0) + b.^-(nn).*(nn<0));


Plot the sequences on the same figure.

hold on
plot(nn,cn)

xlabel('Lag')
legend('xcorr','Analytic')


Verify that switching the order of the operands reverses the sequence.

figure

stem(-(N-1):(N-1),xcorr(xb,xa))

hold on
stem(-(N-1):(N-1),fliplr(r),'--*')

xlabel('Lag')
legend('xcorr(x_b,x_a)','fliplr(xcorr(x_a,x_b))')


Generate the 20-sample exponential sequence . Compute and plot its cross-correlations with and . Output the lags to make the plotting easier. xcorr appends zeros at the end of the shorter sequence to match the length of the longer one.

xc = 0.77.^(0:20-1);

[xca,la] = xcorr(xa,xc);
[xcb,lb] = xcorr(xb,xc);

figure

subplot(2,1,1)
stem(la,xca)
subplot(2,1,2)
stem(lb,xcb)
xlabel('Lag')


### GPU Acceleration for Autocorrelation Sequence Estimation

This example requires Parallel Computing Toolbox™ software and a CUDA-enabled NVIDIA GPU with compute capability 1.3 or above. See GPU System Requirements for details.

Create a signal consisting of a 10 Hz sine wave in additive noise, sampled at 1 kHz. Use gpuArray to create a gpuArray object stored on your computer's GPU.

t = 0:0.001:10-0.001;
x = cos(2*pi*10*t) + randn(size(t));
X = gpuArray(x);

Compute the normalized autocorrelation sequence to lag 200.

[xc,lags] = xcorr(X,200,'coeff');

The output, xc, is a gpuArray object.

Use gather to transfer the data from the GPU to the MATLAB® workspace as a double-precision vector.

xc = gather(xc);

## Input Arguments

expand all

### x — Input arrayvector | matrix | gpuArray object

Input array, specified as a vector, a matrix, or a gpuArray object.

See GPU Computing and GPU System Requirements for details on using xcorr with gpuArray objects.

Example: sin(2*pi*(0:9)/10) + randn([1 10])/10 specifies a noisy sinusoid as a row vector.

Example: sin(2*pi*[0.1;0.3]*(0:39))' + randn([40 2])/10 specifies a two-channel noisy sinusoid.

Example: gpuArray(sin(2*pi*(0:9)/10) + randn([1 10])/10) specifies a noisy sinusoid as a gpuArray object.

Data Types: single | double
Complex Number Support: Yes

### y — Input arrayvector | gpuArray object

Input array, specified as a vector or a gpuArray object.

Data Types: single | double
Complex Number Support: Yes

### maxlag — Maximum laginteger scalar

Maximum lag, specified as an integer scalar. If you specify maxlag, the returned cross-correlation sequence ranges from –maxlag to maxlag. If you do not specify maxlag, the lag range equals 2N – 1, where N is the greater of the lengths of x and y.

Data Types: single | double

### scaleopt — Normalization option'none' (default) | 'biased' | 'unbiased' | 'coeff'

Normalization option, specified as one of these strings:

• 'none' — Raw, unscaled cross-correlation. This is the only allowed option when x and y have different lengths.

• 'biased' — Biased estimate of the cross-correlation:

${\stackrel{^}{R}}_{xy,\text{biased}}\left(m\right)=\frac{1}{N}{\stackrel{^}{R}}_{xy}\left(m\right).$

• 'unbiased' — Unbiased estimate of the cross-correlation:

${\stackrel{^}{R}}_{xy,\text{unbiased}}\left(m\right)=\frac{1}{N-|m|}{\stackrel{^}{R}}_{xy}\left(m\right).$

• 'coeff' — Normalizes the sequence so that the autocorrelations at zero lag equal 1.

Data Types: char

## Output Arguments

expand all

### r — Cross-correlation or autocorrelation sequencevector | matrix | gpuArray object

Cross-correlation sequence, returned as a vector, a matrix, or a gpuArray object.

If x is an M × N signal matrix representing N channels in its columns, then xcorr(x) returns a (2M – 1) × N2 matrix with the autocorrelations and mutual cross-correlations of the channels of x. If you specify maxlag, then r has size (2 × maxlag – 1) × N2.

For example, if S is a three-channel signal, $\text{S}=\left(\begin{array}{ccc}{x}_{1}& {x}_{2}& {x}_{3}\end{array}\right)$, then the result of R = xcorr(S) is organized as

$\text{R}=\left(\begin{array}{lllllllll}{R}_{{x}_{1}{x}_{1}}\hfill & {R}_{{x}_{1}{x}_{2}}\hfill & {R}_{{x}_{1}{x}_{3}}\hfill & {R}_{{x}_{2}{x}_{1}}\hfill & {R}_{{x}_{2}{x}_{2}}\hfill & {R}_{{x}_{2}{x}_{3}}\hfill & {R}_{{x}_{3}{x}_{1}}\hfill & {R}_{{x}_{3}{x}_{2}}\hfill & {R}_{{x}_{3}{x}_{3}}\hfill \end{array}\right).$

### lags — Lag indicesvector

Lag indices, returned as a vector.

expand all

### Cross-Correlation and Autocorrelation

The result of xcorr can be interpreted as an estimate of the correlation between two random sequences or as the deterministic correlation between two deterministic signals.

The true cross-correlation sequence of two jointly stationary random processes, xn and yn, is given by

${R}_{xy}\left(m\right)=E\left\{{x}_{n+m}{y}_{n}^{*}\right\}=E\left\{{x}_{n}{y}_{n-m}^{*}\right\},$

where − < n < , the asterisk denotes complex conjugation, and E is the expected value operator. xcorr can only estimate the sequence because, in practice, only a finite segment of one realization of the infinite-length random process is available.

By default, xcorr computes raw correlations with no normalization:

${\stackrel{^}{R}}_{xy}\left(m\right)=\left\{\begin{array}{ll}\sum _{n=0}^{N-m-1}{x}_{n+m}{y}_{n}^{\ast },\hfill & m\ge 0,\hfill \\ {\stackrel{^}{R}}_{yx}^{*}\left(-m\right),\hfill & m<0.\hfill \end{array}$

The output vector, c, has elements given by

$\text{c}\left(\text{m}\right)={\stackrel{^}{R}}_{xy}\left(m-N\right),\text{ }\text{ }m=1,\text{\hspace{0.17em}}2,\dots ,2N-1.$

In general, the correlation function requires normalization to produce an accurate estimate. You can control the normalization of the correlation by using the input argument scaleopt.

## References

[1] Buck, John R., Michael M. Daniel, and Andrew C. Singer. Computer Explorations in Signals and Systems Using MATLAB. 2nd Edition. Upper Saddle River, NJ: Prentice Hall, 2002.

[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.