Calculating Cross Correlation With MatLAB Code

BASIC

Correlation or is a measure of similarity/ relationship between two signals.

If x[n] & h[n] are two discrete-time signals, then the correlation of x[n] with respect to h[n] is given by,

Correlation mathematically is just Convolution with the second sequence time-reversed.

USAGE

Cross Correlation is necessary to compare one reference signal with one or more signals to determine the similarity between the pair and to determine additional information based on the similarity.

TABLE METHOD

Table method to find cross correlation has he steps below.

• Step 1: List the index ‘k’ covering a sufficient range
• Step 2: List the input x[k]
• Step 3: List the sequence h[k], and align the rightmost element of h[k-n] to the leftmost element of x[k]
• Step 4: Cross-multiply and sum the nonzero overlap terms to produce y[n]
• Step 5: Slide h[k-n] to the right by one position
• Step 6: Repeat step 4; stop if all the output values are zero or if
• required.

Example

Find the cross correlation of the given sequence, x[k]=[3, 1, 2], h[k]=[3, 2, 1] Where x and h both have there first value at 0’th position.

Solution

Start of k = –(length of h – 1) = –(3-1) = – 2

End of k = (length of h + length of x – 1) = (3+3–1) = 5

Now, draw a table. At first row put all values of k from start of k to end of k.

k

–2

–1

0

1

2

3

4

5

At second row put all values of x[k]

x[k]

0

0

3

1

2

0

0

0

Put the values of h(k-(-2)) in third row.

h[k+2]

3

2

1

0

0

0

0

0

Shift the non zero values of h(k+2) to right by 1 bit and put them in fourth row.

h[k+1]

0

3

2

1

0

0

0

0

Similarly for other rows. Shift the values of h(k+2) as long as it is inside the boundary.

h[k]

0

0

3

2

1

0

0

0

h[k-1]

0

0

0

3

2

1

0

0

h[k-2]

0

0

0

0

3

2

1

0

h[k-3]

0

0

0

0

0

3

2

1

Full Table

k	-2	-1	0	1	2	3	4	5
x(k)	0	0	3	1	2	0	0	0
h(k+2)	3	2	1	0	0	0	0	0
h(k+1)	0	3	2	1	0	0	0	0
h(k)	0	0	3	2	1	0	0	0
h(k-1)	0	0	0	3	2	1	0	0
h(k-1)	0	0	0	0	3	2	1	0
h(k-3)	0	0	0	0	0	3	2	1

Now, cross multiply and sum the overlapped terms.

That is multiply each term of x(k) with respective term of h(k+2) and sum them to find X(-2). (N.B. X(n)).

Similarly multiply and sum h(k+1), h(k), h(k-1) with x(k) to get X(-1), X(0), X(1) etc.

Here,

X(-2)=3×1=3

X(-1)=3×2+1×1=7

X(0)=3×3+1×2+2×1=13

X(1)=1×3+2×2=7

X(2)=2×3=6

