Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How do I do convolution in F#?

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) 
like image 314
Hallgrim Avatar asked Dec 10 '22 21:12

Hallgrim


1 Answers

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)
like image 96
Nathan Shively-Sanders Avatar answered Dec 25 '22 13:12

Nathan Shively-Sanders