【 資料工程與訊號分析】Hampel 濾波器以及中位數絕對離差MAD檢測離群值

#include <pybind11/pybind11.h>

#include <pybind11/stl.h>

#include <algorithm>

#include <numeric>

#include <vector>

#include <cmath>

#include <iostream>

template <typename T = double, typename C>

inline const T median(const C &container) {

 // find median using std::nth_element

 // edit from https://stackoverflow.com/a/39487448

 std::vector<T> arr(std::begin(container), std::end(container));

 size_t middle_ix = arr.size() / 2;

 std::nth_element(arr.begin(), arr.begin() + middle_ix, arr.end());

 if (arr.size() % 2) {

   return arr[middle_ix];

 } else {

   auto iter_max = std::max_element(arr.begin(), arr.begin() + middle_ix);

   return (*iter_max + arr[middle_ix]) / 2.0;

 }

}

inline const double absolute_deviation(double a, double b) {

 return std::fabs(a - b);

}

double hampel_window(const std::vector<double> &sampled_raw_data,

                    int window_size = 5,

                    double n_times_of_sigma_threshold = 1) {

 const unsigned int window_center_ix = window_size / 2;

 const double factor_k = 1.482602218505602;

 double output = 0;

 double median_value = median<double>(sampled_raw_data);

 // compute absolute deviations

 // Xt_ad = |Xt - median (X_window) |

 std::vector<double> absolute_deviations;

 for (auto &value : sampled_raw_data) {

   double AD = absolute_deviation(value, median_value);

   absolute_deviations.push_back(AD);

 }

 // MAD =  median (Xt_ads)

 double MAD = median<double>(absolute_deviations);

 // estimate sigma with MAD

 double sigma_estimate = MAD * factor_k;

 double target_raw_data = (sampled_raw_data)[window_center_ix];

 double target_AD = absolute_deviation(target_raw_data, median_value);

 double threshold = n_times_of_sigma_threshold * sigma_estimate;

 // replace raw data by median or not

 if (target_AD >= threshold) {

   output = median_value;

 } else {

   output = target_raw_data;

 }

 return output;

}

std::vector<double> filter(std::vector<double> &raw_data, int window_size = 5) {

 std::vector<double> raw_data_in_window, data_output;

 for (int i = 0; i < window_size / 2; i++) {

   raw_data_in_window.push_back(0);

 }

 int pos = 0;  // position

 while (pos < raw_data.size() - window_size / 2) {

   // add next data into window at last position

   raw_data_in_window.push_back(raw_data[pos]);

   // if raw data number equal to window size

   if (raw_data_in_window.size() == window_size) {

     // filter data in window

     double filtered = hampel_window(raw_data_in_window, window_size);

     data_output.push_back(filtered);

     // remove first element of raw_data_in_window

     raw_data_in_window.erase(raw_data_in_window.begin());

   }

   pos++;  // move to next data

 }

 return data_output;

}

PYBIND11_MODULE(hampel_filter, m) {

m.def("filter", &filter, "hampel filter");

}

cmake_minimum_required(VERSION 3.18)

set(NAME hampel_filter)

project(${NAME})

set(CMAKE_CXX_FLAGS "-std=c++11" )

set(CMAKE_BUILD_TYPE "Release")

find_package(PythonLibs 3.9.5 REQUIRED)

find_package(PythonInterp 3.9.5 REQUIRED)

set(pybind11_DIR /media/whuang022/DATA/anaconda3/lib/python3.7/site-packages/pybind11/share/cmake/pybind11)

find_package(pybind11 REQUIRED CONFIG )

pybind11_add_module(module ${NAME}.cpp)

set_target_properties(module PROPERTIES OUTPUT_NAME ${NAME})

set_target_properties(module PROPERTIES SUFFIX ".so")

mkdir build

cd build

cmake ..

make

import hampel_filter

import matplotlib.pyplot as plt

input = [10, 2, 9, 4, 2, 6, 20, 11, 3, 10, 11, 6, 20, 11, 3, 1, 5, 9, 10, 10,

        3, 1, 5, 9, 10, 22, 30, 2, 30, 22, 11, 30, 11, 10, 15, 13, 1, 3, 20, 1, 2]

out = hampel_filter.filter(input, 5)

print(out)

plt.plot(input)

plt.plot(out)

plt.show()