Try to calculate the position of the transmitter from the radio wave propagation model with python [Wi-Fi, Beacon]

• This blog is the first entry of the technical blog of Our Datawise.

Overview

Purpose: In general, position calculation using Wi-Fi or Beacon is something that you can understand what you are doing behind the scenes ~~.

Method: I will explain the basic idea and demonstrate it using Python code.

environment: MacOS Mojave 10.14.6, Python 3.7.4

basic way of thinking

Things necessary

1. Radio receivers (3) and their positions
2. Transmitter (smartphone, etc.) and its position
3. Strength of radio waves reaching each receiver
4. Okumura-Hata Curve

procedure

Step1. Calculate the distance with ** Okumura-Hata curve ** from the signal strength received from the transmitter and the position information of the receiver. Step2. Prepare a circle with the calculated distance as the radius and the position of each receiver as the center. Step3. Find the intersection of circles with ** 3-point positioning **. ** Coordinates of intersection = Position coordinates of transmitter **

A concrete image of three-point positioning is at the bottom of this entry. If you are interested, please see that first.

Code and results

Step0. Preparation of correct answer data

Prepare the correct answer data for the actual calculation. This time, we targeted Toranomon Hills and its surrounding facilities where we are located. The specific position is like this.

Now that we have these latitudes / longitudes Distance from the three-square theorem, Find the (theoretical) signal strength at each receiver from the distance.

This time, we are using Urban environments, small or middle size city from Hata-model: Wikipedia.

# -*- coding: utf-8 -*-
import math
import numpy as np
import pandas as pd
import sympy.geometry as sg

# latitude 1 digree = 0.0090133729745762 km
# longitutde 1 degree = 0.010966404715491394

def eucl_dist(pointA, pointB):
	return np.sqrt((pointA[0]-pointB[0])**2+(pointA[1]-pointB[1])**2)
	# euclidean distance

def Oku_Hata(f, h_b, C_H, d):
  Loss = (69.55 + 26.16*math.log10(f)
              - 13.82*math.log10(h_b) - C_H 
              + (44.9 - 6.55*math.log10(h_b))*math.log10(d))
  return Loss
  # Okumura_Hata Curve

def inv_Oku_Hata(f, h_b, C_H, Loss):
	numerator = Loss - 69.55 - 26.16*math.log10(f) + 13.82*math.log10(h_b) + C_H
	denominator = 44.9 - 6.55 * math.log10(h_b)
	d = 10**(numerator/denominator)
	return round(d, 10)
  # Inversed Okumura_Hata Curve

# variables 
f = 150  # frequency of transmissiton: Megahertz
h_b = 1  # height of base station antenna : meter
h_M =  1  # heigth of mobile station antenna ; meter
C_H =0.8 + (1.1 *math.log10(f) - 0.7) * h_M -1.56*math.log10(f)
d = 500

# Test whether or not the functions function?
euDist = eucl_dist([0,0],[3,4])
Loss = Oku_Hata(f, h_b, C_H, d)
d = inv_Oku_Hata(f, h_b, C_H, Loss) 

print('euDist should be 5: ', euDist)
print('Loss from distance :  ', Loss)
print('d should be 500: ', d)  # Should be 500

Check result

euDist should be 5: 5.0 Loss from distance : 248.5613025107495 d should be 500: 500.0

The ceremony seems to be done properly. Next, calculate the distance from the actual location information and convert it to Loss.

# Set of lattitude/ longitude each GeoData
tmitter = [35.6662344,139.7496597]  # Transmitter from Torano-Mon

sensorLocation = {'S1':['GoodMorningCafe',35.6664959,139.7500406],
                  'S2':['LeParc',35.6666179,139.7479538],
                  'S3':['GrandSuite',35.6675194,139.7489692]}
                  # SensorLocation around Trano-Mon

# Compute distance (= radius) from the transmitter
r1 = eucl_dist(tmitter, sensorLocation['S1'][1:])
r2 = eucl_dist(tmitter, sensorLocation['S2'][1:])
r3 = eucl_dist(tmitter, sensorLocation['S3'][1:])
radii = [r1, r2, r3]  # plural of radius

# Compute Loss
L1 = Oku_Hata(f, h_b, C_H, r1)
L2 = Oku_Hata(f, h_b, C_H, r2)
L3 = Oku_Hata(f, h_b, C_H, r3)

print('List of distance: ', radii)
print('List of Loss: ', L1, L2, L3)

result

List of distance: [0.0004620249560447818, 0.0017484756389397347, 0.0014587718293119327] List of Loss: -22.378972679965557 3.572964717331658 0.040582137429893805

