Monthly Archives: November 2019

Circular Convolution, Discrete Fourier Transforms and Toeplitz Matrix Multiplication

Circular Convolution

It’s fairly well known that multiplying the results of two discrete Fourier transforms and taking the inverse results in a circular convolution. This is fairly straight-forward to prove given the definitions of the DFT and its inverse.

Given the DFT:

    \[ X \left[ k \right] = \sum_{n=0}^{N-1} e^{-\frac{j 2 \pi n k}{N}} x \left[ n \right] \]

And the inverse DFT:

    \[ x \left[ n \right] = \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n k}{N}} X \left[ k \right] \]

We can derive convolution as:

    \[ \begin{array}{@{}ll@{}} y \left[ n \right] &= \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n k}{N}} \left( \sum_{m=0}^{N-1} e^{-\frac{j 2 \pi m k}{N}} a \left[ m \right] \right) \left( \sum_{l=0}^{N-1} e^{-\frac{j 2 \pi l k}{N}} b \left[ l \right] \right) \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m \right] b \left[ l \right] \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi \left( n - l - m \right) k}{N}} \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m N + \left( n - l \right) \mathrm{ mod } N \right] b \left[ l \right] \end{array} \]

The last step in the above recognises that the summation over k is only non-zero for certain values of l, m and n and we make a variable swap of m to attain the result. We can write the above in matrix form as:

    \[ \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{pmatrix} = \begin{pmatrix} a_0 & a_{N-1} & a_{N-2} & \dots & a_1 \\ a_1 & a_0 & a_{N-1} & & \\ a_2 & a_1 & a_0 & \ddots & \\ \vdots & & \ddots & \ddots & a_{N-1} \\ a_{N-1} & & & a_1 & a_0 \end{pmatrix} \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{N-1} \end{pmatrix} \]

The matrix of a coefficients is a circulant matrix. Each row is a shifted copy of the preceeding row. Given that there exist \mathcal{O} \left( n \log n \right) algorithms for computing the DFT, we have shown that multiplying a vector by a circulant matrix has an efficient algorithm (note – this is only a computational reality for large N).

Circular Convolution with a Generalised DFT

Let’s redefine our DFT as:

    \[ X \left[ k \right] = \sum_{n=0}^{N-1} e^{-\frac{j 2 \pi n \left( k + \alpha \right) }{N}} x \left[ n \right] \]

Which has an inverse of:

    \[ x \left[ n \right] = \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n \left( k + \alpha \right)}{N}} X \left[ k \right] \]

This generalisation gives us some control over the boundary conditions of the DFT and hence the assumed data periodicity i.e. the DFT assumes the transformed data continues forever being repeated verbatim over and over – we can change this using \alpha. Let’s derive the convolution again:

    \[ \begin{array}{@{}ll@{}} y \left[ n \right] &= \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi n \left( k + \alpha \right)}{N}} \left( \sum_{m=0}^{N-1} e^{-\frac{j 2 \pi m \left( k + \alpha \right)}{N}} a \left[ m \right] \right) \left( \sum_{l=0}^{N-1} e^{-\frac{j 2 \pi l \left( k + \alpha \right)}{N}} b \left[ l \right] \right) \\ &= \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m \right] b \left[ l \right] \frac{1}{N} \sum_{k=0}^{N-1} e^{\frac{j 2 \pi \left( n - l - m \right) \left( k + \alpha \right)}{N}} \\ \end{array} \]

Given this, we can now draw up some matrices for various values of \alpha. Some interesting (and perhaps useful) values are \frac{1}{4}, \frac{3}{4} and \frac{1}{2}. We will restrict our attention to \frac{1}{2} which has a slightly different form once we obliterate the exponential:

    \[ Y[n] = \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} a \left[ m N + \left( n - l \right) \mathrm{ mod } 2 N \right] b \left[ l \right] \left( -1 \right)^m \]

We find that the matrix operation for this looks like:

    \[ \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \end{pmatrix} = \begin{pmatrix} a_0 & -a_{N-1} & -a_{N-2} & \dots & -a_1 \\ a_1 & a_0 & -a_{N-1} & & \\ a_2 & a_1 & a_0 & \ddots & \\ \vdots & & \ddots & \ddots & -a_{N-1} \\ a_{N-1} & & & a_1 & a_0 \end{pmatrix} \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{N-1} \end{pmatrix} \]

This matrix is no-longer circulant; all entries to the right of the main diagonal have been negated. This convolution might not have practical value by itself, but the symmetry suggests that it might have value when combined with another.

Toeplitz Matrix Vector Multiplication

A Toeplitz matrix has the form:

    \[ \begin{pmatrix} t_0 & t_{-1} & t_{-2} & \dots & t_{-(N-1)} \\ t_1 & t_0 & t_{-1} & & \\ t_2 & t_1 & t_0 & \ddots & \\ \vdots & & \ddots & \ddots & t_{-1} \\ t_{N-1} & & & t_1 & t_0 \end{pmatrix} \]

There are efficient algorithms for performing Toeplitz matrix by vector multiplication that use a circular convolution algorithm. These algorithms end up throwing away much of the computed result (this can be seen in the previous link in the multiplication by a zero matrix). We can avoid this by using the symmetry defined in the previous section.

If we take the sum of a regular DFT convolution of A and X and the previously defined convolution of B and X, we are effectively computing the following matrix operation:

    \[ \begin{small} \begin{pmatrix} a_0+b_0 & a_{N-1}-b_{N-1} & a_{N-2}-b_{N-2} & \dots & a_1-b_1 \\ a_1+b_1 & a_0+b_0 & a_{N-1}-b_{N-1} & & \\ a_2+b_2 & a_1+b_1 & a_0+b_0 & \ddots & \\ \vdots & & \ddots & \ddots & a_{N-1}-b_{N-1} \\ a_{N-1}+b_{N-1} & & & a_1+b_1 & a_0+b_0 \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{N-1} \end{pmatrix} \end{small} \]

We can select values for a and b to create a multiplication of any matrix with equal elements on each diagonal – such as our Toeplitz matrix – using only DFTs of the length of the input data sequence. There are some interesting optimisations that can be made when the matrix is Hermitian to eliminate some of the DFTs entirely.