diff --git a/230124-sampen/Makefile b/230124-sampen/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..2e6dd9a958fa686828413109a0294adce931b1b0 --- /dev/null +++ b/230124-sampen/Makefile @@ -0,0 +1,20 @@ +cxx = g++ +src = $(wildcard *.cpp) +obj = $(patsubst %.cpp, %.o, $(src)) +target = test + +all: $(target) + +$(target): $(obj) + $(cxx) -o $@ $^ + +%.o: %.cpp + $(cxx) -c -O2 $< + +# Windows +clean: + @ del $(obj) $(target).exe + +# Linux +# clean: +# rm -rf $(obj) $(target) diff --git a/230124-sampen/Readme.md b/230124-sampen/Readme.md new file mode 100644 index 0000000000000000000000000000000000000000..679855df7a3fa63f9ff989d987be40457c7fe108 --- /dev/null +++ b/230124-sampen/Readme.md @@ -0,0 +1,21 @@ +# 样本熵 + +## 程序说明 + +test.cpp是用来测试的，要用c++编译器编译 + +具体代码在sampen.cpp，和sampen.h对应，可以用c编译器编译，所以支持C和C++ + +## 编译说明 + +mingw64的g++版本为8.5.0，理论兼容更新的版本 +下载地址 + +vs用的2022版(17.4)，不确定是不是兼容老版本 +下载地址 + +对于mingw64，可以使用mingw32-make命令编译，也可以用g++ test.cpp -o test -D ALL_IN_ONE -O3编译 + +对于vs，打开工程，选择架构等，然后点击运行就可以了 + +如果你用Linux，那就要改一下Makefile里的clean，然后运行make就行了，其余和mingw64相同，不过具体我没测试过 diff --git a/230124-sampen/sampen.cpp b/230124-sampen/sampen.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4cb84428f61abdad2ea5d10bd1625ee308b11ffd --- /dev/null +++ b/230124-sampen/sampen.cpp @@ -0,0 +1,95 @@ +#include "sampen.h" +#include + +static double step(double *X, int N, int m, double r) +{ + int Bi = 0; + for (int i = 0; i <= N - m; i++) { + for (int j = 0; j <= N - m; j++) { + if (i != j) { + double D = fabs(X[i] - X[j]); + for (int k = 1; k < m; k++) { + double t = fabs(X[i + k] - X[j + k]); + if (D < t) + D = t; + } + if (D <= r) + Bi++; + } + } + } + return 1.0 * Bi / (N - m) / (N - m + 1); +} + +double SampEn(double *X, int N, int m, double r) +{ + double B = step(X, N, m, r); + if (B == 0) + return 0; + double A = step(X, N, m + 1, r); + if (A == 0) + return 0; + return -log(A / B); +} + +double FastSampEn(double *X, int N, int m, double r) +{ + int Ai = 0, Bi = 0; + int LoopsSub1 = N - m; + for (int i = 0; i <= LoopsSub1; i++) { + for (int j = 0; j <= LoopsSub1; j++) { + if (i != j) { + double D = fabs(X[i] - X[j]); + for (int k = 1; k < m; k++) { + double t = fabs(X[i + k] - X[j + k]); + if (D < t) + D = t; + } + if (D <= r) + Bi++; + if (i != LoopsSub1 && j != LoopsSub1) { + double t = fabs(X[i + m] - X[j + m]); + if (D < t) + D = t; + if (D <= r) + Ai++; + } + } + } + } + double B = 1.0 * Bi / (N - m) / (N - m + 1); + double A = 1.0 * Ai / (N - m - 1) / (N - m); + if (B == 0 || A == 0) + return 0; + return -log(A / B); +} + +double FastSampEn_m2(double *X, int N, double r) +{ + int Ai = 0, Bi = 0; + int LoopsSub1 = N - 2; + for (int i = 0; i <= LoopsSub1; i++) { + for (int j = 0; j <= LoopsSub1; j++) { + if (i != j) { + double D = fabs(X[i] - X[j]); + double t = fabs(X[i + 1] - X[j + 1]); + if (D < t) + D = t; + if (D <= r) + Bi++; + if (i != LoopsSub1 && j != LoopsSub1) { + double t = fabs(X[i + 2] - X[j + 2]); + if (D < t) + D = t; + if (D <= r) + Ai++; + } + } + } + } + double B = 1.0 * Bi / (N - 2) / (N - 1); + double A = 1.0 * Ai / (N - 3) / (N - 2); + if (B == 0 || A == 0) + return 0; + return -log(A / B); +} diff --git a/230124-sampen/sampen.h b/230124-sampen/sampen.h new file mode 100644 index 0000000000000000000000000000000000000000..2d39bdf1b2df13b837dfd522270022b39535d976 --- /dev/null +++ b/230124-sampen/sampen.h @@ -0,0 +1,8 @@ +#pragma once + +// 正常简单优化的样本熵 +double SampEn(double *X, int N, int m, double r); +// 经过优化的样本熵 +double FastSampEn(double *X, int N, int m, double r); +// 经过优化的且m取2的样本熵 +double FastSampEn_m2(double *X, int N, double r); diff --git a/230124-sampen/sampen.sln b/230124-sampen/sampen.sln new file mode 100644 index 0000000000000000000000000000000000000000..a90cb2e9473012628cbb009c831a260b56af3e1f --- /dev/null +++ b/230124-sampen/sampen.sln @@ -0,0 +1,31 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio Version 17 +VisualStudioVersion = 17.4.33213.308 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "sampen", "sampen.vcxproj", "{D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|x64 = Debug|x64 + Debug|x86 = Debug|x86 + Release|x64 = Release|x64 + Release|x86 = Release|x86 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Debug|x64.ActiveCfg = Debug|x64 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Debug|x64.Build.0 = Debug|x64 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Debug|x86.ActiveCfg = Debug|Win32 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Debug|x86.Build.0 = Debug|Win32 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Release|x64.ActiveCfg = Release|x64 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Release|x64.Build.0 = Release|x64 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Release|x86.ActiveCfg = Release|Win32 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5}.Release|x86.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection + GlobalSection(ExtensibilityGlobals) = postSolution + SolutionGuid = {FD2BB481-ABF1-4C26-A8BF-9129B48D0CD4} + EndGlobalSection +EndGlobal diff --git a/230124-sampen/sampen.vcxproj b/230124-sampen/sampen.vcxproj new file mode 100644 index 0000000000000000000000000000000000000000..62260cb110da0dd1092d17a0b899fec152e599c9 --- /dev/null +++ b/230124-sampen/sampen.vcxproj @@ -0,0 +1,119 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + Debug + x64 + + + Release + x64 + + + + 17.0 + {D1999EF9-B18C-41EA-BF88-CB7AD1F4F5C5} + Win32Proj + 10.0 + + + + Application + true + v143 + + + Application + false + v143 + + + Application + true + v143 + + + Application + false + v143 + + + + + + + + + + + + + + + + + + + + + true + $(ProjectName)test + + + true + $(ProjectName)test + + + $(ProjectName)test + + + $(ProjectName)test + + + + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDebugDLL + Level3 + ProgramDatabase + Disabled + + + MachineX86 + true + Console + + + + + WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions) + MultiThreadedDLL + Level3 + ProgramDatabase + + + MachineX86 + true + Console + true + true + + + + + + + + + + + + + \ No newline at end of file diff --git a/230124-sampen/sampen.vcxproj.filters b/230124-sampen/sampen.vcxproj.filters new file mode 100644 index 0000000000000000000000000000000000000000..336f4ad4a087eceb94778098b7059846afefb7cb --- /dev/null +++ b/230124-sampen/sampen.vcxproj.filters @@ -0,0 +1,30 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hh;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav + + + + + Source Files + + + Source Files + + + + + Header Files + + + \ No newline at end of file diff --git a/230124-sampen/test.cpp b/230124-sampen/test.cpp new file mode 100644 index 0000000000000000000000000000000000000000..171d6fac15601d2dd59ed2e024a0aafa11e3f199 --- /dev/null +++ b/230124-sampen/test.cpp @@ -0,0 +1,48 @@ +#include +#include +#include +#ifdef ALL_IN_ONE // 编译时加上-D ALL_IN_ONE可以不用Makefile，便于测试 +#include "sampen.cpp" +#else +#include "sampen.h" +#endif + +int main() +{ + std::vector x; + for (int i = 0; i < 17; i++) { + x.push_back(85); + x.push_back(80); + x.push_back(89); + } + + /* 必须是0.0008507018803128114舍入后的结果 */ + std::cout << SampEn(x.data(), x.size(), 2, 3) << '\n'; + std::cout << FastSampEn(x.data(), x.size(), 2, 3) << '\n'; + std::cout << FastSampEn_m2(x.data(), x.size(), 3) << '\n'; + + // 继续添加30000个数据 + for (int i = 0; i < 10000; i++) { + x.push_back(85); + x.push_back(80); + x.push_back(89); + } + + double se; + clock_t t; + /* 进行速度比较 */ + t = clock(); + se = SampEn(x.data(), x.size(), 2, 3); + t = clock() - t; + std::cout << se << ", time = " << t << " ms\n"; + + t = clock(); + se = FastSampEn(x.data(), x.size(), 2, 3); + t = clock() - t; + std::cout << se << ", time = " << t << " ms\n"; + + t = clock(); + se = FastSampEn_m2(x.data(), x.size(), 3); + t = clock() - t; + std::cout << se << ", time = " << t << " ms\n"; +}