频域以及时-频域健康指标原理以及c与Python实现

频域健康指标

频域健康特征指标如下表：

部分代码实现如下：

def spectral_gini_index(data, sampling_rate = 1.0):
    # Compute the power spectral density
    freqs, psd = welch(data)

    # Sort PSD in ascending order
    sorted_psd = np.sort(psd)
    # Compute cumulative sum and CDF
    data = np.arange(0, len(sorted_psd))
    shang_sort = data[::-1]

    number_son = np.sum(sorted_psd*shang_sort)
    # Compute spectral Gini index
    
    spectral_gini = number_son / ((len(data)-1)*np.sum(sorted_psd))

    return spectral_gini
#谱峭度
def spectral_kurtosis(data):
    fs = len(data)
    Pxx = []
    nfft = 4
    for i in range(fs//nfft):                                      # 窗与窗之间数据不重叠
        w = np.hanning(nfft)                                      # 加汉尼窗
        p = np.abs(fft(w*data[i*nfft:(i+1)*nfft]))**2/fs/norm(w)**2  # 计算功率谱
        Pxx.append(p) 
    #print(Pxx)
    Pxx1 = np.mean(Pxx,axis=0)     #功率谱密度
    #print(Pxx1)
    f1   = np.array([k*fs/nfft for k in range(nfft//2)])
    #print(f1)
    f_mean =  np.mean(f1)
   # print(f_mean)
    sk = np.sum(((f1 - f_mean)**4)*Pxx1[nfft//2:]) / np.power(np.sum(Pxx1[nfft//2:]),2)
    #sk = np.sum(((f1[nfft//2:] - f_mean)**4)*Pxx1[nfft//2:]) / np.power(np.sum(Pxx1[nfft//2:]),2)
    return sk

C语言部分代码如下：

//傅里叶变换
void change()
{
	complex temp;
	unsigned short i=0,j=0,k=0;
	double t;
	for(i=0;i<size;i++)
	{
		k=i;
		j=0;