X(3=0 [No overlap]

Hence, X(n) = [3, 7, 13, 7, 6, 0]

MATLAB CODE

To calculate cross correlation first of all we need to know the value of x(n) and h(n). Moreover we need to know which one value of them is located at zero position. We will ask for this to the user.

x = input('Enter values of x(n): ');
zerox = input('Zero Position of x: ');
h = input('Enter values of h(n): ');
zeroh = input('Zero Position of h: ');

Then we need to know the starting value of k and end value of k. Carefully listen that we will say these values of k as actual index in later parts of the tutorial. We will find the start value of k and end value of k from the formula above.

k1 = -(length(h)-1);
k2 = length(x)+length(h)-1;

If we see here the values of k is started from negative. But MatLAB cannot take negative index. So we will take a positive k where the value of k will be started from 1 and will run till (k2-k1+1). That is the positive part, equal part of negative and 1 for 0 position. See the fraction.

kstart = 1;
kend = k2 - k1 + 1;

As we need to multiply value of x and values of h so they cannot be empty. Empty array or variable cause an error in matlab. So, we will declare xk and hk array with 0 from index kstart to kend (positive index).

xk = zeros(kstart,kend);
hk = zeros(kstart,kend);

Now we will put the values of x to xk with a loop as the zero positioned value of x is placed at (-k1+1) position. That is our current zero in positive index respective to actual index.

for i=1:length(x)
  xk(-k1+i-zerox+1) = x(i);
endfor

Put the values of h to the array hk as the zero positioned value is currently at (-k1+1) position(current zero).

for i=1:length(h)
  hk(i+zeroh-1) = h(i);
endfor

Now take a variable shift that will count the shifting of h(k-n). Here shift denote the n.

shift = 0;

Now run two loop where one is nested in another. One loop (external) will find the values of X(m). The loop inside will find the sum of h(k) and x(k).

The external loop will run from 1 to (k2-zeroh+1+1). Here k2 is the last limit of our actual index. zeroh is the position at what h denote index zero. zeroh is subtracted because when we flip h zeroh bit goes to right side of the array. And that cross the boundary quickly from x(k). One 1 is to cover zeroh and another is to cover k2.

We set X(shift+1) to 0 at each turn because it is added as previous value (summed). So, it cannot be empty.

The internal for loop runs from 1 to (length(h)+zeroh-1). Because we have no need to multiply and add values where h(k)=0. We assign the summed value in X(shift+1). The cross multiplication at each turn of this loop is hk(i)*xk(i+shift). That is for each turn we don’t shifting h(-k) to right rather we are shifting x(k) to left. And that is same.

You need to increase the value of shift as index and shifting size increase. That will be done after the internal loop.

for j=1:k2-zeroh+1+1
  X(shift+1) = 0;
  for i=1:length(h)+zeroh-1
    X(shift+1) = X(shift+1) + hk(i)*xk(i+shift);
  endfor
  shift = shift+1;
endfor

Now, Everything is done. Display the value with corresponding title.

disp('Value of x(k): ');
disp(xk);
disp('Value of h(-k): ');
disp(hk);
disp('Convolution Sum: ');
disp(X);

Plot the result at discrete format with respective title at (3,1) subplot.

subplot(3,1,1);
stem(1:length(xk),xk);
title('x(n)');
subplot(3,1,2);
stem(1:length(hk),hk);
title('h(n)');
subplot(3,1,3);
stem(1:length(X),X);
title('X(m)');

FULL CODE

x = input('Enter values of x(n): ');
zerox = input('Zero Position of x: ');
h = input('Enter values of h(n): ');
zeroh = input('Zero Position of h: ');
k1 = -(length(h)-1);
k2 = length(x)+length(h)-1;
kstart = 1;
kend = k2 - k1 + 1;
xk = zeros(kstart,kend);
hk = zeros(kstart,kend);
for i=1:length(x)
  xk(-k1+i-zerox+1) = x(i);
endfor
h = flip(h);
for i=1:length(h)
  hk(i+zeroh-1) = h(i);
endfor
shift = 0;
for j=1:k2-zeroh+1+1
  X(shift+1) = 0;
  for i=1:length(h)+zeroh-1
    X(shift+1) = X(shift+1) + hk(i)*xk(i+shift);
  endfor
  shift = shift+1;
endfor
disp('Value of x(k): ');
disp(xk);
disp('Value of h(-k): ');
disp(hk);
disp('Convolution Sum: ');
disp(X);
subplot(3,1,1);
stem(1:length(xk),xk);
title('x(n)');
subplot(3,1,2);
stem(1:length(hk),hk);
title('h(n)');
subplot(3,1,3);
stem(1:length(X),X);
title('X(m)');

I think your experience with tutorial is quite good. If you have to seek anything ask in comment or report by form.

Thanks.