2020-04-09

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

(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)

(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

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

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"

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

// 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; }; };

- 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
- in numpy convolve is implemented using correlate