【搬运】基于WiFi RSSI定位的酒店房间查找-Arduino中文社区 - Powered by Discuz!

Arduino中文社区

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 2818|回复: 0

【搬运】基于WiFi RSSI定位的酒店房间查找

[复制链接]
发表于 2020-6-8 10:55 | 显示全部楼层 |阅读模式
本帖最后由 vany5921 于 2020-6-8 11:59 编辑

      是否可以使用ESP32或ESP8266设备以足够的精度确定酒店房间的位置?答案很简单:是的,如果酒店周围有足够的WiFI接入点,有足够的时间在地板上走来走去几个小时,而且不会被他人起疑心,就可以做到。
       背景:在酒店房间里闲逛会很无聊。电视节目没有什么可看的,在网络视频上所有的东西都已经被看了。好吧,然后我发现我可以在房间里检测到14个无线接入点。哇,太多了。在我楼层的楼梯上,我检测到16个,在走廊上,我能够检测到11到15个WiFi接入点之间的东西。这真的是太棒了,这个程序就是在这么无聊的情况下写出来的。
floor_rvJ8iolsyh.jpg
你知道吗:酒店的走廊又长又单调。如果你经常更换酒店,你往往忘记房间在哪个位置。所以问题来了:我的M5Stack和大量的WiFi接入点能帮助我快速方便地找到我的房间吗?下图是我房间的位置
map_MLdrncb7LL.jpg
答案取决于你如何定义“快速”和“轻松”,但基本上是这样的:

1.绘制地图
作为第一步,系统必须了解您房间在走廊上的位置。要执行此操作,请从房间门前开始,向左(正)或向右(负)移动。在每一步之后,让M5Stack扫描WiFi接入点,让他了解位置和WiFi信号强度之间的关系。为了尽可能宽地表示走廊,向左和向右移动是很有用的。最好在走廊上走来走去两次,但要提防其他酒店客人:他们宁愿相信你想犯罪,也不不相信你是在做实验。

GUI界面将引导您通过走廊。你只需要一步一步,让设备在每个位置扫描。重要的是你总是迈同样大小的步子。否则结果会不准确!我个人认为,横向搜索更容易,即使在搜索过程中我不得不再次横向搜索以避免测量差异。侧向行走更容易保持相同的步长。记住:你不想在浴室里用别人的牙刷在错误的房间里醒来!所以你最好花点时间做个详细的测量。
2.IILTM:逆向强度查找表
在内部,软件将5阶多项式拟合到学习到的数据。此功能映射走廊中的位置与每个检测到的WiFi接入点的信号强度(RSSI)之间的关系。
usable_aps_Ft1eY2DW7x.jpeg
然而,对于走廊沿线的定位,则需要公式的另一个方向:位置必须根据WiFi接入点的信号强度确定。为了能够根据信号强度计算位置,将在走廊长度上创建反向强度查找表地图:
iiltm_JaYQMhZhCr.jpeg
别担心,只要你在走廊里闲逛完,软件就会自动完成这一切(详情如下)。
3.出去喝一杯
现在什么都不会发生!你可以关掉M5Stack,去最近的酒吧转转。这个系统知道你的床在哪里,并帮助你在酒店无尽的走廊中找到自己的房间。
4.让设备带你回房间
好吧,你应该已经知道你的房间在几楼了。否则,你将不得不经历所有这些,你将失去从你必要的就寝时间。最好的办法是把一只旧袜子放在楼梯上标出正确的楼层(或任何不显眼的地方)。在走廊里,你可以问M5Stack你离房间有多远。设备扫描WiFi接入点,并通过反向强度查找表映射(IILTM)。对于每个条目,计算期望信号强度和测量信号强度之间差的平方和。此曲线的最小值表示您最有可能位于走廊中的位置。但是,如果显示的是一个很大的正数,那么您仍然离房间很远。

