# hints for implementing linear convolution

convolve :: a b -> output

convolution multiplies each sample of one array with each sample of another array and sums resulting samples into output starting at the offset of each input sample.
the second array is often called the impulse response kernel. 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.
convolution in the time domain (sample values) corresponds to multiplication in the frequency domain (spectrum). in other words, adding two sample arrays pointwise gives you a sum of frequencies but convolving them gives you a multiplication of frequencies

# links
* [applications in various fields](https://www.allaboutcircuits.com/technical-articles/dsp-applications-of-convolution-part-2/)
* [explanation on dspguide.com](http://www.dspguide.com/ch6.htm)
* [explanations on stackoverflow](https://dsp.stackexchange.com/a/43289)
* on wikipedia: [convolution](https://en.wikipedia.org/wiki/Convolution), [convolution theorem](https://en.wikipedia.org/wiki/Convolution_theorem)

# 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](https://dsp.stackexchange.com/a/43289)

# 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;
  };
};
~~~

# notes
* convolution can be made seamless for frames of continuous input data by writing input-length samples into output and the additional impulse-response-length minus one samples that would occur into an additional carryover buffer that is written into the beginning of the output of the next call. there are some additional caveats to that when the impulse response is bigger than the input array. for this case, the carryover result can be written partially into output and shifted. generally, it is possible to convolve seamlessly sample by sample
* c code for an advanced implementation can be found in [sph-sp](https://github.com/sph-mn/sph-sp)
* in numpy convolve is implemented using correlate