Given two finite sequences, $x$ and $y$, of length $n$, their $n$-point circular convolution can be written in matrix form [1] as follows

$\left[\begin{array}{c} z_0\\ z_1\\ z_2\\ \vdots\\ z_{n-2}\\ z_{n-1}\end{array}\right] = \left[\begin{array}{cccccc} y_0 & y_{n-1} & y_{n-2} & \ldots & y_2 & y_1\\ y_1 & y_0 & y_{n-1} & \ldots & y_3 & y_2\\ y_2 & y_1 & y_0 & \ldots & y_4 & y_3\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ y_{n-2} & y_{n-3} & y_{n-4} & \ldots & y_0 & y_{n-1}\\ y_{n-1} & y_{n-2} & y_{n-3} & \ldots & y_1 & y_0\end{array}\right] \left[\begin{array}{c} x_0\\ x_1\\ x_2\\ \vdots\\ x_{n-2}\\ x_{n-1}\end{array}\right]$

where $z$ is also a sequence of length $n$. Note that the $n \times n$ matrix above is a circulant matrix, as each row is a circularly right-shifted version of the previous one. Each $z_k$ is obtained by taking the inner product of the $k$-th row of the matrix with the $x$ vector.

In Digital Signal Processing (DSP), we often think of finite sequences as vectors, which allows us to compute convolutions using matrix operations. Taking a more Computer Science-ish approach, one can think of finite sequences as finite lists, and compute convolutions using list operations. The matrix above can be thought of as a list of lists (where each row is a list). We know how to multiply a matrix by a vector, but how do we “multiply” a list of lists by a list?

I started teaching myself Haskell this Summer, and I am (obviously) looking for interesting problems at which to throw some functional firepower. Solutions looking for a problem. But, first, let us “begin at the beginning”…

__________

Reverse

Take another look at the matrix above and note that its last row, when viewed as a list, is the reversal of list $y$, which means that we need a function that reverses lists. Fortunately, we do not need to implement such a function, as there is the function reverse in Prelude that uses foldl and flip:

reverse :: [a] -> [a]
reverse = foldl (flip (:)) []

__________

Circular shifts

Yet once again, let us look at the matrix at the top of this post. Note that the 2nd row is a circularly right-shifted version of the 1st one, and that the 3rd row is a circularly right-shifted version of the 2nd one, et cetera. Also, note that the $(n-1)$-th row is a circularly left-shifted version of the $n$-th row, and so and so on. Thus, we need to implement functions that perform circular shifts.

The following function circularly right-shifts a list:

circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

This function merely takes the last element of the list and moves it to the beginning. We can circularly left-shift a list using the function:

circShiftL :: [a] -> [a]
circShiftL [] = []
circShiftL xs = (tail xs) ++ [head xs]

which moves the first element of a list to the end. Let us now test these functions on some simple lists using the GHCi interpreter:

*Main> circShiftR [0..9]
[9,0,1,2,3,4,5,6,7,8]
*Main> circShiftR (circShiftR [0..9])
[8,9,0,1,2,3,4,5,6,7]
*Main> (circShiftR . circShiftR) [0..9]
[8,9,0,1,2,3,4,5,6,7]
*Main> circShiftL [0..9]
[1,2,3,4,5,6,7,8,9,0]
*Main> (circShiftL . circShiftL) [0..9]
[2,3,4,5,6,7,8,9,0,1]

Note that in the 3rd and 5th inputs we composed the circular shifting functions to perform double right- and left-shifts.

__________

Iterate

Since each row of the matrix is a right-shifted version of the previous one or a left-shifted version of the next one, we can generate the matrix from its first or last rows. Putting it in terms of lists, we can generate a list of lists that represents the matrix from a single list using the iterate function in Prelude:

iterate :: (a -> a) -> a -> [a]
iterate f x =  x : iterate f (f x)

and the circular shifting functions we implemented before. To make things concrete, we can generate all four circular right-shifts of a list of length 4, as follows:

*Main> take 4 (iterate circShiftR [0..3])
[[0,1,2,3],[3,0,1,2],[2,3,0,1],[1,2,3,0]]

A word of caution is in order. Note that we take the first four lists, as

iterate f x

does generate an infinite list containing $x$, $f (x)$, $(f \circ f ) (x)$, $(f \circ f \circ f ) (x)$, et cetera. We do not need an infinite list of lists.

We now know how to generate a list of lists representing the matrix.