position_11_zz79oAQRXz.JPG
如果显示器显示的是负数,则表示您已经走过头了;如果显示器显示的是零,则表示您应该站在房间门的正前方
工作原理:
传统的室内定位方法是基于训练调节的方法。但要做到这一点,你必须知道发射器的确切位置。在酒店,很难精确测量所有WiFi接入点的位置。如果侧身在走廊游荡已经让人可疑,那么想想如果你把天花板的盖子取下来开始寻找电子设备会是什么样子?此外,为了得到距离和信号强度之间的可计算关系,信号必须是直线连接的。这与现代建筑相矛盾,在现代建筑中,这些是不可见的。
另一种常用的方法是使用指纹数据库根据信号强度的模式从信号强度推断位置。我的方法类似,但不同:
在现代酒店里,有许多WiFi接入点,为每一位客人提供互联网接入,无论他们身在何处。在我的酒店走廊里,我可以检测到许多接入点,当我沿着走廊走下去时,到接入点的距离会发生变化,检测到的信号强度也随之改变。换言之:走廊沿线的每个位置或多或少都有一个无障碍接入点信号强度的独特组合。然而,正是在WiFi信号的性质中,信号强度(RSSI)是强噪声的,这使得精确定位变得困难。此外,信号强度模式的基于快照的指纹也会受到信号噪声的影响。现在让数学帮我们实现。

信号形状
一般来说,这听起来很简单:离接入点越远,信号强度越低。所以我应该能从信号强度推断出距离,对吧?不幸的是,它有点复杂,信号的噪声并不能使它更容易。让我们看看信号强度和距离之间的相关性,以了解其复杂性:
当一个发射机在我的视线范围内时,距离和信号强度之间有明显的相关性。这种相关性不是线性的,但仍然是一条简单的曲线:
simple_distance_QdkapDRB8I.jpeg

这看起来很有希望,因为我可以清楚地找到给定RSSI值的距离。然而,当我通过发射器到达我的房间时,就变得更加困难了。信号强度和我在走廊里的位置不再有明显的关系。我可能同时在走廊的两个地方
passing_access_point_oEzV6GJUtM.jpeg
这两种说法挺有趣,但是这确实是个令人头痛的问题
如果发射机位于不同楼层或柱后或其他一些结构部件后面,则沿走廊的信号强度可能会变得更加复杂。这些东西不再是一个简单的函数,而且使你更难确定自己的位置。但是,如果走廊区域中有许多接入点,这些接入点希望随机分布,那么可以从大量数据中计算出更精确的位置信号。因此,我们的想法是,不要指望一个接入点,而是指望图片中许多的接入点。类似于以下示例,其中-5的距离由“红色”访问点的-72dBm值和“蓝色”访问点的-45dBm值明确给出:
two_access_points_02_W4Zd3Gi7Qe.jpeg
但要做到这一点,对于酒店楼层沿线的接入点数量,需要对每个接入点的测量数据进行功能表示。
高阶多项式可用于确定位置与测量信号强度值之间的函数关系。像我的M5Stack中的ESP32这样的小型微控制器不太擅长执行高精度的数学计算,但我已经尝试创建一个库,它可以以可接受的精度拟合6次多项式。在高于6或7时,由于浮点运算的限制,结果变得越来越不准确。好吧,很公平。我最喜欢的酒店的测试数据表明,5阶多项式是足够的。

信号质量
曲线拟合库的拟合类允许您在运行时添加新的数据对。这允许多项式表示随时间“增长”,而无需将所有数据临时存储在RAM内存中。数据点越多,拟合越好。如果数据点太少,结果可能会变得混乱。因此,在接下来的步骤中,在走廊中的数据收集完成后,将忽略收集到的数据点少于6个的所有访问点的数据记录。此外,沿走廊行走时,所有信号电平振幅小于15dBm的接入点也将被忽略。
ap1_LZG9i5MSka.jpeg
软件还记住了走廊中的最大位置。由此,现在创建了一个位置增量小于所做步骤的表。使用这些值可以创建具有所有可用访问点的道路的反向强度查找表地图。由于该图基于走廊沿线信号强度的函数表示,因此有效地滤除了信号强度的典型噪声。
iiltm_02_r6bOzG2JpV.jpeg
此映射已保存,稍后将用于本地化。地图还可以传输到其他M5Stack设备,这样我的朋友就可以找到我的房间,例如,不用告诉他们我的房间号。好吧,好吧,也许告诉他们房间号会更容易,但没那么有趣!
最小二乘位置检测法

