## Julia Community 🟣

Steven Siew

Posted on • Updated on

# Convolution in the form of a binary infix operator

## Basic introduction to Convolution

Here is a quick mathematic problem:

What is the result of the following polynomial multiplication?

(1 + 2*X + 3*X^2) * (4 + 5*X + 6*X^2 + 7*X^3 + 8*X^4 +9*X^5)

well the answer is (using mathematica) 4 + 13*X + 28*X^2 + 34*X^3 + 40*X^4 + 46*X^5 + 42*X^6 + 27*X^7

You can obtain the answer by using the convolution function in Julia. Unfortunately, the convolution function is NOT in base Julia but is found in the DSP.jl package

So after you have added DSP to your manifest using package manager, you can do the following.

``````julia> import DSP: conv

julia> conv([1,2,3],[4,5,6,7,8,9])
8-element Vector{Int64}:
4
13
28
34
40
46
42
27
``````

You can perform the convolution as a binary infix function if you define the operator \odot as the convolution function.

``````# help?>  ⊙
# "⊙" can be typed by \odot<tab>
``````

Now you can do this

``````julia> ⊙(x,y) = conv(x,y)
⊙ (generic function with 1 method)

julia> [1,2,3] ⊙ [4,5,6,7,8,9]
8-element Vector{Int64}:
4
13
28
34
40
46
42
27
``````

### 3blue1brown explaination of what convolution is

Here is the link to the youtube video about convolution by 3blue1brown. He explains what a convolution really is.

### How to use convolution to predict the outcome of an election for a political party

Suppose the Lion political party has three candidates (Adam,Bob and Charles) in three different electorates in an election such that

• Adam has 80% chance of winning his seat
• Bob has 70% chance of winning his seat
• Charles has 60% chance of winning his seat

What are the chances of the Lion political party winning 0,1,2,3 seats?

``````julia> adam = [(1.0 - 0.8),0.8]
2-element Vector{Float64}:
0.19999999999999996
0.8

julia> bob = [(1.0 - 0.7),0.7]
2-element Vector{Float64}:
0.30000000000000004
0.7

julia> charles = [(1.0-0.6),0.6]
2-element Vector{Float64}:
0.4
0.6

julia> adam ⊙ bob ⊙ charles
4-element Vector{Float64}:
0.02400000000000002
0.18800000000000006
0.45199999999999996
0.33599999999999997
``````

• 2.4 % chance of winning no seats
• 18.8 % chance of winning only one seats
• 45.2 % chance of winning only two seats
• 33.6 % chance of winning all three seats

### Using convolution to do ordinary multiplication

``````julia> 1234 * 5678
7006652

'''
(4 + 3 * 10^1 + 2 * 10^2 + 1 * 10^3) * (8 + 7 * 10^1 + 6 * 10^2 + 5 * 10^3)
'''

julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
32
52
61
60
34
16
5

julia> for k = 2:length(result)
result[k] = result[k] + (result[k-1] ÷ 10)
result[k-1] = result[k-1] % 10
end

julia> result
7-element Vector{Int64}:
2
5
6
6
0
0
7

julia> result = reverse(result)
7-element Vector{Int64}:
7
0
0
6
6
5
2

julia> string.(result) |> join
"7006652"
``````

You may think that it is a bit silly to use convolution to do multiplication but if you do not have access to BigInt then your only choice left is to use convolution to multiple two gigantic integers together.

### Implementing convolution using Fast Fourier Transform

We can implement convolution using fast fourier transform as this article below shows

``````julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
32
52
61
60
34
16
5

# Now we do the same thing using FFT

julia> using FFTW

julia> a=Float64.([4,3,2,1])
4-element Vector{Float64}:
4.0
3.0
2.0
1.0

julia> b=Float64.([8,7,6,5])
4-element Vector{Float64}:
8.0
7.0
6.0
5.0

julia> len_final = length(a) + length(b) - 1
7

julia> extra_padding_for_a = len_final - length(a)
3

julia> aa = vcat(a,[0,0,0])
7-element Vector{Float64}:
4.0
3.0
2.0
1.0
0.0
0.0
0.0

