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

D id you like this article? Make sure that you do not miss a new article
by subscribing to RSS feed OR subscribing to e-mail newsletter.
* Note: Subscribing via e-mail entitles you to download the free e-Book on BER of BPSK/QPSK/16QAM/16PSK in AWGN.*

{ 10 comments… read them below or add one }

dear krishna,

please help me..i am really need ur help. i am doing research about SCFDMA and OFDMA, can u give me a Matlab coding for PAPR reduction technique using clipping method and the graph should be BER against SNR for both SCFDMA and OFDMA. The coding should be:

(1) PAPR calculation of SCFDMA without clipping and SCFDMA with clipping.

(2)PAPR calculation of OFDMA without clipping and OFDMA with clipping.

I am going to prove that PAPR will be reduce using clipping method instead of no clipping. The parameter should be :

* 5Mhz system bandwidth

* 512 no. of subcarrier

thank you.

@adinda: I have not tried exploring the PAPR subject much. Will do so in future.

Hi can you shade some light on MMSE in SC-FDMA technique?

@sourabh: will add to the TO-DO’s

hi,

I wonder why are you use y1 = (1/N)*sum(xM); in Implementing moving average sum.Why didnt use y1 = (1/N)*ones(1,N)*xM;

Thank you

@tellegant: Did you mean (1/N)*ones(N,1)*xM ? as long as they are mathematically equivalent, the way we code should not matter much.

How could you use Toeplitz matrices for 2D convolutions?

@shenge86: I have not tried 2D convolutions. However, am guessing that we can use toeplitz matrix and do the convolution on each dimension using for loops.

Hola de parte de parejaspareja.es, encontre tu blog navegando por la red buscando matriz en google. Me parece super interesante la información que tienes en tu blog y sin lugar a dudas regresare a leerlo. Tengo una pregunta, si podria traducir tu blog “Using Toeplitz matrices in MATLAB” y añadirlos a un de mis blogs en italiano? Y por supuesto con el link direccionando a tu blog. Estare esperando tu respuesta. parejaspareja.es

@Parejas: Rather than reproduce the content in your blog, plz provide a link to this page. I will add a translate to Italian option.