为了计算到我房间的距离,这个设备扫描了我站着的地方的无障碍接入点的信号强度。现在如果我站在相应的位置,对于反向强度查找表映射中的每个条目来说,它将检查访问点必须具有的信号强度。计算与测量信号强度之差并将其平方。对所有访问点都执行同样的操作,并对值进行求和。
sq_gif_owIwU7Qwd7.gif.2020-06-08 10_50_10.gif
如果我离强度图上的某个点很远,那么差的平方和会给出一个很高的值。但对于我所处的位置,总和很小,因为测量信号强度和理想信号强度之间的差异很小。因此,强度图上这个计算的最小值表明我很可能在哪里。
pos-1_3_7bt2paz7Sk.jpeg
就这样。仅此而已:一个带有WiFi芯片的设备,一点数学知识,一个小GUI,几个小时就可以在大厅里上下走动。没别的了!在我最喜欢的酒店里,它工作得很好,但在一个只有酒吧才有WiFi的老棚子里,最好拍张房间门的照片,或者只是看看房间钥匙就知道号码。或者有人在接待处询问,并在一张纸上写下房间号。但是在一个有很多接入点的现代化酒店里,这个系统非常容易使用,即使这完全没有意义。
[mw_shl_code=arduino,true]/**************************************************************************
* 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[max_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[100];  
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("[OK] %i Networks found\n", n_WiFi_networks);
            else
                M5.Lcd.println("[ERR] 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("[OK] data deleted");
            else
              M5.Lcd.println("\n\n[ERR] 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("[ERR] 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("[OK] %i Networks found\n", n_WiFi_networks);
            else
              M5.Lcd.println("\n\n[ERR] 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("[ERR] 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[ERR] Failed to open file");
        return;
    } else {
        if(!file.println(message))
            M5.Lcd.println("\n\n[ERR] 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[0] ... newx_array[n_newx-1]
      // BSSIDLT[0] ... BSSIDLT[n_usable_APs-1]
      // IILTM[0] ... 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("[ERR] 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[line_count++] = 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[line_count++], 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[AP_index].learn(0.0, RSSI);
            fits[AP_index].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[x] = 0.0;
    for(int AP_index = 0; AP_index < n_usable_APs; ++AP_index){
      if(fits[AP_index].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[AP_index].predict(0.0) - IILTM[(AP_index*n_newx)+x]);
          square_sum_array[x] += (RSSI_diff*RSSI_diff);
        }
      }
    }
    // Finally find the minimum of the square sums:
    double best_pos = newx_array[0];
    double min_sum = square_sum_array[0];
    for(int x = 0; x < n_newx; ++x){
      if(square_sum_array[x] < min_sum){
        min_sum = square_sum_array[x];
        best_pos = newx_array[x];
      }
    }
    // 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[0] || best_pos == newx_array[n_newx-1])
      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("[ERR] 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[AP_count], 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[x], -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[x]);
        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[0] ... newx_array[n_newx-1]
        // BSSIDLT[0] ... BSSIDLT[n_usable_APs-1]
        // IILTM[0] ... 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[x]);
        }
        // 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;
}[/mw_shl_code]

[mw_shl_code=cpp,true]

/***************************************************
*
* 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[n]*x^n + ... + a[1]*x + a[0]
        double *a;
        // Number of learned x, y pairs
        uint32_t N = 0;
        double max_x_, min_x_;
};[/mw_shl_code]

[mw_shl_code=cpp,true]/**************************************************************************
* 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 rule  for 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[2];
*          fit_1.get_coefficients(coeff_values);
*          Serial.printf("y= a1*x + a0\n\na1= %.2f\na0= %.2f",
*                        coeff_values[1], coeff_values[0]);
*
* 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[Mindex(i,j)] = 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[0];
        if(order == 1)
            // a 2x2 matrix can be calculated directly:
            //       |         |                       |            |
            //       | a11 a12 |                       | M[0]  M[1] |
            //   det | a21 a22 | = a11*a22 – a12*a21 = | M[2]  M[3] | = M[0]*M[3] - M[1]*M[2]
            //       |         |                       |            |
            det = ((mainmatrix[0] * mainmatrix[3]) - (mainmatrix[1] * mainmatrix[2]));
    } 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[Mindex(i,j)] = mainmatrix[Mindex(i,j)];
        // 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[Mindex(i,i)];
            int pivotRow = i;
            for (int row = i + 1; row < order+1; ++row) {
                if (fabs(workmatrix[Mindex(row,i)]) > fabs(pivotElement)) {
                    pivotElement = workmatrix[Mindex(row,i)];
                    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[Mindex(i,k)];
                    workmatrix[Mindex(i,k)] = workmatrix[Mindex(pivotRow,k)];
                    workmatrix[Mindex(pivotRow,k)] = 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[Mindex(row,col)] -= workmatrix[Mindex(row,i)] * workmatrix[Mindex(i,col)] / 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[0] M[1] M[2]
        //  M = |  SUM(xi)   SUM(xi^2)  SUM(xi^3) |  = M[3] M[4] M[5]
        //      | SUM(xi^2)  SUM(xi^3)  SUM(xi^4) |    M[6] M[7] M[8]
        //      |                                 |  
        // Matrix notation:
        // aij for each element of the matrix:
        //       |                         |   |             |   |                |
        //       | a_i1_j1 a_i1_j2 a_i1_j3 |   | a11 a12 a13 |   | M[0] M[1] M[2] |
        //   M = | a_i2_j1 a_i2_j2 a_i2_j3 | = | a21 a22 a23 | = | M[3] M[4] M[5] |
        //       | a_i3_j1 a_i3_j2 a_i3_j3 |   | a31 a32 a33 |   | M[6] M[7] M[8] |
        //       |                         |   |             |   |                |  
        // 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[Mindex(i, j)] += pow(x, i);
                } else {
                    if(i < order){
                        // fill with already calculated values from
                        // the previous column to save calulation power
                        M[Mindex(i, j)] = M[Mindex(i+1, j-1)];
                    } else {
                        // calculate only last element of the column with pow()
                        M[Mindex(i, j)] += pow(x, i+j);
                    }
                }
            }
        }
        // fill the first element with the number of samples
        M[0] = N;

        // Secondly add the new y value to the vector b[]
        //      |              |
        //      |    SUM(yi)   |   b[0]
        //  b = |  SUM(xi*yi)  | = b[1]
        //      | SUM(xi^2*yi) |   b[2]
        //      |              |
        for(int n = 0; n < order+1; n++) {
            b[n] += (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[Mindex(i, j)] = b;
                    else
                        Mn[Mindex(i, j)] = M[Mindex(i, j)];
                }
            }
            // calculate the coefficients by dividing the determinants
            double det_Mn = determinant(Mn);
            a[n] = 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[0]
    // order 1 -->  y = a[1]*x + a[0]
    // order 2 -->  y = a[2]*x^2 + a[1]*x + a[0]
    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[n]*x^n + ... + a[1]*x + a[0]
    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[0][0] = M[0]
// M[0][3] = M[3]
// M[j] = 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_;
}[/mw_shl_code]

源文件地址https://github.com/electricidea/M5Stack-Hotel-room-finder

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|Archiver|手机版|Arduino中文社区

GMT+8, 2024-12-1 01:11 , Processed in 0.084903 second(s), 19 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表