Skip to content

Instantly share code, notes, and snippets.

@vthriller
Forked from fredo-dedup/wifi-field
Last active December 24, 2020 06:10
Show Gist options
  • Save vthriller/db8560f5b358d86a67a6f32b570acc4c to your computer and use it in GitHub Desktop.
Save vthriller/db8560f5b358d86a67a6f32b570acc4c to your computer and use it in GitHub Desktop.
Calculating WIFI propagation with IJulia
### A Pluto.jl notebook ###
# v0.12.17
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
el
end
end
# ╔═╡ 1cf76d74-450a-11eb-3900-5dbbe000047a
begin
using Colors
using ColorTypes
using Images
using SparseArrays
using DSP: conv2
using Printf
using PlutoUI
end
# ╔═╡ 14033bfc-4506-11eb-345d-b968dbdd669f
md"# Where's my Wi-Fi signal?"
# ╔═╡ b2375ede-4509-11eb-02c1-fde861ed36ad
md"## Solving the Helmhotz equation to find the good spots in your home"
# ╔═╡ b76f6990-4509-11eb-1f2c-8933290c6292
md"[*Helmhurts* by Jason Cole](http://jasmcole.com/2014/08/25/helmhurts/) explains how to apply the Maxwell equations, and in particular their smaller sibling the Helmholtz equations, to calculate the propagation of Wi-Fi signals."
# ╔═╡ bc937a5e-4509-11eb-0cfa-c5634086f187
md"What follows is an attempt to use Julia to go through the whole process described in the blog entry."
# ╔═╡ 1718edba-450a-11eb-2377-113207096427
md"### Setting the up the data"
# ╔═╡ 18908fcc-450a-11eb-18fb-53d6a09bd211
md"We'll start first by importing a floor plan, in that case a bitmap image that I produced with an image editor, but you can of course create your own image of your own place (or an imaginary one)."
# ╔═╡ 22808b72-450a-11eb-25ae-15d190b9c876
img = load("plan.png")
# ╔═╡ 3432951a-450a-11eb-0a87-05625cdca2f4
size(img)
# ╔═╡ 87240306-450f-11eb-2b4b-93419692a783
typeof(img)
# ╔═╡ 3531c2c2-450a-11eb-3e62-e31c5ad04546
md"Here we have a 2D-array of gray pixels.
For the purpose of the simulation we only need to know if a given point is air or concrete. Let's extract that material information from the pixels:"
# ╔═╡ 3d210916-450a-11eb-0ab4-4bde681a8fa2
plan = gray.(img);
# ╔═╡ 434cbaa6-450a-11eb-1e22-d50acd4eaabc
md"And set the necessary variables:"
# ╔═╡ 476333fe-450a-11eb-1b11-5bf8767c8847
begin
dimx, dimy = size(plan) # spatial dimensions
fsize = length(plan) # full size of the problem
δ = 0.03 # spatial resolution of each pixel
end;
# ╔═╡ 4efa869e-450a-11eb-222f-79c99c8ee078
md"Note that the spatial resolution is given by the actual size of the appartment divided by the number of pixel making the plan picture. Two important considerations come into play here:
1. Beware of not making the problem size (the number of pixels) too important, I found on my setup that you can go quickly from a problem that can be solved in seconds (100 x 100) to a never ending calculation or one that throws out of memory exceptions (1000 x 1000 in my case).
2. You can't on the other hand keep the pixel count low by increasing the spatial step δ (i.e. decreasing the resolution): when it gets too close to the wavelength (about 10 cm for Wi-Fi signals) the mesh becomes too coarse to give meaningful results.
Here, with a spatial resolution of 3 cm we should be safe for the resolution while still simulating a normally sized appartment of about $(@sprintf(\"%.2fm x %.2fm\", dimx * δ, dimy * δ))"
# ╔═╡ 55198d04-450a-11eb-0529-8f14d12e67cf
begin
η_air = 1. # refraction index for air
# refraction index for concrete
# the imaginary part conveys the absorption
# https://www.itu.int/dms_pubrec/itu-r/rec/p/R-REC-P.2040-1-201507-I!!PDF-E.pdf
# §3 Compilations of electrical properties of materials
η_concrete = begin
local f = 2.4 # GHz
local a, b, c, d = 5.31, 0, 0.0326, 0.8095
local η_re = a * (b != 0 ? f ^ b : 1)
local σ = c * (d != 0 ? f ^ d : 1)
local η_im = 17.98im * σ / f
η_re + η_im
end
λ = 0.12 # for a 2.5 GHz signal, wavelength is ~ 12cm
k = 2π / λ # k is the wavenumber
end;
# ╔═╡ 6af3e020-450a-11eb-09ea-392a4fe9d03e
md"We'll now create the Complex Array containing the material information:"
# ╔═╡ 73c7354e-450a-11eb-121d-0b694833a20b
begin
μ = similar(plan, Complex)
μ[plan .!= 0] .= (k / η_air)^2
μ[plan .== 0] .= (k / η_concrete)^2
end;
# ╔═╡ 77c4f118-450a-11eb-336c-1d76291a980c
md"Note that we directly calculate the $\left( \frac{k}{n} \right) ^2$ value appearing in the Helmholtz equation below."
# ╔═╡ 7c03f8aa-450a-11eb-2006-83dfdb3c3527
md"### Solving the equation"
# ╔═╡ 85ae4306-450a-11eb-02dd-7d43dc23b9e0
md"The Helmholtz equation for the amplitude of the electric field is: $$\Delta A + \left( \frac{k}{n} \right) ^2 A= f$$, where $\Delta$ is the Laplacian operator, $A$ the amplitude field and $f$ the wave source that will be equal to zero everywhere except where the antenna is."
# ╔═╡ 8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
md"The first step is to rephrase that problem in one of the ways (there are several) a computer can tackle it: by discretizing the space and using finite differences.
In that framework the Laplacian can be approximated by the neighboring values, giving us:
$$\frac{A(x+δ,y) - 2 A(x,y) + A(x-δ,y)}{δ^2}+\frac{A(x,y+δ) - 2 A(x,y) + A(x,y-δ)}{δ^2}+ \left( \frac{k}{n} \right) ^2 A(x,y) = f(x,y)$$"
# ╔═╡ ab2a9cc4-450a-11eb-23bf-1f194db3142e
md"We have here a linear combination of elements of A equal to elements of f. In others words an ordinary system of linear equations that computers are good at solving.
There is however a roadblock: that system has as many equations as there are values in A: $dimx x $dimy = $fsize. The matrix describing the problem would occupy an huge amount of memory if stored as a dense array: $fsize x $fsize x 8 x 2 (don't forget these are complex numbers) = $(@sprintf(\"%.1f\", fsize ^ 2 * 16 / 2 ^ 30)) GiB!
This memory would be mostly wasted though as the matrix is almost empty (only 5 elements appear in each equation not the whole $fsize). Sparses matrices are the adapted structure in these situations as they only store the non zero elements in memory. Julia provides a sparse matrix type on which most dense matrix methods can apply."
# ╔═╡ b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
md"Let's now build the three vectors describing the non zero elements of the sparse matrix `S` (the `LinearIndices` is used to relate the (x,y) positions in the real space to an index in the matrix):"
# ╔═╡ cbf3f270-450a-11eb-0eba-0be7069aafc4
begin
xs = Array{Int}(undef, 5*dimx*dimy)
ys = Array{Int}(undef, 5*dimx*dimy)
vs = Array{Complex}(undef, 5*dimx*dimy)
local i = 1
local ind = LinearIndices(plan)
for x in 1:dimx, y in 1:dimy # x=1; y=1
xm = (x+dimx-2) % dimx + 1
xp = x % dimx + 1
ym = (y+dimy-2) % dimy + 1
yp = y % dimy + 1
xs[i] = ind[x, y]; ys[i] = ind[ x, y]; vs[i] = μ[x,y] - 2*δ^-2
i += 1
xs[i] = ind[x, y]; ys[i] = ind[xp, y]; vs[i] = δ^-2
i += 1
xs[i] = ind[x, y]; ys[i] = ind[xm, y]; vs[i] = δ^-2
i += 1
xs[i] = ind[x, y]; ys[i] = ind[ x, yp]; vs[i] = δ^-2
i += 1
xs[i] = ind[x, y]; ys[i] = ind[ x, ym]; vs[i] = δ^-2
i += 1
end
end
# ╔═╡ ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
md"We can now declare `S` by calling the `sparse()` function:"
# ╔═╡ d466118e-450a-11eb-0260-e931d1fe5edb
S = sparse(xs, ys, vs, fsize, fsize);
# ╔═╡ d76ed16a-450a-11eb-1ff6-13a784896b1f
md"There is only the emiter field $f$ left to define:"
# ╔═╡ 8a9acc22-4582-11eb-10d1-af27fb9fa40e
# TODO re-render this cell with each update of `img`
@bind antenna_pos html"""
<span class="hack-input">
<!--
the first element is the one that's being watched,
even it it's a text node,
hence this span
-->
Click to pick antenna position: <br>
<canvas width="200" height="200" style="position: relative"></canvas>
</span>
<script>
const emitter = currentScript.closest('pluto-output').querySelector(".hack-input")
const canvas = emitter.querySelector("canvas")
const ctx = canvas.getContext("2d")
canvas.clear = () => {
// XXX UGLY HACK
// I was not able to figure out how to interpolate `img` into HTML literals,
// and HTML in Markdown strings is not supported,
// so instead I clone already existing image from the cell that defined the `img`
let img = document.getElementById("22808b72-450a-11eb-25ae-15d190b9c876").querySelector('pluto-output img')
canvas.width = img.width
canvas.height = img.height
ctx.drawImage(img, 0, 0)
}
canvas.onclick = e => {
// send the value to Pluto
// N.B. coords are getting swapped:
// JS: x→ y↓
// Julia: x↓ y→
emitter.value = [e.layerY, e.layerX]
emitter.dispatchEvent(new CustomEvent("input"))
canvas.clear();
// draw the crosshair
ctx.strokeStyle = '#c00'
ctx.beginPath()
ctx.moveTo(e.layerX - 3, e.layerY)
ctx.lineTo(e.layerX + 3, e.layerY)
ctx.stroke()
ctx.beginPath()
ctx.moveTo(e.layerX, e.layerY - 3)
ctx.lineTo(e.layerX, e.layerY + 3)
ctx.stroke()
}
// show the image for the first time
canvas.clear()
</script>
"""
# ╔═╡ a6105b1a-45a1-11eb-37bc-25e709d2aa61
antenna_pos
# ╔═╡ dd620024-450a-11eb-0d26-39444e86dbc7
begin
f = fill(0 + 0im, (dimx, dimy))
f[antenna_pos...] = 1.0 # our Wi-Fi emitter antenna will be here
end;
# ╔═╡ e13d4758-450a-11eb-21fa-21ef2b053b3c
md"Now, with the above notations, we have the linear system $S.A=f$ to solve for A:"
# ╔═╡ e53559e8-450a-11eb-10df-2f9bd44c829b
A = reshape(S \ vec(f), dimx, dimy);
# ╔═╡ ea215d76-450a-11eb-34c4-8f840a4c07e0
md"A contains the amplitude field we are looking for in the real part of the complex number. Since the relevant physical value is rather the Energy received not the amplitude, we'll calculate the square of the amplitude."
# ╔═╡ ed05c25c-450a-11eb-3ca4-fd98065261f4
E_raw = real(A) .* real(A);
# ╔═╡ f13d360c-450a-11eb-3ea5-a3c6738d3f7b
md"And since the receiving antennas of laptops and smartphones have a size of at least 10cm, we'll convolve the energy field to have a more accurate representation of the received signal strength"
# ╔═╡ 8faa5ac4-4574-11eb-113f-49565b459121
@bind E_convolve html"<input type='checkbox' checked> convolve energy field"
# ╔═╡ f4454b66-450a-11eb-1429-f3e5cd06ae3b
E = E_convolve ? conv2(E_raw, ones(5,5)/25 )[3:end-2, 3:end-2] : E_raw;
# ╔═╡ f7acc81a-450a-11eb-3c0c-61badd6b8835
md"### Plotting the result"
# ╔═╡ faf621ce-450a-11eb-0158-398569973267
md"We can now convert the E field back to an image to visualize where the good spots are for Wi-Fi reception in our appartment.
First we'll translate the energy field to 1-100 interger scale:"
# ╔═╡ ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
Ei = min.(100, max.(1, Int.(trunc.( 1 .+ 100 .* (E .- minimum(E)) / (maximum(E) - minimum(E)) ) )));
# ╔═╡ 039ef896-450b-11eb-231d-e1a15b86df0c
md"Then we'll create a color scale using the `colormap()` function of the `Colors` package:"
# ╔═╡ 07031f76-450b-11eb-1824-8367aa6aca36
cm = reverse(colormap("oranges"));
# ╔═╡ 14232eda-450b-11eb-23a1-9731174997ad
md"The last step is to create a color image by picking the RGB components of the color scale:"
# ╔═╡ 02251d32-4575-11eb-31de-cb4d694dd5b3
@bind overlay_walls html"<input type='checkbox' checked> show walls"
# ╔═╡ 176b18dc-450b-11eb-20bd-c3ab5982c5b6
begin
field = similar(plan, RGB{Float64})
local ind = CartesianIndices(plan)
for i in 1:fsize
field[ind[i]] = cm[Ei[i]]
end
if overlay_walls
field[plan .== 0] .= RGB(0)
end
# show antenna position with a green crosshair
# first, the outline to make it more visible on top of the bright background
local ax, ay = antenna_pos
for d in -3:3
field[ax + d, ay - 1] = RGB(0, .2, 0)
field[ax + d, ay + 1] = RGB(0, .2, 0)
field[ax - 1, ay + d] = RGB(0, .2, 0)
field[ax + 1, ay + d] = RGB(0, .2, 0)
end
# next, the crosshair itself
for d in -3:3
field[ax + d, ay] = RGB(0, 1, 0)
field[ax, ay + d] = RGB(0, 1, 0)
end
field
end
# ╔═╡ e85b093c-4500-11eb-3441-c1d7431c634e
md"## et voilà!"
# ╔═╡ Cell order:
# ╟─14033bfc-4506-11eb-345d-b968dbdd669f
# ╟─b2375ede-4509-11eb-02c1-fde861ed36ad
# ╟─b76f6990-4509-11eb-1f2c-8933290c6292
# ╟─bc937a5e-4509-11eb-0cfa-c5634086f187
# ╟─1718edba-450a-11eb-2377-113207096427
# ╟─18908fcc-450a-11eb-18fb-53d6a09bd211
# ╠═1cf76d74-450a-11eb-3900-5dbbe000047a
# ╠═22808b72-450a-11eb-25ae-15d190b9c876
# ╠═3432951a-450a-11eb-0a87-05625cdca2f4
# ╠═87240306-450f-11eb-2b4b-93419692a783
# ╟─3531c2c2-450a-11eb-3e62-e31c5ad04546
# ╠═3d210916-450a-11eb-0ab4-4bde681a8fa2
# ╟─434cbaa6-450a-11eb-1e22-d50acd4eaabc
# ╠═476333fe-450a-11eb-1b11-5bf8767c8847
# ╟─4efa869e-450a-11eb-222f-79c99c8ee078
# ╠═55198d04-450a-11eb-0529-8f14d12e67cf
# ╟─6af3e020-450a-11eb-09ea-392a4fe9d03e
# ╠═73c7354e-450a-11eb-121d-0b694833a20b
# ╟─77c4f118-450a-11eb-336c-1d76291a980c
# ╟─7c03f8aa-450a-11eb-2006-83dfdb3c3527
# ╟─85ae4306-450a-11eb-02dd-7d43dc23b9e0
# ╟─8dbcf3a8-450a-11eb-10a5-1f0f89515bbb
# ╟─ab2a9cc4-450a-11eb-23bf-1f194db3142e
# ╟─b0e9f36a-450a-11eb-1bcf-af7ce83a9a0c
# ╠═cbf3f270-450a-11eb-0eba-0be7069aafc4
# ╟─ce7c4ef2-450a-11eb-0120-4dd2b5b6ebe2
# ╠═d466118e-450a-11eb-0260-e931d1fe5edb
# ╟─d76ed16a-450a-11eb-1ff6-13a784896b1f
# ╟─8a9acc22-4582-11eb-10d1-af27fb9fa40e
# ╠═a6105b1a-45a1-11eb-37bc-25e709d2aa61
# ╠═dd620024-450a-11eb-0d26-39444e86dbc7
# ╟─e13d4758-450a-11eb-21fa-21ef2b053b3c
# ╠═e53559e8-450a-11eb-10df-2f9bd44c829b
# ╟─ea215d76-450a-11eb-34c4-8f840a4c07e0
# ╠═ed05c25c-450a-11eb-3ca4-fd98065261f4
# ╟─f13d360c-450a-11eb-3ea5-a3c6738d3f7b
# ╟─8faa5ac4-4574-11eb-113f-49565b459121
# ╠═f4454b66-450a-11eb-1429-f3e5cd06ae3b
# ╟─f7acc81a-450a-11eb-3c0c-61badd6b8835
# ╟─faf621ce-450a-11eb-0158-398569973267
# ╠═ff2a96bc-450a-11eb-1d55-c1a1aab6dfc4
# ╟─039ef896-450b-11eb-231d-e1a15b86df0c
# ╠═07031f76-450b-11eb-1824-8367aa6aca36
# ╟─14232eda-450b-11eb-23a1-9731174997ad
# ╟─02251d32-4575-11eb-31de-cb4d694dd5b3
# ╠═176b18dc-450b-11eb-20bd-c3ab5982c5b6
# ╟─e85b093c-4500-11eb-3441-c1d7431c634e
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment