From b4683be1f8a2fe6b691860017b86b0f1640bf330 Mon Sep 17 00:00:00 2001 From: beve Date: Fri, 12 May 2023 13:47:22 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8F=91=E5=B8=83=E7=AC=AC=E4=B8=80=E7=89=88?= =?UTF-8?q?=E6=9C=ACSavitzky-Golay=E5=AE=9E=E7=8E=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .clang-format | 136 +++++++++ .gitignore | 3 + CMakeLists.txt | 33 ++ README.en.md | 61 ++-- README.md | 62 ++-- include/filters.hpp | 323 ++++++++++++++++++++ include/polynomial.hpp | 97 ++++++ include/savitzky_golay.hpp | 340 +++++++++++++++++++++ tests.cpp | 602 +++++++++++++++++++++++++++++++++++++ 9 files changed, 1605 insertions(+), 52 deletions(-) create mode 100755 .clang-format create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 include/filters.hpp create mode 100644 include/polynomial.hpp create mode 100644 include/savitzky_golay.hpp create mode 100644 tests.cpp diff --git a/.clang-format b/.clang-format new file mode 100755 index 0000000..b886e6d --- /dev/null +++ b/.clang-format @@ -0,0 +1,136 @@ +Language: Cpp +# BasedOnStyle: LLVM +AccessModifierOffset: -4 +AlignAfterOpenBracket: Align +AlignConsecutiveMacros: true +AlignConsecutiveAssignments: true +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Right +AlignOperands: true +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: false +AllowAllConstructorInitializersOnNextLine: false +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: Empty +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Never +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: Yes +BinPackArguments: false +BinPackParameters: false +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Always + AfterEnum: false + AfterFunction: true + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: true + AfterUnion: true + AfterExternBlock: true + BeforeCatch: false + BeforeElse: true + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Custom +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: true +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 100 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DeriveLineEnding: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: true +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + - Regex: '.*' + Priority: 1 + SortPriority: 0 +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentCaseLabels: false +IndentGotoLabels: true +IndentPPDirectives: None +IndentWidth: 4 +IndentWrappedFunctionNames: false +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 110 +PointerAlignment: Right +ReflowComments: true +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +Standard: Latest +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 4 +UseCRLF: false +UseTab: Never +... + diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ed38492 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +build/* +.vscode + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..cb34c9f --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required(VERSION 2.8.12) + +PROJECT(savgol) + +SET(CMAKE_BUILD_TYPE Debug) + +find_package(PkgConfig) +pkg_check_modules( EIGEN3 REQUIRED eigen3 ) +include_directories( ${EIGEN3_INCLUDE_DIRS} ) + +add_definitions(${EIGEN_DEFINITIONS}) + +OPTION(ENABLE_CXX17 "Enable C++17 support" ON) + +IF(ENABLE_CXX17) + INCLUDE(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG("-std=c++17" COMPILER_SUPPORTS_CXX17) + CHECK_CXX_COMPILER_FLAG("-std=c++0x" COMPILER_SUPPORTS_CXX0X) + IF(COMPILER_SUPPORTS_CXX17) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17") + ELSEIF(COMPILER_SUPPORTS_CXX0X) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x") + ELSE() + MESSAGE(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++17 support. Please use a different C++ compiler.") + ENDIF() +ENDIF(ENABLE_CXX17) + +INCLUDE_DIRECTORIES(EIGEN3_INCLUDE_DIRS + AND "${PROJECT_SOURCE_DIR}/include") + + +#Declare the executable target built from your sources +add_executable(savgol tests.cpp) diff --git a/README.en.md b/README.en.md index b73d073..039ef57 100644 --- a/README.en.md +++ b/README.en.md @@ -1,36 +1,47 @@ -# Savitzky-Golay +## Savitzky-Golay Filter in C++ -#### Description -{**When you're done, you can delete the content in this README and update the file with details for others getting started with your repository**} +Author: [Dong anakin](henudwj@126.com) -#### Software Architecture -Software architecture description -#### Installation +### Introduction -1. xxxx -2. xxxx -3. xxxx +Implementation of the Savitzky-Golay Filter. Initial testing of this code was on a Ubuntu 18.04 but will work on any other Linux/Windows/Mac OS machine with little effort. -#### Instructions +### Dependencies -1. xxxx -2. xxxx -3. xxxx ++ [Eigen3 Library](http://eigen.tuxfamily.org/index.php?title=Main_Page) for the linear algebra of the matrices, vectors and related algorithms. -#### Contribution +### Example -1. Fork the repository -2. Create Feat_xxx branch -3. Commit your code -4. Create Pull Request +* `./savgol` +### Compilation -#### Gitee Feature +There is a `CMakeLists.txt` file in the project root folder. From the project root directory: -1. You can use Readme\_XXX.md to support different languages, such as Readme\_en.md, Readme\_zh.md -2. Gitee blog [blog.gitee.com](https://blog.gitee.com) -3. Explore open source project [https://gitee.com/explore](https://gitee.com/explore) -4. The most valuable open source project [GVP](https://gitee.com/gvp) -5. The manual of Gitee [https://gitee.com/help](https://gitee.com/help) -6. The most popular members [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) +1. Create a build directory: `mkdir build && cd build` +2. Compile the `cpp` code: `cmake ../` +3. Build your executable: `make` +4. Run the executable: `./savgol` + + +### Citation + +If you have used `Savitzky-Golay` in your work, please cite it. + +```tex +@misc{Savitzky-Golay, + author = {Dong anakin}, + title = {{Savitzky-Golay Filter in C++}}, + year = {2023}, + howpublished = {\url{ https://gitee.com/beve/savitzky-golay}}, + note = {Accessed April, 2023} +} +``` +#### Issues + +If you have issues running the files, please use the issues tab to open a bug. I will generally respond within a 24-hour period. + + +### Reference + scipy.signal.savgol_filter https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html \ No newline at end of file diff --git a/README.md b/README.md index 068b219..039ef57 100644 --- a/README.md +++ b/README.md @@ -1,39 +1,47 @@ -# Savitzky-Golay +## Savitzky-Golay Filter in C++ -#### 介绍 -{**以下是 Gitee 平台说明，您可以替换此简介** -Gitee 是 OSCHINA 推出的基于 Git 的代码托管平台（同时支持 SVN）。专为开发者提供稳定、高效、安全的云端软件开发协作平台 -无论是个人、团队、或是企业，都能够用 Gitee 实现代码托管、项目管理、协作开发。企业项目请看 [https://gitee.com/enterprises](https://gitee.com/enterprises)} +Author: [Dong anakin](henudwj@126.com) -#### 软件架构 -软件架构说明 +### Introduction -#### 安装教程 +Implementation of the Savitzky-Golay Filter. Initial testing of this code was on a Ubuntu 18.04 but will work on any other Linux/Windows/Mac OS machine with little effort. -1. xxxx -2. xxxx -3. xxxx +### Dependencies -#### 使用说明 ++ [Eigen3 Library](http://eigen.tuxfamily.org/index.php?title=Main_Page) for the linear algebra of the matrices, vectors and related algorithms. -1. xxxx -2. xxxx -3. xxxx +### Example -#### 参与贡献 +* `./savgol` -1. Fork 本仓库 -2. 新建 Feat_xxx 分支 -3. 提交代码 -4. 新建 Pull Request +### Compilation +There is a `CMakeLists.txt` file in the project root folder. From the project root directory: -#### 特技 +1. Create a build directory: `mkdir build && cd build` +2. Compile the `cpp` code: `cmake ../` +3. Build your executable: `make` +4. Run the executable: `./savgol` -1. 使用 Readme\_XXX.md 来支持不同的语言，例如 Readme\_en.md, Readme\_zh.md -2. Gitee 官方博客 [blog.gitee.com](https://blog.gitee.com) -3. 你可以 [https://gitee.com/explore](https://gitee.com/explore) 这个地址来了解 Gitee 上的优秀开源项目 -4. [GVP](https://gitee.com/gvp) 全称是 Gitee 最有价值开源项目，是综合评定出的优秀开源项目 -5. Gitee 官方提供的使用手册 [https://gitee.com/help](https://gitee.com/help) -6. Gitee 封面人物是一档用来展示 Gitee 会员风采的栏目 [https://gitee.com/gitee-stars/](https://gitee.com/gitee-stars/) + +### Citation + +If you have used `Savitzky-Golay` in your work, please cite it. + +```tex +@misc{Savitzky-Golay, + author = {Dong anakin}, + title = {{Savitzky-Golay Filter in C++}}, + year = {2023}, + howpublished = {\url{ https://gitee.com/beve/savitzky-golay}}, + note = {Accessed April, 2023} +} +``` +#### Issues + +If you have issues running the files, please use the issues tab to open a bug. I will generally respond within a 24-hour period. + + +### Reference + scipy.signal.savgol_filter https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html \ No newline at end of file diff --git a/include/filters.hpp b/include/filters.hpp new file mode 100644 index 0000000..d73d9d9 --- /dev/null +++ b/include/filters.hpp @@ -0,0 +1,323 @@ +#pragma once + +#include + +namespace scipy { +namespace ndimage { +typedef enum { + NI_EXTEND_FIRST = 0, + NI_EXTEND_NEAREST = 0, + NI_EXTEND_WRAP = 1, + NI_EXTEND_REFLECT = 2, + NI_EXTEND_MIRROR = 3, + NI_EXTEND_CONSTANT = 4, + NI_EXTEND_GRID_WRAP = 5, + NI_EXTEND_GRID_CONSTANT = 6, + NI_EXTEND_LAST = NI_EXTEND_GRID_WRAP, + NI_EXTEND_DEFAULT = NI_EXTEND_MIRROR +} NI_ExtendMode; + +NI_ExtendMode extend_mode_to_code(std::string mode) +{ + // Convert an extension mode to the corresponding integer code. + + if (mode == "nearest") + { + return NI_EXTEND_NEAREST; + } + else if (mode == "wrap") + { + return NI_EXTEND_WRAP; + } + else if (mode == "reflect") + { + return NI_EXTEND_REFLECT; + } + else if (mode == "grid-mirror") + { + return NI_EXTEND_REFLECT; + } + else if (mode == "mirror") + { + return NI_EXTEND_MIRROR; + } + else if (mode == "constant") + { + return NI_EXTEND_CONSTANT; + } + else if (mode == "grid-wrap") + { + return NI_EXTEND_GRID_WRAP; + } + else if (mode == "grid-constant") + { + return NI_EXTEND_GRID_CONSTANT; + } + else + { // raise RuntimeError('boundary mode not supported') + return NI_EXTEND_REFLECT; + } +} +/// @brief Calculate a 1-D correlation along the given axis. +/// The lines of the array along the given axis are correlated with the +/// given weights. +/// @param input +/// @param weights +/// @param mode +/// @param cval +/// @param origin +/// @return +Eigen::RowVectorXd correlate1d(Eigen::RowVectorXd input, + Eigen::RowVectorXd weights, + std::string mode = "reflect", + float cval = 0.0, + int origin = 0) +{ + eigen_assert((weights.size()) && "no filter weights given"); + if (!weights.size()) + { + return {}; + } + int weights_size = weights.size(); + eigen_assert(((-weights_size / 2) <= origin) && (origin <= ((weights_size - 1) / 2)) && + "Invalid origin; origin must satisfy " + "-(len(weights) // 2) <= origin <= " + "(len(weights)-1) // 2"); + if (!(((-weights_size / 2) <= origin) && (origin <= ((weights_size - 1) / 2)))) + { + return {}; + } + int symmetric = 0; + /* test for symmetry or anti-symmetry: */ + const int filter_size = weights.size(); + const int size1 = filter_size / 2; + const int size2 = filter_size - size1 - 1; + double *fw = weights.data(); + if (filter_size & 0x1) + { + symmetric = 1; + for (int ii = 1; ii <= (filter_size / 2); ++ii) + { + if (std::fabs(fw[ii + size1] - fw[size1 - ii]) > __DBL_EPSILON__) + { + symmetric = 0; + break; + } + } + if (symmetric == 0) + { + symmetric = -1; + for (int ii = 1; ii <= (filter_size / 2); ++ii) + { + if (std::fabs(fw[size1 + ii] + fw[size1 - ii]) > __DBL_EPSILON__) + { + symmetric = 0; + break; + } + } + } + } + const int size1_extern = size1 + origin; + const int size2_extern = size2 - origin; + const int line_length = input.size(); + Eigen::RowVectorXd input_extend = + Eigen::RowVectorXd::Zero(line_length + size1_extern + size2_extern); + for (int cc = 0; cc < input.size(); ++cc) + { + input_extend[cc + size1_extern] = input[cc]; + } + if (size1_extern + size2_extern > 0) + { + NI_ExtendMode extend_mode = extend_mode_to_code(mode); + + int size_before = size1_extern; + int size_after = size2_extern; + double *buffer = input_extend.data(); + double *first = buffer + size1_extern; + double *last = first + line_length; + double *src, *dst, val; + + if ((line_length == 1) && (extend_mode == NI_EXTEND_MIRROR)) + { + extend_mode = NI_EXTEND_NEAREST; + } + + switch (extend_mode) + { + /* aaaaaaaa|abcd|dddddddd */ + case NI_EXTEND_NEAREST: + src = first; + dst = buffer; + val = *src; + while (size_before--) + { + *dst++ = val; + } + src = last - 1; + dst = last; + val = *src; + while (size_after--) + { + *dst++ = val; + } + break; + /* abcdabcd|abcd|abcdabcd */ + case NI_EXTEND_WRAP: + case NI_EXTEND_GRID_WRAP: + src = last - 1; + dst = first - 1; + while (size_before--) + { + *dst-- = *src--; + } + src = first; + dst = last; + while (size_after--) + { + *dst++ = *src++; + } + break; + /* abcddcba|abcd|dcbaabcd */ + case NI_EXTEND_REFLECT: + src = first; + dst = first - 1; + while (size_before && src < last) + { + *dst-- = *src++; + --size_before; + } + src = last - 1; + while (size_before--) + { + *dst-- = *src--; + } + src = last - 1; + dst = last; + while (size_after && src >= first) + { + *dst++ = *src--; + --size_after; + } + src = first; + while (size_after--) + { + *dst++ = *src++; + } + break; + /* cbabcdcb|abcd|cbabcdcb */ + case NI_EXTEND_MIRROR: + src = first + 1; + dst = first - 1; + while (size_before && src < last) + { + *dst-- = *src++; + --size_before; + } + src = last - 2; + while (size_before--) + { + *dst-- = *src--; + } + src = last - 2; + dst = last; + while (size_after && src >= first) + { + *dst++ = *src--; + --size_after; + } + src = first + 1; + while (size_after--) + { + *dst++ = *src++; + } + break; + /* kkkkkkkk|abcd]kkkkkkkk */ + case NI_EXTEND_CONSTANT: + case NI_EXTEND_GRID_CONSTANT: + val = cval; + dst = buffer; + while (size_before--) + { + *dst++ = val; + } + dst = last; + while (size_after--) + { + *dst++ = val; + } + break; + default: + break; + } + } + + Eigen::RowVectorXd output = Eigen::RowVectorXd::Zero(line_length); + + /* get lines: */ + double *iline = input_extend.data() + size1; + double *oline = output.data(); + fw = weights.data() + size1; + /* the correlation calculation: */ + if (symmetric > 0) + { + for (int ll = 0; ll < line_length; ll++) + { + oline[ll] = iline[0] * fw[0]; + for (int jj = -size1; jj < 0; jj++) + oline[ll] += (iline[jj] + iline[-jj]) * fw[jj]; + ++iline; + } + } + else if (symmetric < 0) + { + for (int ll = 0; ll < line_length; ll++) + { + oline[ll] = iline[0] * fw[0]; + for (int jj = -size1; jj < 0; jj++) + oline[ll] += (iline[jj] - iline[-jj]) * fw[jj]; + ++iline; + } + } + else + { + for (int ll = 0; ll < line_length; ll++) + { + oline[ll] = iline[size2] * fw[size2]; + for (int jj = -size1; jj < size2; jj++) + oline[ll] += iline[jj] * fw[jj]; + ++iline; + } + } + + return output; +} +/// @brief Calculate a 1-D convolution along the given axis. +/// The lines of the array along the given axis are convolved with the +/// given weights. +/// @param input ndarray +/// 1-D sequence of numbers. +/// @param weights ndarray +/// 1-D sequence of numbers. +/// @param mode +/// @param cval +/// @param origin +/// @return convolve1d : ndarray +/// Convolved array with same shape as input +Eigen::RowVectorXd convolve1d(Eigen::RowVectorXd input, + Eigen::RowVectorXd weights, + std::string mode = "reflect", + float cval = 0.0, + int origin = 0) +{ + // weights = weights[::-1] + weights.reverseInPlace(); + origin = -origin; + if (!(weights.size() & 1)) + { + origin -= 1; + } + return correlate1d(input, weights, mode, cval, origin); +} +} // namespace ndimage + +} // namespace scipy diff --git a/include/polynomial.hpp b/include/polynomial.hpp new file mode 100644 index 0000000..e925d49 --- /dev/null +++ b/include/polynomial.hpp @@ -0,0 +1,97 @@ +#pragma once +#include + +#ifndef savgol_print_matrix +#define savgol_print_matrix(x) \ + do \ + { \ + Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]"); \ + std::cout << #x << " " << __func__ << ":" << __LINE__ << " " << x.rows() << "," \ + << x.cols() << std::endl; \ + std::cout << x.format(HeavyFmt) << std::endl; \ + } while (false) +#endif + +namespace numpy { +/// @brief Least squares polynomial fit. +/// .. note:: +/// This forms part of the old polynomial API. Since version 1.4, the +/// new polynomial API defined in `numpy.polynomial` is preferred. +/// A summary of the differences can be found in the +/// :doc:`transition guide `. +/// +/// Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg` +/// to points `(x, y)`. Returns a vector of coefficients `p` that minimises +/// the squared error in the order `deg`, `deg-1`, ... `0`. +/// +/// The `Polynomial.fit ` class +/// method is recommended for new code as it is more stable numerically. See +/// the documentation of the method for more information. +/// @param x array_like, shape (M,) +/// x-coordinates of the M sample points ``(x[i], y[i])``. +/// @param y array_like, shape (M,) or (M, K) +/// y-coordinates of the sample points. Several data sets of sample +/// points sharing the same x-coordinates can be fitted at once by +/// passing in a 2D-array that contains one dataset per column. +/// @param deg int +/// Degree of the fitting polynomial +/// @return ndarray, shape (deg + 1,) +/// Polynomial coefficients, highest power first. If `y` was 2-D, the +/// coefficients for `k`-th data set are in ``p[:,k]``. +/// +inline Eigen::VectorXd polyfit(Eigen::RowVectorXd x, Eigen::RowVectorXd y, int deg) +{ + if (deg < 0) + { + return {}; + } + int order = deg + 1; + // lhs = vander(x, order) + Eigen::MatrixXd lhs = Eigen::MatrixXd::Zero(x.cols(), order); + for (int r = 0; r < order; r++) + { + lhs.col(order - 1 - r) = x.array().pow(r); + } + Eigen::VectorXd rhs = y; + // Find the least-squares solution of lhs*c = rhs + Eigen::VectorXd coeffs = lhs.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(rhs); + return coeffs; +} +inline Eigen::VectorXd polyder(Eigen::VectorXd p) +{ + if (p.size() <= 1) + { + return Eigen::VectorXd::Zero(1); + } + Eigen::VectorXd res = Eigen::VectorXd::Zero(p.size() - 1); + for (int i = 0; i < res.size(); ++i) + { + res(i) = p(i) * (res.size() - i); + } + return res; +} +inline Eigen::VectorXd polyder(Eigen::VectorXd p, int m) +{ + if (m <= 0) + { + return p; + } + Eigen::VectorXd res = p; + for (int k = 0; k < m; k++) + { + res = polyder(res); + } + return res; +} +inline Eigen::VectorXd polyval(Eigen::VectorXd p, Eigen::VectorXd x) +{ + Eigen::VectorXd y = Eigen::VectorXd::Zero(x.size()); + for (int r = 0; r < x.rows(); r++) + { + double res = 0; + std::for_each(p.cbegin(), p.cend(), [&](auto c) { res = std::fmal(res, x(r), c); }); + y(r) = res; + } + return y; +} +} // namespace numpy diff --git a/include/savitzky_golay.hpp b/include/savitzky_golay.hpp new file mode 100644 index 0000000..078df4c --- /dev/null +++ b/include/savitzky_golay.hpp @@ -0,0 +1,340 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace scipy { +namespace signal { +/// @brief The factorial of a number or array of numbers. +/// The factorial of non-negative integer n is the product of all +/// positive integers less than or equal to n: +/// n! = n * (n - 1) * (n - 2) * ... * 1 +/// @param n +/// @return +inline int factorial(int n) +{ + int fact = 1; + for (int i = 1; i <= n; i++) + { + fact *= i; + } + return fact; +} +#ifndef savgol_print_matrix +#define savgol_print_matrix(x) \ + do \ + { \ + Eigen::IOFormat HeavyFmt(Eigen::FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]"); \ + std::cout << #x << " " << __func__ << ":" << __LINE__ << " " << x.rows() << "," \ + << x.cols() << std::endl; \ + std::cout << x.format(HeavyFmt) << std::endl; \ + } while (false) +#endif +/// @brief Given an N-d array `x` and the specification of a slice of `x` from +/// `window_start` to `window_stop` along `axis`, create an interpolating +/// polynomial of each 1-D slice, and evaluate that polynomial in the slice +/// from `interp_start` to `interp_stop`. Put the result into the +/// corresponding slice of `y`. +/// @param x +/// @param window_start +/// @param window_stop +/// @param interp_start +/// @param interp_stop +/// @param polyorder +/// @param deriv +/// @param delta +/// @param y +inline void fit_edge(Eigen::RowVectorXd x, + int window_start, + int window_stop, + int interp_start, + int interp_stop, + int polyorder, + int deriv, + float delta, + Eigen::RowVectorXd &y) +{ + // Get the edge into a(window_length, -1) array. + // x_edge = axis_slice(x, start = window_start, stop = window_stop, axis = axis); + Eigen::RowVectorXd x_edge = x.segment(window_start, window_stop - window_start); + Eigen::VectorXd xx_edge = x_edge; + bool swapped = false; + // Fit the edges.poly_coeffs has shape(polyorder + 1, -1), + // where '-1' is the same as in xx_edge. + Eigen::VectorXd poly_coeffs = + numpy::polyfit(Eigen::RowVectorXd::LinSpaced(window_stop - window_start, + 0, + window_stop - window_start - 1), + xx_edge, + polyorder); + + if (deriv > 0) + { + poly_coeffs = numpy::polyder(poly_coeffs, deriv); + } + + // # Compute the interpolated values for the edge. + auto i = Eigen::VectorXd::LinSpaced(interp_stop - interp_start, + interp_start - window_start, + interp_stop - window_start - 1); + + Eigen::RowVectorXd values = numpy::polyval(poly_coeffs, i) / (std::pow(delta, deriv)); + + // # Now put the values into the appropriate slice of y. + // # First reshape values to match y. + // shp = list(y.shape) + // shp[0], shp[axis] = shp[axis], shp[0] + // values = values.reshape(interp_stop - interp_start, *shp[1:]) + // if swapped: + // values = values.swapaxes(0, axis) + // # Get a view of the data to be replaced by values. + // auto y_edge = y.segment(interp_start, interp_stop - interp_start - 1); + for (int i = interp_start; i < std::min(int(y.size()), interp_stop); i++) + { + y(i) = values(i - interp_start); + } +} + +/// @brief Use polynomial interpoly_edgeation of x at the low and high ends of the axis +/// to fill in the halflen values in y. +/// This function just calls _fit_edge twice, once for each end of the axis. +/// @param x +/// @param window_length +/// @param polyorder +/// @param deriv +/// @param delta +/// @param y +inline void fit_edges_polyfit(Eigen::RowVectorXd x, + int window_length, + int polyorder, + int deriv, + float delta, + Eigen::RowVectorXd &y) +{ + int halflen = window_length / 2; + fit_edge(x, 0, window_length, 0, halflen, polyorder, deriv, delta, y); + int n = x.size(); + fit_edge(x, n - window_length, n, n - halflen, n, polyorder, deriv, delta, y); +} + +/// @brief Compute the coefficients for a 1-D Savitzky-Golay FIR filter. +/// References +/// https://github.com/scipy/scipy.git scipy/signal/_savitzky_golay.py 0.14.0 +/// +/// A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of Data by +/// Simplified Least Squares Procedures. Analytical Chemistry, 1964, 36 (8), +/// pp 1627-1639. +/// +/// Jianwen Luo, Kui Ying, and Jing Bai. 2005. Savitzky-Golay smoothing and +/// differentiation filter for even number data. Signal Process. +/// 85, 7 (July 2005), 1429-1434. +/// @param window_length +/// The length of the filter window (i.e., the number of coefficients). +/// @param polyorder +/// The order of the polynomial used to fit the samples. +/// `polyorder` must be less than `window_length`. +/// @param deriv +/// The order of the derivative to compute. This must be a +/// nonnegative integer. The default is 0, which means to filter +/// the data without differentiating +/// @param delta +/// The spacing of the samples to which the filter will be applied. +/// This is only used if deriv > 0. +/// @param pos +/// If pos is not nullptr, it specifies evaluation position within the +/// window. The default is the middle of the window. +/// @param use +/// Either 'conv' or 'dot'. This argument chooses the order of the +/// coefficients. The default is 'conv', which means that the +/// coefficients are ordered to be used in a convolution. With +/// use='dot', the order is reversed, so the filter is applied by +/// dotting the coefficients with the data set. +/// @return coeffs +/// 1-D ndarray, The filter coefficients. +inline Eigen::RowVectorXd savgol_coeffs(int window_length, + int polyorder, + int deriv = 0, + float delta = 1.0, + int *pos = nullptr, + std::string use = "conv") +{ + + eigen_assert((window_length > polyorder) && + "transposeInPlace() called on a non-square non-resizable " + "matrix"); + if (polyorder >= window_length) + { + return {}; + } + int halflen = window_length / 2; + int rem = window_length % 2; + float fpos = halflen; + if (pos == nullptr) + { + if (rem == 0) + { + fpos = halflen - 0.5f; + } + else + { + fpos = halflen; + } + } + else + { + fpos = *pos; + } + eigen_assert(((0 <= fpos) && (fpos < window_length)) && "pos must be nonnegative and less than " + "window_length."); + if (!((0 <= fpos) && (fpos < window_length))) + { + return {}; + } + const auto use_array = std::array{"conv", "dot"}; + bool use_in_flag = std::find(use_array.begin(), use_array.end(), use) != use_array.end(); + eigen_assert(use_in_flag && "pos must be nonnegative and less than " + "window_length."); + if (!use_in_flag) + { + return {}; + } + + if (deriv > polyorder) + { + return Eigen::VectorXd::Zero(window_length); + } + + Eigen::RowVectorXd x = + Eigen::VectorXd::LinSpaced(window_length, -fpos, window_length - fpos - 1); + if (use == "conv") + { + x.reverseInPlace(); + } + + Eigen::VectorXi order = Eigen::VectorXi::LinSpaced(polyorder + 1, 0, polyorder); + + // python A = x ** order + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(order.rows(), x.cols()); + for (int r = 0; r < order.rows(); r++) + { + A.row(r) = x.array().pow(order(r, 0)); + } + + // y determines which order derivative is returned. + Eigen::VectorXd y = Eigen::VectorXd::Zero(polyorder + 1); + // The coefficient assigned to y[deriv] scales the result to take into + // account the order of the derivative and the sample spacing. + y(deriv) = factorial(deriv) / (std::pow(delta, deriv)); + // Find the least-squares solution of A*c = y + Eigen::RowVectorXd coeffs = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(y); + return coeffs; +} + +/// @brief Apply a Savitzky-Golay filter to an array. +/// This is a 1-D filter. Input type is RowVector, so remove axis param +/// @param x array_like +/// The data to be filtered. If `x` is not a single or double precision +/// floating point array, it will be converted to type ``numpy.float64`` +/// before filtering. +/// @param window_length int +/// The length of the filter window (i.e., the number of coefficients). +/// If `mode` is 'interp', `window_length` must be less than or equal +/// to the size of `x`. +/// @param polyorder int +/// The order of the polynomial used to fit the samples. +/// `polyorder` must be less than `window_length`. +/// @param deriv int, optional +/// The order of the derivative to compute. This must be a +/// nonnegative integer. The default is 0, which means to filter +/// the data without differentiating. +/// @param delta float, optional +/// The spacing of the samples to which the filter will be applied. +/// This is only used if deriv > 0. Default is 1.0. +/// @param mode str, optional +/// Must be 'mirror', 'constant', 'nearest', 'wrap' or 'interp'. This +/// determines the type of extension to use for the padded signal to +/// which the filter is applied. When `mode` is 'constant', the padding +/// value is given by `cval`. See the Notes for more details on 'mirror', +/// 'constant', 'wrap', and 'nearest'. +/// When the 'interp' mode is selected (the default), no extension +/// is used. Instead, a degree `polyorder` polynomial is fit to the +/// last `window_length` values of the edges, and this polynomial is +/// used to evaluate the last `window_length // 2` output values. +/// @param cval scalar, optional +/// Value to fill past the edges of the input if `mode` is 'constant'. +/// Default is 0.0. +/// @return y : ndarray, same shape as `x` +/// The filtered data. +/// @note Details on the `mode` options: +/// +/// 'mirror': +/// Repeats the values at the edges in reverse order. The value +/// closest to the edge is not included. +/// 'nearest': +/// The extension contains the nearest input value. +/// 'constant': +/// The extension contains the value given by the `cval` argument. +/// 'wrap': +/// The extension contains the values from the other end of the array. +/// +/// For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8], and +/// `window_length` is 7, the following shows the extended data for +/// the various `mode` options (assuming `cval` is 0):: +/// +/// mode | Ext | Input | Ext +/// -----------+---------+------------------------+--------- +/// 'mirror' | 4 3 2 | 1 2 3 4 5 6 7 8 | 7 6 5 +/// 'nearest' | 1 1 1 | 1 2 3 4 5 6 7 8 | 8 8 8 +/// 'constant' | 0 0 0 | 1 2 3 4 5 6 7 8 | 0 0 0 +/// 'wrap' | 6 7 8 | 1 2 3 4 5 6 7 8 | 1 2 3 +/// +// .. versionadded:: 0.14.0 + +inline Eigen::RowVectorXd savgol_filter(Eigen::RowVectorXd x, + int window_length, + int polyorder, + int deriv = 0, + float delta = 1.0, + std::string mode = "interp", + double cval = 0.0) +{ + const auto mode_array = + std::array{"mirror", "constant", "nearest", "interp", "wrap"}; + bool mode_in_flag = std::find(mode_array.begin(), mode_array.end(), mode) != mode_array.end(); + eigen_assert(mode_in_flag && "mode must be 'mirror', 'constant', 'nearest' " + "'wrap' or 'interp'."); + if (!mode_in_flag) + { + return {}; + } + + auto coeffs = savgol_coeffs(window_length, polyorder, deriv, delta); + if (mode == "interp") + { + eigen_assert((window_length <= x.cols()) && + "If mode is 'interp', window_length must be less " + "than or equal to the size of x."); + if (window_length > x.cols()) + { + return {}; + } + // Do not pad. Instead, for the elements within `window_length // 2` + // of the ends of the sequence, use the polynomial that is fitted to + // the last `window_length` elements. + Eigen::RowVectorXd y = scipy::ndimage::convolve1d(x, coeffs, "constant"); + fit_edges_polyfit(x, window_length, polyorder, deriv, delta, y); + return y; + } + else + { + Eigen::RowVectorXd y = scipy::ndimage::convolve1d(x, coeffs, mode, cval); + return y; + } +} + +} // namespace signal +} // namespace scipy diff --git a/tests.cpp b/tests.cpp new file mode 100644 index 0000000..1c5d117 --- /dev/null +++ b/tests.cpp @@ -0,0 +1,602 @@ +/* + * Copyright April 2023 + * Author: Dong Anakin + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * + * See the License for the specific language governing permissions and + * limitations under the License. + * + */ + +// Include Files +#include "savitzky_golay.hpp" +#include +#include +#define VERIFY_EQUAL(a, b) \ + do \ + { \ + auto ret = (a - b).array().abs().maxCoeff(); \ + if (ret < 1e-6) \ + std::cout << __func__ << ":" << __LINE__ << " OK" << std::endl; \ + else \ + std::cerr << __func__ << ":" << __LINE__ << " Failed" << std::endl; \ + } while (false) +#define VERIFY_EMPTY(a) \ + do \ + { \ + int ret = a.rows() * a.cols(); \ + if (ret == 0) \ + std::cout << __func__ << ":" << __LINE__ << " OK" << std::endl; \ + else \ + std::cerr << __func__ << ":" << __LINE__ << " Failed" << std::endl; \ + } while (false) +void test_correlate01() +{ + Eigen::RowVectorXd input(2); + input << 1, 2; + Eigen::RowVectorXd weights(1); + weights << 2; + Eigen::RowVectorXd expected(2); + expected << 2, 4; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} +void test_correlate02() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(1); + weights << 1; + Eigen::RowVectorXd expected(3); + expected << 1, 2, 3; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} +void test_correlate03() +{ + Eigen::RowVectorXd input(1); + input << 1; + Eigen::RowVectorXd weights(2); + weights << 1, 1; + Eigen::RowVectorXd expected(1); + expected << 2; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} + +void test_correlate04() +{ + Eigen::RowVectorXd input(2); + input << 1, 2; + Eigen::RowVectorXd weights(2); + weights << 1, 1; + Eigen::RowVectorXd tcor(2); + tcor << 2, 3; + Eigen::RowVectorXd tcov(2); + tcov << 3, 4; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate05() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(2); + weights << 1, 1; + Eigen::RowVectorXd tcor(3); + tcor << 2, 3, 5; + Eigen::RowVectorXd tcov(3); + tcov << 3, 5, 6; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate06() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 3; + Eigen::RowVectorXd tcor(3); + tcor << 9, 14, 17; + Eigen::RowVectorXd tcov(3); + tcov << 7, 10, 15; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate07() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd expected(3); + expected << 5, 8, 11; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} + +void test_correlate08() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(3); + weights << 1, 2, -1; + Eigen::RowVectorXd tcor(3); + tcor << 1, 2, 5; + Eigen::RowVectorXd tcov(3); + tcov << 3, 6, 7; + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate09() +{ + Eigen::RowVectorXd input(0); + Eigen::RowVectorXd weights(2); + weights << 1, 1; + Eigen::RowVectorXd tcor(0); + Eigen::RowVectorXd tcov(0); + + auto r = scipy::ndimage::correlate1d(input, weights); + // savgol_print_matrix(r); + VERIFY_EMPTY(r); + + auto c = scipy::ndimage::convolve1d(input, weights); + // savgol_print_matrix(c); + + VERIFY_EMPTY(c); +} + +void test_correlate17() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(2); + weights << 1, 1; + Eigen::RowVectorXd tcor(3); + tcor << 3, 5, 6; + Eigen::RowVectorXd tcov(3); + tcov << 2, 3, 5; + + auto r = scipy::ndimage::correlate1d(input, weights, "reflect", 0.0f, -1); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "reflect", 0.0f, -1); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate22() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd expected(6); + expected << 10, 8, 10, 11, 16, 17; + + auto r = scipy::ndimage::correlate1d(input, weights, "wrap"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights, "wrap"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, expected); +} + +void test_correlate23() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd expected(6); + expected << 5, 8, 10, 11, 16, 22; + + auto r = scipy::ndimage::correlate1d(input, weights, "nearest"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights, "nearest"); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} + +void test_correlate24() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd tcor(6); + tcor << 8, 10, 11, 16, 22, 24; + Eigen::RowVectorXd tcov(6); + tcov << 4, 5, 8, 10, 11, 16; + + auto r = scipy::ndimage::correlate1d(input, weights, "nearest", 0.0f, -1); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "nearest", 0.0f, -1); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate25() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd tcor(6); + tcor << 4, 5, 8, 10, 11, 16; + Eigen::RowVectorXd tcov(6); + tcov << 8, 10, 11, 16, 22, 24; + + auto r = scipy::ndimage::correlate1d(input, weights, "nearest", 0.0f, 1); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "nearest", 0.0f, 1); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate26() +{ + Eigen::RowVectorXd input(1); + input << 1; + Eigen::RowVectorXd weights(5); + weights << 1, 1, 1, 1, 1; + Eigen::RowVectorXd expected(1); + expected << 5; + + auto r = scipy::ndimage::correlate1d(input, weights, "mirror"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, expected); + + auto c = scipy::ndimage::convolve1d(input, weights, "mirror"); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, expected); +} + +void test_correlate27() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, 1; + Eigen::RowVectorXd tcor(6); + tcor << 9, 8, 10, 11, 16, 21; + Eigen::RowVectorXd tcov(6); + tcov << 9, 8, 10, 11, 16, 21; + + auto r = scipy::ndimage::correlate1d(input, weights, "constant", 5.0f); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "constant", 5.0f); + // savgol_print_matrix(c); + + VERIFY_EQUAL(c, tcov); +} + +void test_correlate28() +{ + Eigen::RowVectorXd input(6); + input << 1, 2, 3, 2, 4, 6; + Eigen::RowVectorXd weights(3); + weights << 1, 2, -1; + Eigen::RowVectorXd tcor(6); + tcor << 5, 2, 6, 3, 4, 11; + Eigen::RowVectorXd tcov(6); + tcov << -1, 6, 6, 5, 12, 13; + + auto r = scipy::ndimage::correlate1d(input, weights, "constant", 5.0f); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "constant", 5.0f); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_correlate30() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(8); + weights << 1, 0, 0, 0, 0, 0, 0, 0; + Eigen::RowVectorXd tcor(3); + tcor << 1, 1, 1; + Eigen::RowVectorXd tcov(3); + tcov << 3, 3, 3; + + auto r = scipy::ndimage::correlate1d(input, weights, "nearest"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "nearest"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_correlate31() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(8); + weights << 1, 0, 0, 0, 0, 0, 0, 0; + Eigen::RowVectorXd tcor(3); + tcor << 3, 1, 2; + Eigen::RowVectorXd tcov(3); + tcov << 2, 3, 1; + + auto r = scipy::ndimage::correlate1d(input, weights, "wrap"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "wrap"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_correlate32() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(8); + weights << 1, 0, 0, 0, 0, 0, 0, 0; + Eigen::RowVectorXd tcor(3); + tcor << 3, 3, 2; + Eigen::RowVectorXd tcov(3); + tcov << 2, 1, 1; + + auto r = scipy::ndimage::correlate1d(input, weights, "reflect"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "reflect"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_correlate33() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(8); + weights << 1, 0, 0, 0, 0, 0, 0, 0; + Eigen::RowVectorXd tcor(3); + tcor << 1, 2, 3; + Eigen::RowVectorXd tcov(3); + tcov << 1, 2, 3; + + auto r = scipy::ndimage::correlate1d(input, weights, "mirror"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "mirror"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_correlate34() +{ + Eigen::RowVectorXd input(3); + input << 1, 2, 3; + Eigen::RowVectorXd weights(8); + weights << 1, 0, 0, 0, 0, 0, 0, 0; + Eigen::RowVectorXd tcor(3); + tcor << 0, 0, 0; + Eigen::RowVectorXd tcov(3); + tcov << 0, 0, 0; + + auto r = scipy::ndimage::correlate1d(input, weights, "constant"); + // savgol_print_matrix(r); + VERIFY_EQUAL(r, tcor); + + auto c = scipy::ndimage::convolve1d(input, weights, "constant"); + // savgol_print_matrix(c); + VERIFY_EQUAL(c, tcov); +} + +void test_savgol_coeffs_01() +{ + Eigen::RowVectorXd expected(5); + expected << 0.77142857, -0.18571429, -0.57142857, -0.38571429, 0.37142857; + + int pos = 4; + auto coeff = scipy::signal::savgol_coeffs(5, 2, 1, 1.0f, &pos, "conv"); + // savgol_print_matrix(coeff); + VERIFY_EQUAL(coeff, expected); +} + +void test_savgol_coeffs_02() +{ + Eigen::RowVectorXd expected(5); + expected << 0.37142857, -0.38571429, -0.57142857, -0.18571429, 0.77142857; + + int pos = 4; + auto coeff = scipy::signal::savgol_coeffs(5, 2, 1, 1.0f, &pos, "dot"); + // savgol_print_matrix(coeff); + VERIFY_EQUAL(coeff, expected); +} + +void test_savgol_filter_01() +{ + Eigen::RowVectorXd x(9); + x << 2, 2, 5, 2, 1, 0, 1, 4, 9; + + Eigen::RowVectorXd expected(9); + expected << 1.65714286, 3.17142857, 3.54285714, 2.85714286, 0.65714286, 0.17142857, 1.0, 4.0, + 9.0; + + auto y = scipy::signal::savgol_filter(x, 5, 2); + // savgol_print_matrix(y); + VERIFY_EQUAL(y, expected); +} + +void test_savgol_filter_02() +{ + Eigen::RowVectorXd x(9); + x << 2, 2, 5, 2, 1, 0, 1, 4, 9; + + Eigen::RowVectorXd expected(9); + expected << 1.48571429, 3.02857143, 3.54285714, 2.85714286, 0.65714286, 0.17142857, 1., + 5.02857143, 6.94285714; + + auto y = scipy::signal::savgol_filter(x, 5, 2, 0, 1.0, "mirror"); + // savgol_print_matrix(y); + VERIFY_EQUAL(y, expected); +} + +void test_savgol_filter_03() +{ + Eigen::RowVectorXd x(9); + x << 2, 2, 5, 2, 1, 0, 1, 4, 9; + + Eigen::RowVectorXd expected(9); + expected << 1.71428571, -0.42857143, -1.14285714, -0.85714286, 1.14285714, 1.42857143, 2., + -1.42857143, -3.14285714; + + auto y = scipy::signal::savgol_filter(x, 5, 3, 2, 1.0, "mirror"); + // savgol_print_matrix(y); + VERIFY_EQUAL(y, expected); +} + +void test_savgol_filter_04() +{ + Eigen::RowVectorXd x(9); + x << 2, 2, 5, 2, 1, 0, 1, 4, 9; + + Eigen::RowVectorXd expected(9); + expected << -0.14285714, -0.64285714, -1.14285714, -0.85714286, 1.14285714, 1.42857143, 2., 2., + 2.; + + auto y = scipy::signal::savgol_filter(x, 5, 3, 2, 1.0); + // savgol_print_matrix(y); + VERIFY_EQUAL(y, expected); +} + +void test_numpy_polyfit_01() +{ + Eigen::RowVectorXd x(6); + x << 0.0, 1.0, 2.0, 3.0, 4.0, 5.0; + Eigen::RowVectorXd y(6); + y << 0.0, 0.8, 0.9, 0.1, -0.8, -1.0; + Eigen::RowVectorXd expected(4); + expected << 0.08703704, -0.81349206, 1.69312169, -0.03968254; + + Eigen::RowVectorXd z = numpy::polyfit(x, y, 3); + // savgol_print_matrix(z); + VERIFY_EQUAL(z, expected); +} + +int main(int argc, char **argv) +{ + test_correlate01(); + test_correlate02(); + test_correlate03(); + test_correlate04(); + test_correlate05(); + test_correlate06(); + test_correlate07(); + test_correlate08(); + test_correlate09(); + test_correlate17(); + test_correlate22(); + test_correlate23(); + test_correlate24(); + test_correlate25(); + test_correlate26(); + test_correlate27(); + test_correlate28(); + test_correlate30(); + test_correlate31(); + test_correlate32(); + test_correlate33(); + test_correlate34(); + test_numpy_polyfit_01(); + test_savgol_coeffs_01(); + test_savgol_coeffs_02(); + test_savgol_filter_01(); + test_savgol_filter_02(); + test_savgol_filter_03(); + test_savgol_filter_04(); + return 0; +} -- Gitee