-
-
Save vthriller/db8560f5b358d86a67a6f32b570acc4c to your computer and use it in GitHub Desktop.
Calculating WIFI propagation with IJulia
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
### 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