__________

Inner product

To obtain the circular convolution of $x$ and $y$, we multiply a matrix by a (column) vector, which consists of computing $n$ inner products. However, in this post we are thinking in terms of lists. We thus need to implement a function that computes the inner product of two (finite) lists, as follows:

innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

where we first use zipWith to obtain a list of the $x_i y_i$ products, and then use sum to compute the sum $\sum_i x_i y_i$. Let us test the inner product function:

*Main> let xs = [1..4]
*Main> let ys = [0..3]
*Main> zipWith (*) xs ys
[0,2,6,12]
*Main> sum(zipWith (*) xs ys)
20

__________

Circular convolution

We finally have all we need to compute the circular convolution.

For example, suppose we want to compute the circular convolution of two 4-point sequences, say, $x = (1,1,1,1)$ and $y = (0,1,2,3)$. Using matrix-vector multiplication, we obtain:

$\left[\begin{array}{cccc} 0 & 3 & 2 & 1\\ 1 & 0 & 3 & 2\\ 2 & 1 & 0 & 3\\ 3 & 2 & 1 & 0\end{array}\right] \left[\begin{array}{c} 1\\ 1\\ 1\\ 1\end{array}\right] = \left[\begin{array}{c} 6\\ 6\\ 6\\ 6\end{array}\right]$

Alternatively, we can compute the circular convolution by representing sequences by lists and using the functions we have implemented before:

*Main> -- define lists
*Main> let xs = [1,1,1,1]
*Main> let ys = [0,1,2,3]
*Main> -- reverse and right-shift ys
*Main> let ys' = (circShiftR . reverse) ys
*Main> ys'
[0,3,2,1]
*Main> -- compute list of lists representing the matrix
*Main> let yss = take 4 (iterate circShiftR ys')
*Main> yss
[[0,3,2,1],[1,0,3,2],[2,1,0,3],[3,2,1,0]]
*Main> -- compute inner product of xs with each list in yss
*Main> map (innerProd xs) yss
[6,6,6,6]

Note that I wrote comments on the GHCi command line in order to improve readability. We are now ready to implement the function that computes the circular convolution:

circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
where
n   = length xs
ys' = (circShiftR . reverse) ys
yss = take n (iterate circShiftR ys')

How does it work? We first take sequence $y$, reverse and right-shift it (to obtain the 1st row of the matrix), generate the remaining rows by right-shifting the previous ones, and finally compute the inner product of each row with sequence $x$ using map. Note that

innerProd xs

is a partial application of the inner product function, as we provide only the first argument. We could define a function $g$ as follows:

*Main> let g = innerProd [1,1,1,1]
*Main> :type g
g :: [Integer] -> Integer
*Main> g [0..3]
6

Note that $g$ takes a list of integers and returns its inner product with the list of four ones. Putting it all together in a .hs file, we have:

-- circular right-shift
circShiftR :: [a] -> [a]
circShiftR [] = []
circShiftR x = last x : init x

-- inner product of two lists
innerProd :: Num a => [a] -> [a] -> a
innerProd xs ys = sum(zipWith (*) xs ys)

-- circular convolution of two lists
circConv :: Num a => [a] -> [a] -> [a]
circConv xs ys = map (innerProd xs) yss
where
n   = length xs
ys' = (circShiftR . reverse) ys
yss = take n (iterate circShiftR ys')

We use several higher-order functions in this implementation: map, iterate, zipWith. We also use function composition and partial function application. A lot of beautiful ideas in only a dozen lines of code!

I am a Haskell neophyte, and I make no claims that the code above cannot be improved. If you have any constructive suggestions, I would be happy to hear them. Please use the HTML tags <pre> and </pre> to post code in the comments.

__________

Testing

Last but not least, let us carry out a couple of tests using examples from Mitra’s book [1]:

*Main> -- problem 5.2 b)
*Main> let xs = [-1,5,3,0,3]
*Main> let hs = [-2,0,5,3,-2]
*Main> circConv xs hs
[1,-1,-2,16,26]
*Main> -- problem 5.45 b)
*Main> let gs = [2,-1,3,0]
*Main> let hs = [-2,4,2,-1]
*Main> circConv gs hs
[3,7,-6,8]

It appears to be working! Please let me know if you find any bugs.

__________

References

[1] Sanjit K. Mitra, Digital Signal Processing: a computer-based approach, 3rd edition, McGraw-Hill, 2005.