问MATLAB用hist函数计算直方图的算法是什么？
EN

Stack Overflow用户

提问于 2022-11-20 02:18:15

回答 1查看 35关注 0票数 -1

我正在把一些旧的MATLAB代码翻译成C++。我已经注意到，我的自定义函数计算直方图，应该是等价于MATLAB [counts,centers]= hist(___)给出不同的结果。我在实现中找不到漏洞，所以我使用MATLAB从MATLAB代码中生成C++函数，并将其与我的C++代码进行比较。下面是我用来生成C++代码的一个简单的MATLAB函数：

function [counts, centers] = my_hist(values, bins)
    [counts, centers] = hist(values, bins);
    disp(centers);
    disp(counts);
end

以及调用它的脚本，所以MATLAB可以定义输入：

values = rand(1,1000);
bins = linspace(0.05, 0.95, 10);

[counts, centers] = my_hist(values, bins);

基于上述，Coder将生成以下函数：

//
// File: my_hist.cpp
//
// MATLAB Coder version            : 5.3
// C/C++ source code generated on  : 17-Nov-2022 15:46:17
//

// Include Files
#include "my_hist.h"
#include "rt_nonfinite.h"
#include <algorithm>
#include <cmath>
#include <cstring>
#include <math.h>

// Function Definitions
//
// MY_HIST Summary of this function goes here
//    Detailed explanation goes here
//
// Arguments    : const double values[1000]
//                const double bins[10]
//                double counts[10]
//                double centers[10]
// Return Type  : void
//
void my_hist(const double values[1000], const double bins[10],
             double counts[10], double centers[10])
{
  double edges[11];
  double nn[11];
  double absx;
  int k;
  int low_i;
  std::copy(&bins[0], &bins[10], &centers[0]);
  for (k = 0; k < 9; k++) {
    absx = bins[k];
    edges[k + 1] = absx + (bins[k + 1] - absx) / 2.0;
  }
  edges[0] = rtMinusInf;
  edges[10] = rtInf;
  for (k = 0; k < 9; k++) {
    double absx_tmp;
    absx_tmp = edges[k + 1];
    absx = std::abs(absx_tmp);
    if ((!std::isinf(absx)) && (!std::isnan(absx))) {
      if (absx <= 2.2250738585072014E-308) {
        absx = 4.94065645841247E-324;
      } else {
        frexp(absx, &low_i);
        absx = std::ldexp(1.0, low_i - 53);
      }
    } else {
      absx = rtNaN;
    }
    edges[k + 1] = absx_tmp + absx;
  }
  std::memset(&nn[0], 0, 11U * sizeof(double));
  low_i = 1;
  int exitg1;
  do {
    exitg1 = 0;
    if (low_i + 1 < 12) {
      if (!(edges[low_i] >= edges[low_i - 1])) {
        for (low_i = 0; low_i < 11; low_i++) {
          nn[low_i] = rtNaN;
        }
        exitg1 = 1;
      } else {
        low_i++;
      }
    } else {
      for (k = 0; k < 1000; k++) {
        low_i = 0;
        absx = values[k];
        if (!std::isnan(absx)) {
          if ((absx >= edges[0]) && (absx < edges[10])) {
            int high_i;
            int low_ip1;
            low_i = 1;
            low_ip1 = 2;
            high_i = 11;
            while (high_i > low_ip1) {
              int mid_i;
              mid_i = (low_i + high_i) >> 1;
              if (values[k] >= edges[mid_i - 1]) {
                low_i = mid_i;
                low_ip1 = mid_i + 1;
              } else {
                high_i = mid_i;
              }
            }
          }
          if (values[k] == edges[10]) {
            low_i = 11;
          }
        }
        if (low_i > 0) {
          nn[low_i - 1]++;
        }
      }
      exitg1 = 1;
    }
  } while (exitg1 == 0);
  std::copy(&nn[0], &nn[10], &counts[0]);
  counts[9] += nn[10];
}

//
// File trailer for my_hist.cpp
//
// [EOF]
//

我不明白在这段代码中发生了什么以及为什么要这样做：

  for (k = 0; k < 9; k++) {
    double absx_tmp;
    absx_tmp = edges[k + 1];
    absx = std::abs(absx_tmp);
    if ((!std::isinf(absx)) && (!std::isnan(absx))) {
      if (absx <= 2.2250738585072014E-308) {
        absx = 4.94065645841247E-324;
      } else {
        frexp(absx, &low_i);
        absx = std::ldexp(1.0, low_i - 53);
      }
    } else {
      absx = rtNaN;
    }
    edges[k + 1] = absx_tmp + absx;
  }