【搬运】基于WiFi RSSI定位的酒店房间查找
本帖最后由 vany5921 于 2020-6-8 11:59 编辑是否可以使用ESP32或ESP8266设备以足够的精度确定酒店房间的位置？答案很简单：是的，如果酒店周围有足够的WiFI接入点，有足够的时间在地板上走来走去几个小时，而且不会被他人起疑心，就可以做到。
背景：在酒店房间里闲逛会很无聊。电视节目没有什么可看的，在网络视频上所有的东西都已经被看了。好吧，然后我发现我可以在房间里检测到14个无线接入点。哇，太多了。在我楼层的楼梯上，我检测到16个，在走廊上，我能够检测到11到15个WiFi接入点之间的东西。这真的是太棒了，这个程序就是在这么无聊的情况下写出来的。
你知道吗：酒店的走廊又长又单调。如果你经常更换酒店，你往往忘记房间在哪个位置。所以问题来了：我的M5Stack和大量的WiFi接入点能帮助我快速方便地找到我的房间吗？下图是我房间的位置
答案取决于你如何定义“快速”和“轻松”，但基本上是这样的：
1.绘制地图
作为第一步，系统必须了解您房间在走廊上的位置。要执行此操作，请从房间门前开始，向左（正）或向右（负）移动。在每一步之后，让M5Stack扫描WiFi接入点，让他了解位置和WiFi信号强度之间的关系。为了尽可能宽地表示走廊，向左和向右移动是很有用的。最好在走廊上走来走去两次，但要提防其他酒店客人：他们宁愿相信你想犯罪，也不不相信你是在做实验。
https://hackster.imgix.net/uploads/attachments/1133993/walk_01_5KBZf7iGQ8.GIF?auto=compress&gifq=35&w=740&h=555&fit=max
GUI界面将引导您通过走廊。你只需要一步一步，让设备在每个位置扫描。重要的是你总是迈同样大小的步子。否则结果会不准确！我个人认为，横向搜索更容易，即使在搜索过程中我不得不再次横向搜索以避免测量差异。侧向行走更容易保持相同的步长。记住：你不想在浴室里用别人的牙刷在错误的房间里醒来！所以你最好花点时间做个详细的测量。
2.IILTM:逆向强度查找表
在内部，软件将5阶多项式拟合到学习到的数据。此功能映射走廊中的位置与每个检测到的WiFi接入点的信号强度（RSSI）之间的关系。
然而，对于走廊沿线的定位，则需要公式的另一个方向：位置必须根据WiFi接入点的信号强度确定。为了能够根据信号强度计算位置，将在走廊长度上创建反向强度查找表地图：
别担心，只要你在走廊里闲逛完，软件就会自动完成这一切（详情如下）。
3.出去喝一杯
现在什么都不会发生！你可以关掉M5Stack，去最近的酒吧转转。这个系统知道你的床在哪里，并帮助你在酒店无尽的走廊中找到自己的房间。
4.让设备带你回房间
好吧，你应该已经知道你的房间在几楼了。否则，你将不得不经历所有这些，你将失去从你必要的就寝时间。最好的办法是把一只旧袜子放在楼梯上标出正确的楼层（或任何不显眼的地方）。在走廊里，你可以问M5Stack你离房间有多远。设备扫描WiFi接入点，并通过反向强度查找表映射（IILTM）。对于每个条目，计算期望信号强度和测量信号强度之间差的平方和。此曲线的最小值表示您最有可能位于走廊中的位置。但是，如果显示的是一个很大的正数，那么您仍然离房间很远。
如果显示器显示的是负数，则表示您已经走过头了；如果显示器显示的是零，则表示您应该站在房间门的正前方
工作原理：
传统的室内定位方法是基于训练调节的方法。但要做到这一点，你必须知道发射器的确切位置。在酒店，很难精确测量所有WiFi接入点的位置。如果侧身在走廊游荡已经让人可疑，那么想想如果你把天花板的盖子取下来开始寻找电子设备会是什么样子？此外，为了得到距离和信号强度之间的可计算关系，信号必须是直线连接的。这与现代建筑相矛盾，在现代建筑中，这些是不可见的。
另一种常用的方法是使用指纹数据库根据信号强度的模式从信号强度推断位置。我的方法类似，但不同：
在现代酒店里，有许多WiFi接入点，为每一位客人提供互联网接入，无论他们身在何处。在我的酒店走廊里，我可以检测到许多接入点，当我沿着走廊走下去时，到接入点的距离会发生变化，检测到的信号强度也随之改变。换言之：走廊沿线的每个位置或多或少都有一个无障碍接入点信号强度的独特组合。然而，正是在WiFi信号的性质中，信号强度（RSSI）是强噪声的，这使得精确定位变得困难。此外，信号强度模式的基于快照的指纹也会受到信号噪声的影响。现在让数学帮我们实现。
信号形状
一般来说，这听起来很简单：离接入点越远，信号强度越低。所以我应该能从信号强度推断出距离，对吧？不幸的是，它有点复杂，信号的噪声并不能使它更容易。让我们看看信号强度和距离之间的相关性，以了解其复杂性：
当一个发射机在我的视线范围内时，距离和信号强度之间有明显的相关性。这种相关性不是线性的，但仍然是一条简单的曲线：
这看起来很有希望，因为我可以清楚地找到给定RSSI值的距离。然而，当我通过发射器到达我的房间时，就变得更加困难了。信号强度和我在走廊里的位置不再有明显的关系。我可能同时在走廊的两个地方
这两种说法挺有趣，但是这确实是个令人头痛的问题
如果发射机位于不同楼层或柱后或其他一些结构部件后面，则沿走廊的信号强度可能会变得更加复杂。这些东西不再是一个简单的函数，而且使你更难确定自己的位置。但是，如果走廊区域中有许多接入点，这些接入点希望随机分布，那么可以从大量数据中计算出更精确的位置信号。因此，我们的想法是，不要指望一个接入点，而是指望图片中许多的接入点。类似于以下示例，其中-5的距离由“红色”访问点的-72dBm值和“蓝色”访问点的-45dBm值明确给出：
但要做到这一点，对于酒店楼层沿线的接入点数量，需要对每个接入点的测量数据进行功能表示。
高阶多项式可用于确定位置与测量信号强度值之间的函数关系。像我的M5Stack中的ESP32这样的小型微控制器不太擅长执行高精度的数学计算，但我已经尝试创建一个库，它可以以可接受的精度拟合6次多项式。在高于6或7时，由于浮点运算的限制，结果变得越来越不准确。好吧，很公平。我最喜欢的酒店的测试数据表明，5阶多项式是足够的。
信号质量
曲线拟合库的拟合类允许您在运行时添加新的数据对。这允许多项式表示随时间“增长”，而无需将所有数据临时存储在RAM内存中。数据点越多，拟合越好。如果数据点太少，结果可能会变得混乱。因此，在接下来的步骤中，在走廊中的数据收集完成后，将忽略收集到的数据点少于6个的所有访问点的数据记录。此外，沿走廊行走时，所有信号电平振幅小于15dBm的接入点也将被忽略。
软件还记住了走廊中的最大位置。由此，现在创建了一个位置增量小于所做步骤的表。使用这些值可以创建具有所有可用访问点的道路的反向强度查找表地图。由于该图基于走廊沿线信号强度的函数表示，因此有效地滤除了信号强度的典型噪声。
此映射已保存，稍后将用于本地化。地图还可以传输到其他M5Stack设备，这样我的朋友就可以找到我的房间，例如，不用告诉他们我的房间号。好吧，好吧，也许告诉他们房间号会更容易，但没那么有趣！
最小二乘位置检测法
为了计算到我房间的距离，这个设备扫描了我站着的地方的无障碍接入点的信号强度。现在如果我站在相应的位置，对于反向强度查找表映射中的每个条目来说，它将检查访问点必须具有的信号强度。计算与测量信号强度之差并将其平方。对所有访问点都执行同样的操作，并对值进行求和。
如果我离强度图上的某个点很远，那么差的平方和会给出一个很高的值。但对于我所处的位置，总和很小，因为测量信号强度和理想信号强度之间的差异很小。因此，强度图上这个计算的最小值表明我很可能在哪里。
就这样。仅此而已：一个带有WiFi芯片的设备，一点数学知识，一个小GUI，几个小时就可以在大厅里上下走动。没别的了！在我最喜欢的酒店里，它工作得很好，但在一个只有酒吧才有WiFi的老棚子里，最好拍张房间门的照片，或者只是看看房间钥匙就知道号码。或者有人在接待处询问，并在一张纸上写下房间号。但是在一个有很多接入点的现代化酒店里，这个系统非常容易使用，即使这完全没有意义。
/**************************************************************************
* WiFi based location system to find home.
*
* A Indoor Positioning System based on WiFi RSSI data to return to the
* right room of a hotel corridor. Maybe useless, but it works just fine!
*
* Hague Nusseck @ electricidea
* v1.1 05.June.2020
* https://github.com/electricidea/M5Stack-Hotel-room-finder
*
* Changelog:
* v1.1 = - first published version
*
* Distributed as-is; no warranty is given.
*
***************************************************************************/
#include <Arduino.h>
#include <M5Stack.h>
// install the library:
// pio lib install "M5Stack"
#include "WiFi.h"
// Free Fonts for nice looking fonts on the screen
#include "Free_Fonts.h"
// library for liniear and nonlinear fits
#include "curve_fit.h"
// maximum number of access points = maximum number of fits
const int max_fits = 40;
// array of of a number of fits
curve_fit fits = curve_fit();
// position on the floor
double min_pos = 99999;
double max_pos = -99999;
// array for the calculated x positions along the floor
double *newx_array;
int n_newx = 0;
// the inverse intensity lookup table Map
double *IILTM;
int n_usable_APs = 0;
// the BSSID lookup table
typedef char cstring;
cstring *BSSIDLT;
// the array for the square sums
double *square_sum_array;
// state machine index to switch between the menu states
int menu_state = 0;
#define STATE_START 1
#define STATE_MEASURE 2
#define STATE_GET_DATA 3
#define STATE_DATA 4
#define STATE_RUN 5
// value for the measurement along the floor
int measure_position = 0;
//==============================================================
// function forward declaration
uint16_t RGB2Color(uint8_t r, uint8_t g, uint8_t b);
void writeFile(fs::FS &fs, const char * path, const char * message);
void Clear_Screen();
void print_menu(int menu_index);
String split(String source, char delimiter, int location);
uint8_t collect_WiFi_data(String filename, bool append = true);
bool load_measurement(String filename);
bool load_floor_data();
String calculate_position();
bool analyze_measurements();
void setup() {
// initialize the M5Stack object
M5.begin();
// configure the Lcd display
M5.Lcd.setBrightness(100); //Brightness (0: Off - 255: Full)
M5.Lcd.setTextColor(TFT_WHITE);
M5.Lcd.setTextSize(1);
Clear_Screen();
// configure centered String output
M5.Lcd.setTextDatum(CC_DATUM);
M5.Lcd.setFreeFont(FF2);
M5.Lcd.drawString("Hotel room Finder", (int)(M5.Lcd.width()/2), (int)(M5.Lcd.height()/2), 1);
M5.Lcd.setFreeFont(FF1);
M5.Lcd.setTextColor(TFT_RED);
M5.Lcd.drawString("Version 1.1 | 05.06.2020", (int)(M5.Lcd.width()/2), (int)(M5.Lcd.height()/2)+50, 1);
// configure Top-Left oriented String output
M5.Lcd.setTextDatum(TL_DATUM);
M5.Lcd.setTextColor(TFT_WHITE);
// Start Menu
menu_state = 1;
print_menu(menu_state);
}
void loop() {
M5.update();
// left Button
if (M5.BtnA.wasPressed()){
switch (menu_state) {
case STATE_START: { //START -> MEASURE
Clear_Screen();
menu_state = STATE_MEASURE;
print_menu(menu_state);
break;
}
case STATE_MEASURE: { //MEASURE -> ADD
Clear_Screen();
M5.Lcd.println("Stand in front of the door\nand face the door.\n");
M5.Lcd.println("got to the LEFT and press (<)");
M5.Lcd.println("or to the RIGHT and press (>)");
measure_position = 0;
menu_state = STATE_GET_DATA;
print_menu(menu_state);
break;
}
case STATE_GET_DATA: { //NEW -> < (-)
Clear_Screen();
++measure_position;
M5.Lcd.printf("measure %i steps away\n", measure_position);
int n_WiFi_networks = collect_WiFi_data("/WiFi_data.txt");
if(n_WiFi_networks > 0)
M5.Lcd.printf(" %i Networks found\n", n_WiFi_networks);
else
M5.Lcd.println(" unable to scan WiFi");
print_menu(menu_state);
break;
}
case STATE_DATA: { //DATA -> DELETE
Clear_Screen();
M5.Lcd.println("Delete all measured data...");
if(SD.remove("/WiFi_data.txt"))
M5.Lcd.println(" data deleted");
else
M5.Lcd.println("\n\n unable to deleted data");
writeFile(SD, "/WiFi_data.txt","pos;n;name;id;RSSI");
measure_position = 0;
n_usable_APs = 0;
n_newx = 0;
min_pos = 99999;
max_pos = -99999;
print_menu(menu_state);
break;
}
case STATE_RUN: { //RUN -> CHECK
// check my position
Clear_Screen();
M5.Lcd.setTextDatum(CC_DATUM);
M5.Lcd.setFreeFont(FF3);
M5.Lcd.drawString("Let me check...", (int)(M5.Lcd.width()/2), (int)(M5.Lcd.height()/2), 1);
String pos_result = calculate_position();
Clear_Screen();
M5.Lcd.setTextDatum(CC_DATUM);
M5.Lcd.setFreeFont(FF4);
M5.Lcd.drawString(pos_result.c_str(), (int)(M5.Lcd.width()/2), (int)(M5.Lcd.height()/2), 1);
print_menu(menu_state);
break;
}
}
}
// center Button
if (M5.BtnB.wasPressed()){
switch (menu_state) {
case STATE_START: { //START -> RUN
Clear_Screen();
M5.Lcd.setTextDatum(TL_DATUM);
M5.Lcd.setFreeFont(FF1);
// load floor data from SD card
if(load_floor_data()) {
Clear_Screen();
M5.Lcd.println("\n\n OK, ready to run");
menu_state = STATE_RUN;
} else
M5.Lcd.println("Failed to load data");
print_menu(menu_state);
break;
}
case STATE_MEASURE: { //MEASURE -> BACK
Clear_Screen();
menu_state = STATE_START;
print_menu(menu_state);
break;
}
case STATE_GET_DATA: { //MEASURE -> DONE
Clear_Screen();
// analyze measured data
M5.Lcd.println("let's analyze the data");
if(analyze_measurements())
M5.Lcd.println("\n\n OK, Success!");
else
M5.Lcd.println("\nSorry\nSomething went wrong.. :-(");
menu_state = STATE_START;
print_menu(menu_state);
break;
}
case STATE_DATA: { //DATA -> BACK
Clear_Screen();
menu_state = STATE_START;
print_menu(menu_state);
break;
}
case STATE_RUN: { //RUN -> DONE
Clear_Screen();
menu_state = STATE_START;
print_menu(menu_state);
break;
}
}
}
// right Button
if (M5.BtnC.wasPressed()){
switch (menu_state) {
case 1: { //START -> DATA
menu_state = STATE_DATA;
print_menu(menu_state);
break;
}
case STATE_MEASURE: { //MEASURE -> NEW
Clear_Screen();
M5.Lcd.println("Delete all measured data...");
if(!SD.remove("/WiFi_data.txt"))
M5.Lcd.println(" unable to deleted data");
writeFile(SD, "/WiFi_data.txt","pos;n;name;id;RSSI");
M5.Lcd.println("\nReady for new measurements");
M5.Lcd.println("\nStand in front of the door\nand face the door.\n");
M5.Lcd.println("got to the LEFT and press (<)");
M5.Lcd.println("or to the RIGHT and press (>)");
measure_position = 0;
menu_state = STATE_GET_DATA;
print_menu(menu_state);
break;
}
case STATE_GET_DATA: { //NEW ->> (+)
Clear_Screen();
--measure_position;
M5.Lcd.printf("measure %i steps away\n", measure_position);
int n_WiFi_networks = collect_WiFi_data("/WiFi_data.txt");
if(n_WiFi_networks > 0)
M5.Lcd.printf(" %i Networks found\n", n_WiFi_networks);
else
M5.Lcd.println("\n\n unable to scan WiFi");
print_menu(menu_state);
break;
}
case STATE_DATA: { //DATA -> INFO
Clear_Screen();
M5.Lcd.printf("Data Info:\n\n");
M5.Lcd.printf("usable APs: %i\n", n_usable_APs);
M5.Lcd.printf("min x pos: %.1f\n", min_pos);
M5.Lcd.printf("max x pos: %.1f\n", max_pos);
print_menu(menu_state);
break;
}
case STATE_RUN: { //RUN -> Nothing
break;
}
}
}
delay(50);
}
//==============================================================
// convert a RGB color into a LCD color value
uint16_t RGB2Color(uint8_t r, uint8_t g, uint8_t b) {
return ((r / 8) << 11) | ((g / 4) << 5) | (b / 8);
}
//==============================================================
// split() returns a substring out of a source string
// defined by a delimiter and a location
// inspired by:
// http://wp.scalesoft.de/arduino-split/
String split(String source, char delimiter, int location) {
String result = "";
int locationCount = 0;
int FromIndex = 0, ToIndex = -1;
while (location >= locationCount) {
FromIndex = ToIndex + 1;
ToIndex = source.indexOf(delimiter,FromIndex);
// if there is no delimiter anymore within the source string
if(ToIndex == -1){
// return the end of the string starting from the last delimiter
// if the actual location is the target location
if (location == locationCount) {
result = source.substring(FromIndex);
} else {
// if not, return a empty string
return "";
}
} else {
// if it is not the last delimiter, but the actual
// loacation is the target location, return the sub String
if (location == locationCount)
result = source.substring(FromIndex,ToIndex);
}
++locationCount;
}
return result;
}
//==============================================================
// Scan for WiFi networks and save the SSID, BSSID and RSSI
// into a file on the SD card
// fix file name: "pos_data.txt"
uint8_t collect_WiFi_data(String filename, bool append){
File file;
if(append)
file = SD.open(filename.c_str(), FILE_APPEND);
else
file = SD.open(filename.c_str(), FILE_WRITE);
int n = WiFi.scanNetworks();
if (n == 0) {
M5.Lcd.println(" no networks found");
} else {
for (int i = 0; i < n; ++i) {
// Print SSID, BSSID and RSSI for each network found
file.printf("%i;%i;%s;%s;%i\n",measure_position , i+1, WiFi.SSID(i).c_str(), WiFi.BSSIDstr(i).c_str(), WiFi.RSSI(i));
delay(10);
}
}
file.close();
return n;
}
//==============================================================
// Write Text into a file
void writeFile(fs::FS &fs, const char * path, const char * message){
File file = fs.open(path, FILE_APPEND);
if(!file){
M5.Lcd.println("\n\n Failed to open file");
return;
} else {
if(!file.println(message))
M5.Lcd.println("\n\n Write failed");
file.close();
}
}
//==============================================================
// Print a small menu at the bottom of the display above the buttons
void print_menu(int menu_index){
M5.Lcd.fillRect(0, M5.Lcd.height()-25, M5.Lcd.width(), 25, RGB2Color(50,50,50));
M5.Lcd.setCursor(0, 230);
M5.Lcd.setFreeFont(FF1);
M5.Lcd.setTextColor(TFT_WHITE);
switch (menu_index) {
case 0: { // never used
M5.Lcd.print(" - - - ");
break;
}
case STATE_START: { // start menu
M5.Lcd.print(" MEASURE RUN DATA");
break;
}
case STATE_MEASURE: { // Submenu 1 for new Data collection
M5.Lcd.print(" ADD BACK NEW ");
break;
}
case STATE_GET_DATA: { // Submenu 2 for new Data collection
M5.Lcd.print(" < DONE > ");
break;
}
case STATE_DATA: { // DATA Submenu
M5.Lcd.print(" DELETE BACK INFO");
break;
}
case STATE_RUN: { // RUN Submenu
M5.Lcd.print(" CHECK DONE ");
break;
}
default: { // should never been called
M5.Lcd.print(" - - - ");
break;
}
}
}
//==============================================================
// Clear the entire screen and add one row
// The added row is important. Otherwise the first row is not visible
void Clear_Screen(){
M5.Lcd.fillScreen(BLACK);
M5.Lcd.setCursor(0, 0);
M5.Lcd.println("");
}
//==============================================================
// loads a stored measurement of positions and BSSID, RSSI data
// the data is used to learn the fits for each WiFi access point
// fix file name: "WiFi_data.txt"
bool load_measurement(String filename){
// reset min and max for new measurements
min_pos = 99999;
max_pos = -99999;
// reset all fits
// and set the tag to -1 = not learned
for(int i = 0; i < max_fits; ++i){
// init as fith order polynomial
fits.init(5);
fits.reset();
fits.tag = -1;
}
File file = SD.open("/WiFi_data.txt");
if(!file){
M5.Lcd.println("Failed to open file");
} else {
String line = "";
while(file.available()){
// file format:
// pos;n;name;id;RSSI
char chread = file.read();
// ASCII printable characters (character code 32-127)
if((chread >= 32) && (chread <= 127))
line = line + String(chread);
// CRLF = \r\n = ASCII code 13 and then ASCII code 10
if((chread == '\r' || chread == '\n') && (line.length() > 0)){
// get the BSSID at index 3
String BSSID = split(line, ';', 3);
// skip the header line
if(BSSID != "id"){
// get the position at index 0
double pos = split(line, ';', 0).toDouble();
if(pos > max_pos)
max_pos = pos;
if(pos < min_pos)
min_pos = pos;
// get the RSSI at index 4
double RSSI = split(line, ';', 4).toDouble();
int i = 0;
bool OK_Flag = false;
// check all fits if the BSSID exist
// then let this fit learn the new data pair
// otherwise, use the next unused fit for this BSSID
do {
if(fits.name == BSSID){
fits.learn(pos, RSSI);
OK_Flag = true;
}
if(fits.tag == -1){
fits.tag = i;
fits.name = BSSID;
Serial.print(i);
Serial.print(": ");
Serial.println(BSSID);
fits.learn(pos, RSSI);
OK_Flag = true;
}
if(++i == max_fits){
OK_Flag = true;
}
} while (!OK_Flag);
}
line = "";
}
}
file.close();
return true;
}
return false;
}
//==============================================================
// loads a stored floor data from SD card
// the data is used to find the room
// fix file name: "floor_data.txt"
bool load_floor_data(){
// load floor data from file
M5.Lcd.printf("loading from file:\n-->/floor_data.txt\n");
File file = SD.open("/floor_data.txt");
if(!file){
return false;
} else {
// File format:
// n_newx
// n_usable_APs
// newx_array ... newx_array
// BSSIDLT ... BSSIDLT
// IILTM ... IILTM[((n_usable_APs-1)*n_newx)+(n_newx-1)]
String line = "";
int File_Block_index = 0;
// 0 = header information with array dimensions
// 1 = newx_array
// 2 = BSSIDLT
// 3 = IILTM
int line_count = 0;
while(file.available()){
char chread = file.read();
if(chread != '\n'){
line = line + String(chread);
} else {
switch (File_Block_index) {
// 0 = header information with array dimensions
case 0:
n_newx = split(line, ';', 0).toInt();
n_usable_APs = split(line, ';', 1).toInt();
// free all daynamic arrays
free(newx_array);
free(square_sum_array);
free(BSSIDLT);
free(IILTM);
M5.Lcd.printf("new_x array size: %i \n", n_newx);
M5.Lcd.printf("n_usable_APs: %i \n", n_usable_APs);
Serial.printf("new_x array size: %i \n", n_newx);
Serial.printf("n_usable_APs: %i \n", n_usable_APs);
// allocate all dynamic arrays with the dimensions from the file
newx_array = (double*) malloc(n_newx*sizeof(double));
square_sum_array = (double*) malloc(n_newx*sizeof(double));
BSSIDLT = (cstring*) malloc(n_usable_APs*sizeof(cstring));
IILTM = (double*) malloc(n_newx*n_usable_APs*sizeof(double));
// stop, if allocation fails
if(!newx_array || !square_sum_array || !IILTM || !BSSIDLT){
M5.Lcd.printf(" unable to allocate memory\n");
delay(5000);
// stop processing readed data from file
File_Block_index = -1;
file.close();
return false;
} else{
Serial.println("done: header!");
Serial.println("read newx_array data...");
File_Block_index = 1;
}
break;
// 1 = newx_array data
case 1:
if(line != ""){
// read the new-x values line by line
newx_array = split(line, ';', 0).toDouble();
if(line_count == n_newx){
Serial.println("done: read newx_array!");
Serial.println("read BSSIDLT data...");
File_Block_index = 2;
line_count = 0;
}
}
break;
// 2 = BSSIDLT data
case 2:
if(line != ""){
// read the BSSID values line by line
strcpy(BSSIDLT, split(line, ';', 0).c_str());
if(line_count == n_usable_APs){
Serial.println("done: read BSSIDLT!");
Serial.println("read IILTM data...");
File_Block_index = 3;
line_count = 0;
}
}
break;
// 3 = IILTM data
case 3:
if(line != ""){
// read the IILTM data
for(int i=0; i < n_usable_APs; ++i){
IILTM[(i*n_newx)+line_count] = split(line, ';', i).toDouble();
}
++line_count;
if(line_count == n_newx){
Serial.println("done: read IILTM!");
File_Block_index = -1;
line_count = 0;
}
}
break;
default:
break;
}
line = "";
}
}
file.close();
}
return true;
}
//==============================================================
// scan fout times the available APS
// Calculate the best fitting positon based on the IILTM
// Return the number as text, or text if the position can't be calculated
String calculate_position(){
String result = "...";
// scan APs and save to File "pos_data.txt"
// four times...
int n_WiFi_networks = collect_WiFi_data("/pos_data.txt", false);
collect_WiFi_data("/pos_data.txt");
collect_WiFi_data("/pos_data.txt");
collect_WiFi_data("/pos_data.txt");
// open the file and go throw all APs...
File file = SD.open("/pos_data.txt");
if(!file || n_newx == 0 || n_usable_APs == 0 || n_WiFi_networks == 0){
// without a file or without any APs, we are unable to find the room
return "No idea :-(";
} else {
// Because we scanned four times, we have to average the RSSI data
// This can be done with the fit-class ()
// reset all fits
// and set the tag to -1 = not learned
for(int i = 0; i < max_fits; ++i){
fits.init(0);
fits.reset();
fits.tag = -1;
}
String line = "";
while(file.available()){
char chread = file.read();
// ASCII printable characters (character code 32-127)
if((chread >= 32) && (chread <= 127))
line = line + String(chread);
// CRLF = \r\n = ASCII code 13 and then ASCII code 10
if((chread == '\r' || chread == '\n') && (line.length() > 0)){
// file format:
// pos;n;name;BSSID;RSSI
String BSSID = split(line, ';', 3);
// skip the header line
if(BSSID != "id"){
double RSSI = split(line, ';', 4).toDouble();
// find AP in BSSIDLT
int AP_index = -1;
for(int i = 0; i < n_usable_APs; ++i){
if(strcmp(BSSIDLT, BSSID.c_str()) == 0)
AP_index = i;
}
if(AP_index > -1){
fits.learn(0.0, RSSI);
fits.tag = AP_index;
}
}
line = "";
}
}
file.close();
// Now, the fits are filled with the average RSSI data from the APs
// Time to calculate the square sum array:
for(int x = 0; x < n_newx; ++x)
square_sum_array = 0.0;
for(int AP_index = 0; AP_index < n_usable_APs; ++AP_index){
if(fits.tag > -1){
for(int x = 0; x < n_newx; ++x){
// predict(0.0) will return the average value because
// the fit is initialized with degree = 0
double RSSI_diff = (fits.predict(0.0) - IILTM[(AP_index*n_newx)+x]);
square_sum_array += (RSSI_diff*RSSI_diff);
}
}
}
// Finally find the minimum of the square sums:
double best_pos = newx_array;
double min_sum = square_sum_array;
for(int x = 0; x < n_newx; ++x){
if(square_sum_array < min_sum){
min_sum = square_sum_array;
best_pos = newx_array;
}
}
// If best pos is the first or the last position of the new x array
// then we can say that we are far away, because we might don't know
// the right value of the distance
if(best_pos == newx_array || best_pos == newx_array)
result = "far away...";
else
result = String(best_pos).c_str();
}
return result;
}
//==============================================================
// load the measurements from SD card (file: /WiFi_data.txt)
// build the new-x array, the BSSIDLT and the IILTM
// return true if the procedure was succesfull
bool analyze_measurements(){
bool result = false;
M5.Lcd.printf("Reading file:\n --> /WiFi_data.txt\n");
// load measured data from file and let the fits learn...
if(load_measurement("/WiFi_data.txt")){
M5.Lcd.printf("Analyze AP data\n");
// build new_x array....
// free the daynamic arrays
free(newx_array);
free(square_sum_array);
// get the position range out of the data
int x_range = round(max_pos - min_pos);
n_newx = x_range *2;
// allocate the memory for the array with the new size
newx_array = (double*) malloc(n_newx*sizeof(double));
square_sum_array = (double*) malloc(n_newx*sizeof(double));
// if the memory allocation failed
if(!newx_array || !square_sum_array)
return false;
else {
// fill the newx_array with the fine position steps
for(int i=0; i < n_newx; ++i){
newx_array = min_pos + (i*((double)x_range / (double)n_newx));
}
// check for usable APs out of the fits
// criteria:
// at least 6 valid data points
// --> fith order polynome should have at least 6 values
// estimated min and max y values should not be out of bounds [-25 .. -95]
// a minimum of 15dBm amplitude over the data range is required
n_usable_APs = 0;
for(int i = 0; i < max_fits; ++i){
// check all fits for criteria
double min_y = fits.estimate_min_y();
double max_y = fits.estimate_max_y();
if(fits.tag > -1){
if((fits.count() < 6) ||
(min_y < -95.0) || (max_y > -25.0) ||
(fabs(max_y - min_y) < 15)) {
fits.reset();
fits.name = "";
fits.tag = -1;
}
}
if(fits.tag > -1){
Serial.printf("%i: N: %i min: %.2f max: %.2f\n", i, fits.count(), fits.estimate_min_y(), fits.estimate_max_y());
++n_usable_APs;
}
}
// build Inverse Intensity Lookup Table Map (IILTM)
// ....
M5.Lcd.printf("Build IILTM and BSSIDLT\n");
Serial.println("Build the IILTM and the BSSIDLT:");
if(n_usable_APs > 0){
Serial.printf("number of usable APs: %i \n", n_usable_APs);
// free the daynamic arrays
free(IILTM);
free(BSSIDLT);
// allocate the memory for the array with the new size
BSSIDLT = (cstring*) malloc(n_usable_APs*sizeof(cstring));
IILTM = (double*) malloc(n_newx*n_usable_APs*sizeof(double));
// if the memory allocation failed
if(!IILTM || !BSSIDLT) {
Serial.println(" malloc failed");
return false;
} else {
int AP_count = 0;
for(int i = 0; i < max_fits; ++i){
//Serial.printf("check: %i: \n", i);
if(fits.tag > -1){
//Serial.printf("tag > -1 --> %i: %s\n", AP_count, fits.name.c_str());
strcpy(BSSIDLT, fits.name.c_str());
for(int x = 0; x < n_newx; ++x){
// -95dBm for x values outside the learned range
IILTM[(AP_count*n_newx)+x] = fits.predict(newx_array, -95.0);
}
++AP_count;
}
}
}
} else {
M5.Lcd.printf("no usable APs found!\n");
Serial.println("no usable APs found!");
return false;
}
Serial.println("the BSSIDLT:");
for(int i = 0; i < n_usable_APs; ++i){
Serial.printf("%i: %s\n", i, BSSIDLT);
}
Serial.println("the IILTM:");
for(int x = 0; x < n_newx; ++x){
Serial.printf("\n%.2f", newx_array);
for(int i = 0; i < n_usable_APs; ++i){
Serial.printf(" %.2f", IILTM[(i*n_newx)+x]);
}
}
// save floor data to file
M5.Lcd.printf("Writing to file:\n --> /floor_data.txt\n");
File file = SD.open("/floor_data.txt", FILE_WRITE);
if(!file){
M5.Lcd.println("Failed to open file");
return false;
} else {
// File format:
// n_newx
// n_usable_APs
// newx_array ... newx_array
// BSSIDLT ... BSSIDLT
// IILTM ... IILTM[((n_usable_APs-1)*n_newx)+(n_newx-1)]
// header with dimensions
file.printf("%i;%i\n",n_newx, n_usable_APs);
// save newx_array
for(int x=0; x < n_newx; ++x){
file.printf("%.6f\n",newx_array);
}
// save BSSIDLT array
for(int i=0; i < n_usable_APs; ++i){
file.printf("%s\n",BSSIDLT);
}
// save IILTM array
for(int x=0; x < n_newx; ++x){
file.printf("\n%.6f",IILTM[(0*n_newx)+x]);
for(int i=1; i < n_usable_APs; ++i){
file.printf(";%.6f",IILTM[(i*n_newx)+x]);
}
}
file.printf("\n");
file.close();
}
M5.Lcd.println("done..");
Serial.println("");
Serial.println("done!");
}
result = true;
}
return result;
}
/***************************************************
*
* Fitting of a Polynomial using Cramer Rule
*
* Hague Nusseck @ electricidea
*
* --> see curve_fit.cpp for more details on how it
* works and how to use it.
*
* Distributed as-is; no warranty is given.
*
* ************************************************/
// class definition
class curve_fit {
public:
curve_fit(uint8_t degree=2);
bool init(uint8_t degree);
void learn(double x, double y);
double predict(double x);
double predict(double x, double outside_value);
void get_coefficients(double values[]);
String get_formula(uint8_t decimals = 6);
uint32_t get_order();
void reset();
double max_x() const { return max_x_; }
double min_x() const { return min_x_; }
int count() const { return N; }
double estimate_max_y(uint32_t steps = 100);
double estimate_min_y(uint32_t steps = 100);
int tag;
String name;
private:
// Mindex generates the index for adressing the matrices
uint32_t Mindex(uint32_t i, uint32_t j);
double determinant(double *mainmatrix);
int order = -1;
double *M;
double *b;
// y = a*x^n + ... + a*x + a
double *a;
// Number of learned x, y pairs
uint32_t N = 0;
double max_x_, min_x_;
};
/**************************************************************************
* Library to calculate polynomial regression in Arduino.
*
* Simple, efficient and memory efficient method to
* calculate mean values, linear regressions or
* parabolas from pairs of values.
*
* Based on Cramers rulefor fitting of a polynomial
* regression model using least square method.
* For a very good explanation, please see:
* https://neutrium.net/mathematics ... ng-of-a-polynomial/
*
* and thanks to Cubiwan:
* https://github.com/cubiwan/LinearRegressino
* and Chandu yadav:
* https://www.tutorialspoint.com/c ... rminant-of-a-matrix
* for the inspiration on how to implement it
* and Martin R.
* https://codereview.stackexchange ... g-gauss-elimination
* for the tip about partial pivoting...
*
*
* Hague Nusseck @ electricidea
* v1.4 24.May.2020
* https://github.com/electricidea/xxxx
*
* Changelog:
* v1.2 = - first version ready with "Laplace expansion" Method to calculate determinants
* v1.3 = - Gaussian triangular method for the determination of determinants (faster!)
* - added name and tag
* - added max_x and min_x
* - established aij Matrix notation
* - String output for formular
* v1.4 = - dynamic matrix sizes according to initialized order
* - functions to estimate max and min y values over the known range of x
* - add function count() to return N
*
*
* Distributed as-is; no warranty is given.
*
**************************************************************************
*
* ==== How to use it: ====
*
* 1.) intitialize the fitting object:
*
* curve_fit fit_1 = curve_fit(1);
* NOTE:
* degree = 0 --> A constant (or: The Average of y values)
* degree = 1 --> A first order polynomial (or: linear regression)
* degree = 2 --> A second order polynomial (or: parabolic curve)
* degree = 3 --> A third order polynomial
* Note:
* Generally the degree is not limited, but from the 6th or 7th
* degree the floating-point arithmetic reaches its limits and
* there may be inaccuracies.
*
* In principle, there is no limit for the degree of polynomials,
* but from the 6th or 7th degree on, first calculation errors appear
* due to the inaccurate floating-point arithmetic.
*
* 2.) Let the function solve the polynomial regression model
* by adding pairs of x and y values:
*
* fit_1.learn(0,5);
* fit_1.learn(1,18.65);
* ...
* fit_1.learn(10,52.4);
*
* During each of these steps the polynomial representation is calculated.
* The value pairs are not stored. As a result, the computing time is
* higher, but the class can be used more flexible.
*
* 3.) print out the calculated coefficients:
*
* double coeff_values;
* fit_1.get_coefficients(coeff_values);
* Serial.printf("y= a1*x + a0\n\na1= %.2f\na0= %.2f",
* coeff_values, coeff_values);
*
* 4.) Calculate the y values for a new x value:
*
* new_y = fit_1.predict(new_x);
*
* Some NOTES:
* At any time additional pairs of values can be added to improve
* the calculation / allow further learning.....
* The value pairs are not stored. They cannot be read out again
* or changed afterwards.
*
*
* ==== How the Math works: ====
*
* Representation of a kth order polynomial
* y = ak*x^k +⋯+ a1*x + a0
*
* with
* a : Polynomial coefficient
* k : The degree of the polynomial
* N : The number of points to be regressed
*
* The coefficients of the polynomial regression model (ak,ak−1,⋯,a1) may be
* determined by solving the following system of linear equations:
*
* | | | | | |
* | N SUM(xi) ... SUM(xi^k)| | a0 | | SUM(yi) |
* |SUM(xi) SUM(xi^2 ) ... SUM(xi^k+1) | | a1 | |SUM(xi*yi)|
* | ... ... ... ... | * | .. | = | ... |
* | SUM(xi^k)SUM(xi^k+1) ... SUM(xi^2k) | | ak | | SUM(xi^k*yi) |
* | | | | | |
*
* Cramer's rule allows you to solve the linear system of equations to find the
* regression coefficients using the determinants of the square matrix M. Each
* of the coefficients ak may be determined using the following equation:
*
* ak = det(Mi) / det(M)
*
* Where Mi is the matrix M with the ith column replaced with the column
* vector b (remembering the system is presented in the form Ma=b).
* For example M0 could be calculated as follows:
*
*
* | |
* | SUM(yi) SUM(xi) ... SUM(xi^k)|
* | SUM(xi*yi) SUM(xi^2 ) ... SUM(xi^k+1) |
* M0 = | ... ... ... ... |
* | SUM(xi^k*yi) SUM(xi^k+1) ... SUM(xi^2k) |
* | |
*
* Example for an linear regression (polynominal 1th order)
*
*y = a1*x + a0
*
* | | | | | |
* | N SUM(xi) | | a0 | | SUM(yi)|
* | SUM(xi)SUM(xi^2)| * | a1 | = | SUM(xi*yi) |
* | | | | | |
*
*
* | |
* | N SUM(xi) |
*M = | SUM(xi)SUM(xi^2)|
* | |
*
* | |
* | SUM(yi) SUM(xi) |
* M0 = | SUM(xi*yi)SUM(xi^2)|
* | |
*
* | |
* | N SUM(yi) |
* M1 = | SUM(xi)SUM(xi*yi)|
* | |
*
* a0 = det(M0) / det(M)
* a1 = det(M1) / det(M)
*
* How to calculate the determinants of a 2x2 matrix M:
*
* | |
* | a11 a12 |
* det | a21 a22 | = a11*a22 – a12*a21
* | |
*
***************************************************************************/
#include "Arduino.h"
#include "curve_fit.h"
//==============================================================
// the constructor
curve_fit::curve_fit(uint8_t degree) {
tag = 0;
name = "";
curve_fit::init(degree);
}
//==============================================================
// Initialization of the fit object
// degree = 0 --> A constant (or: The Average of y values)
// degree = 1 --> A first order polynomial (or: linear regression)
// degree = 2 --> A second order polynomial (or: parabolic curve)
// degree = 3 --> A third order polynomial
// Note:
// Generally the degree is not limited, but from the 6th or 7th
// degree the floating-point arithmetic reaches its limits and
// there may be inaccuracies.
// this function can be used to change the degree at runtime
bool curve_fit::init(uint8_t degree) {
order = degree;
free(a);
free(b);
free(M);
// the matrix dimension need to order + 1 for the calculations
a = (double*) malloc((order+1)*sizeof(double));
b = (double*) malloc((order+1)*sizeof(double));
M = (double*) malloc((order+1)*(order+1)*sizeof(double));
// if malloc fails, set order to -1
// and return false
if(!a || !b || !M) {
order = -1;
N = 0;
return false;
}
curve_fit::reset();
return true;
}
//==============================================================
// clear all calculated coefficients and buffered values
// set all values back to 0.0 to start a new calculation
void curve_fit::reset() {
// initialize the matrices with 0.0
for(int i = 0; i <= order; i++) {
for(int j = 0; j <= order; j++)
M = 0.0;
a = 0.0;
b = 0.0;
}
// reset the number of used x, y pairs
N = 0;
max_x_ = 0.0;
min_x_ = 0.0;
}
//==============================================================
// calculate the determinant by following Gaussian Method:
// Use elementary row operations to bring the matrix into
// a row echelon form, to be able to calculated the determinant
// by using the diagonal values.
double curve_fit::determinant( double *mainmatrix) {
double det = 1;
// order < 2 = matrix dimention up to 2x2
// for these matrix sizes, the determinants can be calculate directly
if (order < 2) {
if(order == -1)
// A matrix without dimension?
// would say the deteminat is zero....
det = 0;
if(order == 0)
// The determinant of a 1x1 matrix is the element of the matrix itself
det = mainmatrix;
if(order == 1)
// a 2x2 matrix can be calculated directly:
// | | | |
// | a11 a12 | | MM |
// det | a21 a22 | = a11*a22 – a12*a21 = | MM | = M*M - M*M
// | | | |
det = ((mainmatrix * mainmatrix) - (mainmatrix * mainmatrix));
} else {
// otherwise, elementary row operations are required to bring the
// matrix into a row echelon form.
// First, create a work-matrix with the same dimensions
// and fill it with the data from the main matrix
double workmatrix[(order+1)*(order+1)];
for(int i = 0; i < order+1; ++i)
for(int j = 0; j < order+1; ++j)
workmatrix = mainmatrix;
// here are some hints about Gaussian Elimination with Partial Pivoting:
// https://youtu.be/euIXYdyjlqo
// https://www.youtube.com/watch?v=SQ98-eImrgA
// Now, go throw each column
for (int i = 0; i < order+1; ++i) {
// find the largest element for partial pivoting
double pivotElement = workmatrix;
int pivotRow = i;
for (int row = i + 1; row < order+1; ++row) {
if (fabs(workmatrix) > fabs(pivotElement)) {
pivotElement = workmatrix;
pivotRow = row;
}
}
// important:
// there is no solution, if the pivot element is zero!
if (pivotElement == 0.0) {
return 0.0;
}
// if the found pivot element is not in the actual row,
// we need to swap the rows..
if (pivotRow != i) {
// swapping
for (int k = 0; k < order+1; ++k) {
double t = workmatrix;
workmatrix = workmatrix;
workmatrix = t;
}
// change the sign of the determinant after swapping a line!
det *= -1.0;
}
// The determinant is calculated by a multiplication of the diagonal
// values of the upper triangle matrix part.
// in this case, the pivot element of each column
det *= pivotElement;
// Gaussian triangular steps to reach a row echelon form of the matrix
for (int row = i + 1; row < order+1; ++row) {
for (int col = i + 1; col < order+1; ++col) {
workmatrix -= workmatrix * workmatrix / pivotElement;
}
}
}
}
return det;
}
//==============================================================
// learn:
// Adding a new pair of x and y values and
// solving the polynomial regression model with Cramers rule
void curve_fit::learn(double x, double y) {
++N;
// find min and max values for x
if(N == 1) {
max_x_ = x;
min_x_ = x;
} else {
if(x < min_x_)
min_x_ = x;
if(x > max_x_)
max_x_ = x;
}
// order is -1 when the matrices are not allocated (undefines)
if(order > -1) {
// Example for matrix values
// A second-order polynomial (or: parabolic curve)
// | |
// | N SUM(xi) SUM(xi^2) | M M M
//M = |SUM(xi) SUM(xi^2)SUM(xi^3) |= M M M
// | SUM(xi^2)SUM(xi^3)SUM(xi^4) | M M M
// | |
// Matrix notation:
// aij for each element of the matrix:
// | | | | | |
// | a_i1_j1 a_i1_j2 a_i1_j3 | | a11 a12 a13 | | M M M |
// M = | a_i2_j1 a_i2_j2 a_i2_j3 | = | a21 a22 a23 | = | M M M |
// | a_i3_j1 a_i3_j2 a_i3_j3 | | a31 a32 a33 | | M M M |
// | | | | | |
// to calculate the matrix index based on i and j, the function Mindex(i, j)
// can be used..
//
// First, add the new x values to the existing values of the matrix M[]:
for(int j = 0; j < order+1; j++) {
for(int i = 0; i < order+1; i++) {
// fill the first column with pow() for each element
if(j==0) {
M += pow(x, i);
} else {
if(i < order){
// fill with already calculated values from
// the previous column to save calulation power
M = M;
} else {
// calculate only last element of the column with pow()
M += pow(x, i+j);
}
}
}
}
// fill the first element with the number of samples
M = N;
// Secondly add the new y value to the vector b[]
// | |
// | SUM(yi) | b
//b = |SUM(xi*yi)| = b
// | SUM(xi^2*yi) | b
// | |
for(int n = 0; n < order+1; n++) {
b += (pow(x, n)*y);
}
// Thirdly calculate the coefficients by dividing the determinants
//a0=det(M0)/det(M)
//a1=det(M1)/det(M)
//...
//an=det(Mn)/det(M)
double det_M = determinant(M);
// create a work-matrix to calculate M0, M1, ... , Mn
// according to Cramer's rule
double Mn[(order+1)*(order+1)];
for(int n = 0; n < order+1; n++) {
for(int j = 0; j < order+1; j++) {
for(int i = 0; i < order+1; i++) {
if(j == n)
Mn = b;
else
Mn = M;
}
}
// calculate the coefficients by dividing the determinants
double det_Mn = determinant(Mn);
a = det_Mn / det_M;
}
}
}
//==============================================================
// predict:
// returning the predicted y values of a given x values
// based on the calculated (learned) polynomial regression model
double curve_fit::predict(double x) {
double y = 0.0;
for(int i = 0; i <= order; i++)
y = y + pow(x, i) * a;
return y;
}
//==============================================================
// predict:
// returning the predicted y values of a given x values
// based on the calculated (learned) polynomial regression model
// if x is outside the learned range, y is replaced with outside_value
double curve_fit::predict(double x, double outside_value) {
double y =0.0;
if(x >max_x_ || x < min_x_){
y = outside_value;
} else {
y = predict(x);
}
return y;
}
//==============================================================
// Fills the given field with the current coefficients.
// the values field need to have the right size!
void curve_fit::get_coefficients(double values[]) {
// order 0 -->y = a
// order 1 -->y = a*x + a
// order 2 -->y = a*x^2 + a*x + a
for(int i = 0; i <= order; ++i)
values = a;
}
//==============================================================
// Return the formula as String with all coefficients.
// The number of decimal places can be set the optional parameter.
String curve_fit::get_formula(uint8_t decimals) {
// order n -->y = a*x^n + ... + a*x + a
String formula = "("+String(N)+") y= ";
for(int i = order; i >= 0; --i) {
// at a '+' if the value is positiv
if(a > 0 && i < order)
formula += "+";
if(i > 1) {
formula += String(a,decimals)+"x^"+i+" ";
} else {
if(i== 1)
formula += String(a,decimals)+"x ";
else
formula += String(a,decimals);
}
}
return formula;
}
//==============================================================
// return the actual order.
uint32_t curve_fit::get_order() {
return order;
}
//==============================================================
// return the matrix index based on i and j
// M = M
// M = M
// M = M[(i*(order+1))+j]
uint32_t curve_fit::Mindex(uint32_t i, uint32_t j) {
return (i*(order+1))+j;
}
//==============================================================
// return an estimation of the max y value over the existing x
// range (min_x .. max_x)
// The optional parameter "steps" can be used to define the
// number of steps between min_x and max_x (default = 100)
double curve_fit::estimate_max_y(uint32_t steps) {
if(steps == 0)
steps = 1;
double stepwidth = (max_x_ - min_x_) / steps;
double max_y_ = predict(min_x_);
for(int i = 1; i <= steps; ++i) {
if(predict(min_x_ + (i*stepwidth)) > max_y_)
max_y_ = predict(min_x_ + (i*stepwidth));
}
return max_y_;
}
//==============================================================
// return an estimation of the min y value over the existing x
// range (min_x .. max_x)
// The optional parameter "steps" can be used to define the
// number of steps between min_x and max_x (default = 100)
double curve_fit::estimate_min_y(uint32_t steps) {
if(steps == 0)
steps = 1;
double stepwidth = (max_x_ - min_x_) / steps;
double min_y_ = predict(min_x_);
for(int i = 1; i <= steps; ++i) {
if(predict(min_x_ + (i*stepwidth)) < min_y_)
min_y_ = predict(min_x_ + (i*stepwidth));
}
return min_y_;
}
源文件地址https://github.com/electricidea/M5Stack-Hotel-room-finder
页:
[1]