Python TsFresh特征的Java实现——spkt_welch_density

 TsFresh（TimeSeries Fresh）是一个Python第三方工具包。它可以方便地对时间序列数据进行处理，获得大量的特征。这些特征可以用以训练分类器，以高效地实现对时间序列数据的分类、识别等。然而，在工程实现时，更多地是采用Java等语言，这需要利用Java实现对TsFresh的特征进行直接计算，故需要对TsFresh的某些特征进行深入地分析，并在Java语言下实现。

特征spkt_welch_density简介

命令格式：spkt_welch_density(x, param)

计算时间序列的功率谱密度在特征频率点上的值作为输出特征。

参数：

• x：时间序列，数据类型：numpy.ndarray
• param：参数字典{"coeff": q}，其中q为表示输出特征的功率谱密度数组的索引值

特征spkt_welch_density计算原理

1. 根据Hanning窗的表达式，形成与时间序列点数相匹配的Hanning权矢量
2. 去除时间序列x中的直流分量，形成新的序列x_detrend
3. 用Hanning权矢量对序列x_detrend加权，形成加权后序列weighted_data
4. 对加权后序列weighted_data做单边Fourier变换，得其单边频谱fx
5. 利用单边频谱fx共轭相乘并进行归一化处理后，获得单边功率谱pxx
6. 利用参数Coeff选择单边功率谱pxx相应索引处的值作为特征输出

Java实现

构造Hanning窗

public static double[] Window_Hann(int width_Window) {
	double dTheta = Math.PI*2.0/width_Window;
	double[] angle = new double[width_Window];
	for (int idx = 0; idx < width_Window; idx++) {
    	angle[idx] = -Math.PI + idx * dTheta;
	}
	System.out.println(Arrays.toString(angle));
	double[] weights = new double[width_Window];
	for (int idx = 0; idx < width_Window; idx++) {
		weights[idx] = 0.5 + 0.5 * Math.cos(angle[idx]);
	}
	return weights;
}

在进行Fourier变换之前需要对时间序列做去直流处理。

由于时间序列的长度不能够保证是，所以无法使用Apache Commons Math 3.6.1 API中的FastFourierTransform，通常在网上能够找到的Fourier变换程序也是需要对数据长度要求为，因此，需要对此进行独立编写。

Complex[] fx = new Complex[n_x/2 + 1];
int n_f = fx.length;
for (int idx_n = 0; idx_n < n_f; idx_n++) {
	fx[idx_n] = new Complex(0);
	for (int idx_k = 0; idx_k < n_x; idx_k++) {
		Complex cplx_exp = new Complex(Math.cos(angle[idx_k]*idx_n), 
            -Math.sin(angle[idx_k]*idx_n));
		fx[idx_n] = fx[idx_n].add(cplx_exp.multiply(new Complex(weighted_data[idx_k])));
	}
}
		

单边功率谱不是简单地等于单边频谱的平方，而应该是除零频之外的应是单边频谱平方的两倍

double[] pxx = new double[n_f]; 
for (int idx_n = 0; idx_n < n_f; idx_n++) {
	if (idx_n == 0) {
		pxx[idx_n] = fx[idx_n].conjugate().multiply(fx[idx_n]).getReal()/sum_wgt;
	} else {
		pxx[idx_n] = fx[idx_n].conjugate().multiply(fx[idx_n]).getReal()*2/sum_wgt;
	}
}