2017-11-12

hints for implementing linear convolution

convolve :: a b -> output

links

applications in various fields

explanation on dspguide.com

explanations on stackoverflow

wikipedia description

input/output examples

(1 2 3) (1) -> (1 2 3)
(1 2 3) (2) -> (2 4 6)
(2 0 0) (1) -> (2 0 0)
(2 0 0) (1 0 0) -> (2 0 0 0 0)
(2 0 0) (0 1 0) -> (0 2 0 0 0)
(2 2 0) (0 1 0) -> (0 2 2 0 0)
(2 0 0) (1 1) -> (2 2 0 0)
(2 2 0) (1 1) -> (2 4 2 0)
(3 2) (2 2) -> (6 10 4)
(3 2) (1 1 1 1 1 1 1) -> (3 5 5 5 5 5 5 2)
(1 1 1 1) (1 1 1 1) -> (1 2 3 4 3 2 1)
(1 1 1 1) (2 2 2 2) -> (2 4 6 8 6 4 2)
(2 2 2 2) (1 1 1 1) -> (2 4 6 8 6 4 2)
(1 1 1 1) (3 3 3 3) -> (3 6 9 12 9 6 3)
(1 1 1 1) (0 1 0) -> (0 1 1 1 1 0)
(1 1 1 1) (1 0 1) -> (1 1 2 2 1 1)
(1 1 1 1) (1 1 1) -> (1 2 3 3 2 1)

implementation in scheme

(define (produce-one f a b) (map (lambda (b) (f a b)) b))

(define (convolve a b)
  (let loop ((products (produce-one * (first a) b)) (a (tail a)))
    (if (null? a) products
      (pair (first products)
        (loop (map + (append (tail products) (list 0)) (produce-one * (first a) b)) (tail a))))))

algorithm: sum the tails of productions for each element of "a"

signature: list list -> list

implementation in coffeescript, javascript

convolve = (a, b, result) ->
  a_index = 0
  while a_index < a.length
    b_index = 0
    while b_index < b.length
      result_index = a_index + b_index
      result[result_index] += a[a_index] * b[b_index]
      b_index += 1
    a_index += 1
  result

example call

convolve [1, 2, 3], [4, 5, 6], [0, 0, 0, 0, 0] -> [4, 13, 28, 27, 18]

the inner loop will be repeated with the following values for result_index, a_index and b_index

0 0 0
1 0 1
2 0 2
1 1 0
2 1 1
3 1 2
2 2 0
3 2 1
4 2 2

result_index is a_index + b_index

a_index advances for each production of one element from "a" with all elements from "b"

alternative implementation

convolve = (a, b, result) ->
  result_index = 0
  while result_index < result.length
    state = 0
    b_length = b.length
    b_index = 0
    while b_index < b_length
      a_index = result_index + b_index - b_length + 1
      if a_index >= 0 and a_index < b_length
        state += a[a_index] * b[b_length - b_index - 1]
      b_index += 1
    result[result_index] = state
    result_index += 1
  result

adapted from source

implementation in c

/** discrete linear convolution.
    result length must be at least a-len + b-len - 1 */
void convolve(double *result, double *a, size_t a_len,
              double *b, size_t b_len) {
  size_t a_index = 0;
  size_t b_index = 0;
  while (a_index < a_len) {
    while (b_index < b_len) {
      result[a_index + b_index] += a[a_index] * b[b_index];
      b_index += 1;
    };
    b_index = 0;
    a_index += 1;
  };
};

other

an impulse response allows each input sample to influence the output for impulse-response-length samples. an impulse response of length one maps sample to sample

in numpy convolve is implemented using correlate


tags: programming start q1 document guide computer signal-processing convolution