How do I do convolution in F#?

Question

I would like convolve a discrete signal with a discrete filter. The signal and filter is sequences of float in F#.

The only way I can figure out how to do it is with two nested for loops and a mutable array to store the result, but it does not feel very functional.

Here is how I would do it non-functional:

conv = double[len(signal) + len(filter) - 1]
for i = 1 to len(signal)
  for j = 1 to len(filter)
    conv[i + j] = conv[i + j] + signal(i) * filter(len(filter) - j) 

Answer 1

I don't know F#, but I'll post some Haskell and hopefully it will be close enough to use. (I only have VS 2005 and an ancient version of F#, so I think it would be more confusing to post something that works on my machine)

Let me start by posting a Python implementation of your pseudocode to make sure I'm getting the right answer:

def convolve(signal, filter):
    conv = [0 for _ in range(len(signal) + len(filter) - 1)]
    for i in range(len(signal)):
        for j in range(len(filter)):
            conv[i + j] += signal[i] * filter[-j-1]
    return conv

Now convolve([1,1,1], [1,2,3]) gives [3, 5, 6, 3, 1]. If this is wrong, please tell me.

The first thing we can do is turn the inner loop into a zipWith; we're essentially adding a series of rows in a special way, in the example above: [[3,2,1], [3,2,1], [3,2,1]]. To generate each row, we'll zip each i in the signal with the reversed filter:

makeRow filter i = zipWith (*) (repeat i) (reverse filter)

(Note: according to a quick google, zipWith is map2 in F#. You might have to use a list comprehension instead of repeat) Now:

makeRow [1,2,3] 1
=> [3,2,1]
makeRow [1,2,3] 2
=> [6,4,2]

To get this for all i, we need to map over signal:

map (makeRow filter) signal
=> [[3,2,1], [3,2,1], [3,2,1]]

Good. Now we just need a way to combine the rows properly. We can do this by noticing that combining is adding the new row to the existing array, except for the first element, which is stuck on front. For example:

[[3,2,1], [6,4,2]] = 3 : [2 + 6, 1 + 4] ++ [2]
// or in F#
[[3; 2; 1]; [6; 4; 2]] = 3 :: [2 + 6; 1 + 4] @ [2]

So we just need to write some code that does this in the general case:

combine (front:combinable) rest =
    let (combinable',previous) = splitAt (length combinable) rest in
    front : zipWith (+) combinable combinable' ++ previous

Now that we have a way to generate all the rows and a way to combine a new row with an existing array, all we have to do is stick the two together with a fold:

convolve signal filter = foldr1 combine (map (makeRow filter) signal)

convolve [1,1,1] [1,2,3]
=> [3,5,6,3,1]

So that's a functional version. I think it's reasonably clear, as long as you understand foldr and zipWith. But it's at least as long as the imperative version and like other commenters said, probably less efficient in F#. Here's the whole thing in one place.

makeRow filter i = zipWith (*) (repeat i) (reverse filter)
combine (front:combinable) rest =
    front : zipWith (+) combinable combinable' ++ previous
    where (combinable',previous) = splitAt (length combinable) rest
convolve signal filter = foldr1 combine (map (makeRow filter) signal)

Edit:

As promised, here is an F# version. This was written using a seriously ancient version (1.9.2.9) on VS2005, so be careful. Also I couldn't find splitAt in the standard library, but then I don't know F# that well.

open List
let gen value = map (fun _ -> value)
let splitAt n l = 
  let rec splitter n l acc =
    match n,l with
    | 0,_ -> rev acc,l
    | _,[] -> rev acc,[]
    | n,x::xs -> splitter (n - 1) xs (x :: acc)
  splitter n l [] 
let makeRow filter i = map2 ( * ) (gen i filter) (rev filter)
let combine (front::combinable) rest =
  let combinable',previous = splitAt (length combinable) rest
  front :: map2 (+) combinable combinable' @ previous
let convolve signal filter = 
  fold1_right combine (map (makeRow filter) signal)

Answer 2

Indeed, you generally want to avoid loops (plain, nested, whatever) and anything mutable in functional programming.

There happens to be a very simple solution in F# (and probably almost every other functional language):

let convolution = Seq.zip seq1 seq2

The zip function simply combines the two sequences into one of pairs containing the element from seq1 and the element from seq2. As a note, there also exist similar zip functions for the List and Array modules, as well as variants for combining three lists into triples (zip3). If you want tom ore generally zip (or "convolute") n lists into a list of n-tuples, then you'll need to write your own function, but it's pretty straightforward.

(I've been going by this description of convolution by the way - tell me if you mean something else.)

Answer 3

In principle, it should be possible to use the (Fast) Fourier Transform, or the related (Discrete) Cosine Transform, to calculate the convolution of two functions reasonably efficiently. You calculate the FFT for both functions, multiply them, and apply the inverse FFT on the result.

mathematical background

That's the theory. In practice you'd probably best find a math library that implements it for you.

Answer 4

Try this function:

let convolute signal filter =
    [|0 .. Array.length signal + Array.length filter - 1|] |> Array.map (fun i ->
        [|0 .. i|] |> Array.sum_by (fun j -> signal.[i] * filter.[Array.length filter - (i - j) - 1]))

It's probably not the nicest function solution, but it should do the job. I doubt there exists a purely functional solution that will match the imperative one for speed however.

Hope that helps.

Note: The function is currently untested (though I've confirmed it compiles). Let me know if it doesn't quite do what it should. Also, observe that the i and j variables do not refer to the same things as is your original post.