julia> extra_padding_for_b = len_final - length(b)
3

julia> bb = vcat(b,[0,0,0])
7-element Vector{Float64}:
8.0
7.0
6.0
5.0
0.0
0.0
0.0

julia> a_fft = fft(aa)
7-element Vector{ComplexF64}:
10.0 + 0.0im
4.524458669761152 - 4.729234010885294im
2.1539892641849523 - 1.275184775842325im
2.3215520660538953 - 0.7129161645984384im
2.3215520660538953 + 0.7129161645984384im
2.1539892641849523 + 1.275184775842325im
4.524458669761152 + 4.729234010885294im

julia> b_fft = fft(bb)
7-element Vector{ComplexF64}:
26.0 + 0.0im
6.524458669761152 - 13.491806545954942im
4.153989264184952 - 0.31203553822726793im
4.321552066053895 - 3.2208368399238467im
4.321552066053895 + 3.2208368399238467im
4.153989264184952 + 0.31203553822726793im
6.524458669761152 + 13.491806545954942im

julia> product_fft = a_fft .* b_fft
7-element Vector{ComplexF64}:
260.0 + 0.0im
-34.2862667915158 - 91.89881294123597im
8.54974531072476 - 5.969225068086821im
7.736521480791037 - 10.558244744191306im
7.736521480791037 + 10.558244744191306im
8.54974531072476 + 5.969225068086821im
-34.2862667915158 + 91.89881294123597im

julia> ff = ifft(product_fft)
7-element Vector{ComplexF64}:
32.0 + 0.0im
52.0 + 0.0im
61.0 + 0.0im
60.00000000000001 + 0.0im
34.00000000000001 + 0.0im
16.0 + 0.0im
5.0 + 0.0im

julia> final_result = real.(ff)
7-element Vector{Float64}:
32.0
52.0
61.0
60.00000000000001
34.00000000000001
16.0
5.0

# compare with
julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
32
52
61
60
34
16
5
``````

### Now put them together to perform integer multiplication using FFT

``````function fft_multiplication(str_a,str_b)
local arr_a, arr_b, len_final
arr_a = Float64.( reverse(map(x->parse(Int64,x),split(str_a,""))) )
arr_b = Float64.( reverse(map(x->parse(Int64,x),split(str_b,""))) )
len_final = length(arr_a) + length(arr_b) - 1
arr_a = vcat(arr_a,[ 0.0 for k in 1:extra_padding_for_a ])
arr_b = vcat(arr_b,[ 0.0 for k in 1:extra_padding_for_b ])
local a_fft = fft(arr_a)
local b_fft = fft(arr_b)
local product_fft = a_fft .* b_fft
local complex_temp = ifft(product_fft)
local result = real.(complex_temp)
result = Int64.(  map(x->round(x,RoundNearest),result)  )
local k
for k = 2:length(result)
result[k] = result[k] + (result[k-1] ÷ 10)
result[k-1] = result[k-1] % 10
end
result = reverse(result)
result = join(string.(result))
return result
end

julia> println( fft_multiplication("12345","67890")  )
838102050

# compare with
julia> 12345 * 67890
838102050
``````

### Finally the algorithm for a very naive implementation of convolution

How can we talk about convolution without a naive implementation of convolution. This is a very inefficient implementation of an algorithm to calculate convolution.

``````function naive_convolution(data, kernel)
local len_data = length(data)
local len_kernel = length(kernel)
local len_ans = len_data + len_kernel - 1
local typeofdata = eltype(data)
local ans = zeros(typeofdata, len_ans)
local j,k,acc,loc

# Next we work out the value for each of the element of the solution
for j = 1:len_ans
acc = zero(typeofdata)
for k = 1:len_kernel
loc = j - (k - 1)
if !( loc < 1 || loc > len_data )
acc += data[loc] * kernel[k]
end
end
ans[j] = acc
end
return ans
end
``````

### Why is the output of FFT symmetrical?

Because there are two identical frequencies. One is positive and travel forward in time and the other is negative and travel backwards in time.

Okay, okay. Please do not take me seriously. I am just joking.