模糊熵（FE）计算的C++ 实现，及参考MATLAB代码

模糊熵（FE）计算的C++ 实现，及参考MATLAB代码

基本原理

算法基本原理

算法基本于参照下图，模糊熵的意义请自行搜索，代码里面有注释，具体看代码就行。

为什么用C++写

时间复杂度基本与MATLAB一样，测试中，MATLAB占用CPU（50%）约为C++（5%左右）的十倍，计算速度为C++的四倍，多核并行计算的时候C++可能会更快一些（未做测试）。

MATLAB代码

souce

function [Out_FuzEn,P] = FuzEn(x,m,r,n,tau)
	%
	% This function calculates fuzzy entropy (FuzEn) of a univariate signal x
	%
	% Inputs:
	%
	% x: univariate signal - a vector of size 1 x N (the number of sample points)
	% m: embedding dimension
	% r: threshold (it is usually equal to 0.15 of the standard deviation of a signal - because we normalize signals to have a standard deviation of 1, here, r is usually equal to 0.15)
	% n: fuzzy power (it is usually equal to 2)
	% tau: time lag (it is usually equal to 1)
	%
	% Outputs:
	%
	% Out_FuzEn: scalar quantity - the FuzEn of x
	% P: a vector of length 2 - [the global quantity in dimension m, the global quantity in dimension m+1]
	%
	%
	% Ref:
	% [1] H. Azami and J. Escudero, "Refined Multiscale Fuzzy Entropy based on Standard Deviation for Biomedical Signal Analysis", Medical & Biological Engineering &
	% Computing, 2016.
	% [2] W. Chen, Z. Wang, H. Xie, and W. Yu,"Characterization of surface EMG signal based on fuzzy entropy", IEEE Transactions on neural systems and rehabilitation engineering, vol. 15, no. 2, pp.266-272, 2007.
	%
	%
	% If you use the code, please make sure that you cite references [1] and [2].
	%
	% Hamed Azami and Javier Escudero Rodriguez
	% hamed.azami@ed.ac.uk and javier.escudero@ed.ac.uk
	%
	%  7-September-16
	%%
	
	
	if nargin == 4, tau = 1; end
	if nargin == 3, n = 2; tau=1; end
	if tau > 1, x = downsample(x, tau); end
	
	N = length(x);
	P = zeros(1,2);
	% 用来暂存x的内容
	xMat = zeros(m+1,N-m);
	
	% 将时间序列分为k = n - m + 1个序列（？
	for i = 1:m+1
	    xMat(i,:) = x(i:N-m+i-1);
	end
	
	% 计算每个序列与所有k个序列之间的距离， 并列表
	for k = m:m+1
	    count = zeros(1,N-m);
	    tempMat = xMat(1:k,:);
	    
	
	    for i = 1:N-k
	        % calculate Chebyshev distance without counting self-matches
	        % dij = max|x_{i + k}(t) - x_{j + k}(t)|
	        dist = max(abs(tempMat(:,i+1:N-m) - repmat(tempMat(:,i),1,N-m-i)));
	        test = repmat(tempMat(:,i),1,N-m-i);
	        % D_{ij}^m = exp(\frac{- (d_{ij}^m)^n})
	        DF=exp((-dist.^n)/r);
	
	        % 对除自身外所有隶属度求平均
	        count(i) = sum(DF)/(N-m);
	    end
	    
	    P(k-m+1) = sum(count)/(N-m);
	end
	
	Out_FuzEn = log(P(1)/P(2));
end

C++代码

其中部分是我自己写的辅助函数，如sum，downsample是写了个简单的降采样方式，需要自取。

#include <iostream>
#include <cmath>
#include <cstring>
#include <fstream>
#include <time.h>
using namespace std;

double sum(double* data , int N){
	// 数组求和
	double result = 0;
	for(int i = 0 ; i < N ; i++){
		result += data[i];
	}
	return result;
}

double* downsample(double* x ,long N ,int tau){
	long new_N = N / tau;
	double* result = new double[new_N];
	for(int i = 0 ; i * tau < new_N ; i++){
		result[i] = x[i * tau];
	}
	return result;
}

double fussyEntropy(double* x ,long N ,int m ,double r ,double n ,int tau)
{
	/*
		x: 信号
		N: x的长度
		m: embedding dimension
		r: threshold (it is usually equal to 0.15 of the standard deviation of a signal - because we normalize signals to have a standard deviation of 1, here, r is usually equal to 0.15)
		n: fuzzy power (it is usually equal to 2)
	*/

	if(tau > 1){
		x = downsample(x , N , tau);
		N = N / tau + 1;
	}

	double P[2] = {0};

	// 计算每个序列与所有k个序列之间的距离并列表
	for(int k = m ; k <= m+1 ; k++){

		double* count = new double[N - m];
		memset(count , 0 , sizeof(double) * (N - m));
		double *DF = new double[N - k];

		for(int i = 0 ; i < N - k ; i++){
			memset(DF , 0 , sizeof(double) * (N - k));
			//比较第i个序列和所有序列
			for(int j = i ; j < N - k ; j++){
				// 逐个序列进行比较
				double dist_temp = 0 , dist_max = 0;
				if(i == j){ DF[j] = 1 ; continue;} // 同一个序列 直接置零
				else{
					for(int a = 0 ; a < k ; a++){
						// 第i段 - 第j段
						dist_temp = abs(x[a + i] - x[a + j]);
						dist_max = dist_max > dist_temp ? dist_max : dist_temp; //找到最大距离
					}
				}
				DF[j] = exp(-pow(dist_max , n)/r);
			}
			count[i] = (sum(DF , N - k) - 1)/ (N - k - 1); // -1 为了抹去对对本身的距离计算
		}
		P[k - m] = sum(count , N - k) / (N - k); 	

		delete []DF;
		delete []count;
	}
	
	double result = log(P[0]) - log(P[1]);
	return result;
}

double *get_signal(int N){
	// 读取测试文件数据
	double *signal = new double[N];
	memset(signal , 0 , sizeof(double) * N);

	ifstream inFile;
	inFile.open("white_noise.csv");
	
	if (!inFile) {
		cerr << "Unable to open file datafile.txt";
		exit(1);   // call system to stop
	}

	double number;
	char c;
	double value;
	int n = 0;
	while(inFile >> number && inFile >> c && inFile >> value){
		signal[n++] = value;
		if(n == N) break;
	}
	inFile.close();
	
	return signal;
}

int main(){
	// 读取测试文件数据
	int N = 10000;
	double *signal = get_signal(N);
	
	// 计算
	cout << fussyEntropy(signal , N , 2 , 0.15 , 2 , 1) << endl;
	delete []signal;
    return 0;  
}