Step1. Calculate the distance with ** Okumura-Hata curve ** from the signal strength received from the transmitter and the position information of the receiver.

Since Loss could be calculated (from distance) Calculate the distance (from Loss). ~~ I'm doing it with correct answer data, so it can't be helped if I go around. Don't tell me to use the distance from the beginning! ~~

d1 = inv_Oku_Hata(f, h_b, C_H, L1)
d2 = inv_Oku_Hata(f, h_b, C_H, L2)
d3 = inv_Oku_Hata(f, h_b, C_H, L3)

print('Computed Value: ', [d1, d2, d3])  # plural of radius
print('Original Value: ', radii)

Check result

Computed Value: [0.000462025, 0.0017484756, 0.0014587718] Original Value: [0.0004620249560447818, 0.0017484756389397347, 0.0014587718293119327]

Since the value is rounded in the middle of the formula, the number of digits is different, but it seems that it can be calculated properly.

Step2. Prepare a circle with the calculated distance as the radius and the position of each receiver as the center.

code. I used sympy because it seems to be fun to use. Qiita: The intersection of circles with sympy

# Step2.
centers = [sg.Point(sensorLocation['S1'][2], sensorLocation['S1'][1]),
           sg.Point(sensorLocation['S2'][2], sensorLocation['S2'][1]),
           sg.Point(sensorLocation['S3'][2], sensorLocation['S3'][1])]
           # x = longitude, y = latitude

circles = [sg.Circle(centers[0], radii[0]), 
           sg.Circle(centers[1], radii[1]),
           sg.Circle(centers[2], radii[2])]
           # define all three circles

print('List of centers (latitude/longitude) :', centers)
print('List of Circles :', circles)

result

Center coordinates

List of centers (latitude/longitude) : [Point2D(698750203/5000000, 356664959/10000000), Point2D(698739769/5000000, 356666179/10000000), Point2D(349372423/2500000, 178337597/5000000)]

Yen information. The radius has been added. It's a list of center coordinates and radii.

List of Circles : [Circle(Point2D(698750203/5000000, 356664959/10000000), 231012478022391/500000000000000000), Circle(Point2D(698739769/5000000, 356666179/10000000), 174847563893973/100000000000000000), Circle(Point2D(349372423/2500000, 178337597/5000000), 145877182931193/100000000000000000)]

Step3. Find the intersection of circles with ** 3-point positioning **. ** Coordinates of intersection = Position coordinates of transmitter **

code. After all, I use sympy. There are 3 circles, but since we are calculating the intersections by 2 each, we can get 6 solutions. Therefore, the logic is to extract the intersections that appear three times out of the six solutions. Therefore, it is rounded and the calculation error is ignored.

# Step3. 
crossPoints12 = sg.intersection(circles[0], circles[1])
crossPoints23 = sg.intersection(circles[1], circles[2])
crossPoints31 = sg.intersection(circles[2], circles[0])
# Compute crosspoints for each two circles

crossPointsList = [
    [float(crossPoints12[i].y),float(crossPoints12[i].x)] for i in [0,1]
    ] + [
    [float(crossPoints23[i].y),float(crossPoints23[i].x)] for i in [0,1]
    ] + [
    [float(crossPoints31[i].y),float(crossPoints31[i].x)] for i in [0,1]
    ]
crossPointsList = [[round(i[0],8), round(i[1],8)]for i in crossPointsList]
df = pd.DataFrame(crossPointsList).groupby([0, 1])[0].count()
print('Computed Value: ', df[df.values == max(df.values)].index[0])
print('Original Value: ', tmitter)
# print the most frequent set of lattitude-longitude

result

Computed Value: (35.6662344, 139.7496597) Original Value: [35.6662344, 139.7496597]

It seems that the correct answer was derived. I'm happy.

Supplement: Image of three-point positioning, etc.

If you google with "three-point positioning", you will see a lot, but for the time being, this blog was the oldest entry I've seen, so I'll introduce respect. Reference: I wish I could get the 2D coordinates of the iPhone with RSSI of iBeacon

So, when I talk about the end of the story, I'm talking about blogs other than the above blogs, In practice, this code, especially the Hata curve, is useless.

In the city, objects are scattered and radio waves are reflected diffusely, so radio wave propagation does not reach the receiver in the ideal state. I don't plan to introduce the detailed story of that area because it goes beyond the technical blog. If you want to know, I think it would be nice if you could google or join us.

that's all.

~~ I don't care, but I feel like I'm going to bite my tongue with Denpa Denpa ~~