Using Toeplitz matrices in MATLAB

The definition of Toeplitz matrix from [1] is:

A matrix is said to be Toeplitz if the elements are determined completely by the difference .

Pictorially, if a line is drawn parallel to the main diagonal then all the elements on this line are equal.

A Toeplitz matrix is completely determined by the 1st row and 1st column of the matrix i.e. by elements.

For example, the matrix

is a Toeplitz matrix.

With this understanding, let us move on to some useful examples in Matlab where the Toeplitz matrix construction is used to implement some standard functions.

Implementing convolution

Convolution of two functions and is the sum of the product of one function with the time reversed copy of other function i.e.

.

Typically, one can visualize implementing this using a for-loop OR use the conv() function in MATLAB. However, if you don’t want to use the conv() operator and do not want the for-loop, you can still implement convolution by as a matrix multiplication by constructing one of input as a Toeplitz matrix.

% implementing convolution
x = [1:10]; % first sequence
h = [1 1 -1 2 3 1 -1]; % second sequence

% Forming the toeplitz matrix from x
% each column of the matrix stores the values of x as they slide in through the tapped delay line.

% the number of elements in the tapped delay line is equal to the number
% of elements in h

xM = toeplitz([x(1) zeros(1,length(h)-1) ], [x zeros(1,length(h)-1) ]);

y1 = conv(x,h); % convolution output

y2 = h*xM; % convolution using toeplitz matrix

diff = (y2 – y1)

Implementing moving average sum

Moving average sum of stores the mean of the latest N samples i.e.

.

Infact, the N point moving average sum can be implemented as convolution of the input sequence with a matrix having N ones. Since, the operation involves convolution, the same can be implemented using a Toeplitz matrix structure for the input.

% implementing Moving Average sum of the latest N samples
x = [1:10]; % input sequence
N = 5; % number of samples to be accumulated
% forming the toeplitz matrix from x
% each column of the matriz stores the latest N samples
xM = toeplitz([x(1) zeros(1,N-1) ], [x ]);
y1 = (1/N)*sum(xM);

y2 = conv(x,ones(1,N))/N; % implementing MA as convolutions

diff = y1 – y2(1:length(x))

Implementing auto correlation of the signal

The auto-correlation of the signal is defined as

.

Obviously as this involves convolution, the operation can be implemented as the product of two matrices, where one of them is a Toeplitz matrix.

x = cos(2*pi*1/32*[1:320]); % input sequence
N = length(x)

% padded x with N-1 zeros on both sides
xM = [ zeros(1,N-1) x zeros(1,N-1)];
% modified x to form a toeplitx matrix. The first column of x stores
% the elements of x followed by 2*N-2 zeros. The second column stores
% zero followed by x followed by 2*N-1 zeros and so on…
% the matrix stores x
xDelay_M =[ toeplitz([x zeros(1,2*N-2)], [ x(1) zeros(1,2*N-2) ] ) ];

y1 = xM*conj(xDelay_M);
y2 = xcorr(x);

% mean square error
diff = (y2-y1)*(y2-y1)’/length(y2-y1)

As can be observed, the output of the Toeplitx matrix implementation matches the xcorr() function provided in Matlab.

Implementing cross correlation of the signal

Once auto-correlation implementation is shown to be correct, it might be reasonably easy to extend it for the cross correlation case. The cross-correlation of two sequence and is

.

This can be implemented using the same steps performed for computing auto-correlation. However, to address the case where the size of the input sequence differ, the shorter sequence is padded with zeros.

x = 10*cos(2*pi*1/32*[1:320]) ;
y = randn(1,310) + j*randn(1,310)
nX = length(x);
nY = length(y);

% when the length of x, y are different, padding the smaller matrix with
% zeros
m = max(nX,nY);
x1 = zeros(1,m);
y1 = zeros(1,m);
x1(1:length(x)) = x;
y1(1:length(y)) = y;

% padded x with m-1 zeros on both sides
xM = [ zeros(1,m-1) x1 zeros(1,m-1)];
% modified y1 to form a toeplitx matrix.
yDelay_M =[ toeplitz([y1 zeros(1,2*m-2)],[ y1(1) zeros(1,2*m-2) ] ) ];
op1 = xM*conj(yDelay_M);

op2 = xcorr(x,y);
diff = (op1-op2)*(op1-op2)’/length(op1-op2);

It can be observed that the output of the xcorr() function in Matlab matches with the output from the Toeplitz implementation.

Implementing auto-correlation of the signal delayed by D samples and the accumulated by N samples

Typically, for identifying noisy sequence with a known periodic properties, it might be desirable to perform autocorrelation for a given delay D accumulated for given number of samples N. The operation is:

.

The sequence x is converted into two Toeplitz matrices – one storing the current samples and the other storing the delayed samples. The output is computed by taking the column sum of dot product multiplication of the matrix storing the current samples with the conjugate of the delayed sample matrix.

ip = cos(2*pi*1/10*[1:10]);

x = [ ip ip ip ip ip zeros(1,40) ] ; % input sequence
D = 10; % delay
N = 20; % accumulation

% forming the toeplitz matrix storing the delayed elements
% The first D columns of the matrix is all zeros. From D+1 column onwards,
% the elements on x starts sliding in. Each column is of length N
xDelay_M = toeplitz([zeros(1,N) ], [zeros(1,D) x zeros(1,N-1)]);

% forming the toeplitz matrix storing the current elements
% From 1st column onwards, the elements on x starts sliding in. Each column
% is of length N. The last D columns are all zeros.
xM = toeplitz([x(1) zeros(1,N-1) ], [x zeros(1,D + N-1)]);

% dot product of the current and the conjugate of the delayed sequence
y1 = sum(conj(xDelay_M).*xM);
plot(abs(y1));

As expected it can be seen that, due to the D delayed periodicity of the input, the auto-correlation output rises and hits the plateau value at the (D+N)th sample. The autocorrelation output starts falling from the plateau value at the end of the periodic seqeunce, which in this example is the 50th sample.

References:

[1] Multirate Systems And Filter Banks, P. P. Vaidyanathan

Technorati tags: Toeplitz matrix