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