(5 votes, average: 5.00 out of 5)

# Using Toeplitz matrices in MATLAB

by on April 21, 2007

The definition of Toeplitz matrix from [1] is:

A $N\mbox{x}N$ matrix is said to be Toeplitz if the elements $A_{i,j}$ are determined completely by the difference $i-j$.

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 $2N-1$ elements.

For example, the matrix

$A = \left[ \begin{array}{ccc} a_{1,1} & a_{1,2} & a_{1,3} \\ \\ a_{2,1} & a_{1,1} & a_{1,2} \\ \\ a_{3,1} & a_{2,1} & a_{1,1}\end{array}\right]$ 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 $x$ and $h$ is the sum of the product of one function with the time reversed copy of other function i.e.

$x(t)\ast h(t) = \int x(\tau)h(t-\tau) d\tau$ .

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.

$y(n) = \frac{1}{N}\sum_{k=0}^{N-1}x(n-k)$ .

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 $x$ is defined as

$R_{xx}(\tau) = x^{\ast}(-\tau) \ast x(\tau) = \int x(t) x^\ast(t-\tau) dt$.

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 $x$ and $y$is

$R_{xy}({\tau}) = x^{\ast}(-\tau) \ast y(\tau) = \int x(\tau) y^{\ast} (t-\tau) d\tau$.

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:

$y(k) = \sum_{n=0}^{N-1} x_{(k+n)}x^{\ast}_{(k+n+D)}$.

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:

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.

adinda ainaa April 18, 2012 at 2:09 pm

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.

Krishna Sankar April 19, 2012 at 5:02 am

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

sourabh April 2, 2012 at 11:39 am

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

Krishna Sankar April 3, 2012 at 3:48 am

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

tellegant February 10, 2010 at 10:08 pm

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

Krishna Sankar April 4, 2010 at 4:05 am

@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.

shenge86 July 31, 2009 at 1:18 am

How could you use Toeplitz matrices for 2D convolutions?

Krishna Sankar July 31, 2009 at 4:57 am

@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.

Parejas July 3, 2009 at 7:15 am

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