Given two finite sequences, and , of length , their -point circular convolution can be written in matrix form  as follows
where is also a sequence of length . Note that the matrix above is a circulant matrix, as each row is a circularly right-shifted version of the previous one. Each is obtained by taking the inner product of the -th row of the matrix with the 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”…
Take another look at the matrix above and note that its last row, when viewed as a list, is the reversal of list , 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 (:)) 
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 -th row is a circularly left-shifted version of the -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.
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 , , , , et cetera. We do not need an infinite list of lists.
We now know how to generate a list of lists representing the matrix.
To obtain the circular convolution of and , we multiply a matrix by a (column) vector, which consists of computing 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)
*Main> let xs = [1..4] *Main> let ys = [0..3] *Main> zipWith (*) xs ys [0,2,6,12] *Main> sum(zipWith (*) xs ys) 20
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, and . Using matrix-vector multiplication, we obtain:
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 , 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 using map. Note that
is a partial application of the inner product function, as we provide only the first argument. We could define a function as follows:
*Main> let g = innerProd [1,1,1,1] *Main> :type g g :: [Integer] -> Integer *Main> g [0..3] 6
Note that 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.
Last but not least, let us carry out a couple of tests using examples from Mitra’s book :
*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.
 Sanjit K. Mitra, Digital Signal Processing: a computer-based approach, 3rd edition, McGraw-Hill, 2005.