CEEMDAN-多尺度熵（MSE）-SVM滚动轴承故障诊断-CFANZ编程社区

1、完全集合经验模态分解CEEMDAN本文就不做多介绍了，此算法是EEMD和EMD改进后得到的，一定程度上有效的抑制了模态混叠等问题。

加入经 EMD 分解后含辅助噪声的 IMF 分量，而不是将高斯白噪声信号直接添加在原始信号中；
EEMD 分解和 CEEMD 分解是将经验模态分解后得到的模态分量进行总体平均，CEEMDAN 分解则在得到的第一阶 IMF分量后就进行总体平均计算，得到最终的第一阶 IMF分量，然后对残余部分重复进行如上操作，这样便有效地解决了白噪声从高频到低频的转移传递问题。

2、多尺度熵

多尺度熵MSE的基本原理就是对时间序列进行粗粒化或下采样，即主要是在越来越粗略的时间分辨率下分析时间序列。

def MSE(signal , max_scale):
    result = []
    length = len(signal)
    std = np.std(signal)

    for scale in range(1 , max_scale + 1):
        # 确定截取的长度
        length = int(len(signal) / scale) - 1
        # 分段取平均
        scale_i = signal[ : len(signal) : scale][:length]
        for i in range(1,scale):
            scale_i = scale_i + signal[i: len(signal) : scale][:length]
        scale_i = scale_i / scale
        #计算样本熵
        result.append(Sample_Entropy(scale_i, std ,r = 0.15))
        print("scale:" , scale, 'SampEn' , result[-1])
    return result

再附上样本熵代码

def Sample_Entropy(x, m, r=0.15):
    
    
    x = np.array(x)
    # 检查x是否为一维数据
    if x.ndim != 1:
        raise ValueError("x的维度不是一维")
    # 计算x的行数是否小于m+1
    if len(x) < m+1:
        raise ValueError("len(x)小于m+1")
    # 将x以m为窗口进行划分
    entropy = 0  
    for temp in range(2):
        X = []
        for i in range(len(x)-m+1-temp):
            X.append(x[i:i+m+temp])
        X = np.array(X)
        # 计算X任意一行数据与所有行数据对应索引数据的差值绝对值的最大值
        D_value = []  
        for index1, i in enumerate(X):
            sub = []
            for index2, j in enumerate(X):
                if index1 != index2:
                    sub.append(max(np.abs(i-j)))
            D_value.append(sub)
        # 计算阈值
        F = r*np.std(x, ddof=1)
        # 判断D_value中的每一行中的值比阈值大的个数除以len(x)-m+1的比例
        num = np.sum(D_value>F, axis=1)/(len(X)-m+1-temp)
        # 计算num的对数平均值
        Lm = np.average(np.log(num))
        entropy = abs(entropy) - Lm

    return entropy

注：r你可以自己设置，具体范围通常取0.1-0.25.

3.利用相关系数对CEEMDAN分解的信号进行重构。

def P(x,y):
    x1 = np.sum(x) / len(x)
    y1 = np.sum(y) / len(y)
    p = np.sum((x-x1)*(y-y1)) / (np.sqrt(np.sum((x-x1)**2))*np.sqrt(np.sum((y-y1)**2)))
    return p

k1 = []
k2 = []
k3 = []
k4 = []
k5 = []
k6 = []
k7 = []
k8 = []
k9 = []
k10 = []
for i in range(400):
    data = data1[i]
    ceemdan = CEEMDAN()
    eIMFs = ceemdan.ceemdan(data)
    nIMFs = eIMFs.shape[0]
    Y = []
    for j in range(nIMFs):
        p = P(data,eIMFs[j])
        if p > 0.1:
           Y.append(eIMFs[j])
    y_c = np.zeros(2048)
    for k in range(len(Y)):
        y_c = y_c+Y[k]
    print(y_c)
    mse = MSE(y_c,10)
    k1.append(mse[0])
    k2.append(mse[1])
    k3.append(mse[2])
    k4.append(mse[3])
    k5.append(mse[4])
    k6.append(mse[5])
    k7.append(mse[6])
    k8.append(mse[7])
    k9.append(mse[8])
    k10.append(mse[9])

P为计算得到的相关系数，当它大于0.1就保留下来。

注意：本次是计算10个熵值，还没学会怎么将多维数组直接保存到表格中，于是就保存了十个文本，然后手动操作了，哈哈哈，哪天学会了再分享，就不会这么麻烦了。

4、选取实验数据损伤为0.018、转速1797r/min的西储大学轴承故障诊断数据，外圈故障、内圈故障、滚动体故障和正常状态分别取100个数据，每个数据长度为2048个点。

5、完整代码

from PyEMD import CEEMDAN
import numpy as np
import matplotlib.pyplot as plt
import os
import glob


path = r'E:/CEEMDAN-SSA-SVM/0.018_0hp/'
def read_txt(path):
    cate = [path + x for x in os.listdir(path) if os.path.isdir(path + x)]
    imgs = []
    labels = []
    for idx, folder in enumerate(cate):
        for im in glob.glob(folder + '/*.txt'):
            img = np.loadtxt(im)
            imgs.append(img)
            labels.append(idx)
    return np.asarray(imgs, np.float32), np.asarray(labels, np.int32)


data1, label = read_txt(path)
print(label)  


def Sample_Entropy(x, m, r=0.15):
    """
    样本熵
    m 滑动时窗的长度
    r 阈值系数 取值范围一般为：0.1~0.25
    """
    # 将x转化为数组
    x = np.array(x)
    # 检查x是否为一维数据
    if x.ndim != 1:
        raise ValueError("x的维度不是一维")
    # 计算x的行数是否小于m+1
    if len(x) < m+1:
        raise ValueError("len(x)小于m+1")
    # 将x以m为窗口进行划分
    entropy = 0  # 近似熵
    for temp in range(2):
        X = []
        for i in range(len(x)-m+1-temp):
            X.append(x[i:i+m+temp])
        X = np.array(X)
        # 计算X任意一行数据与所有行数据对应索引数据的差值绝对值的最大值
        D_value = []  # 存储差值
        for index1, i in enumerate(X):
            sub = []
            for index2, j in enumerate(X):
                if index1 != index2:
                    sub.append(max(np.abs(i-j)))
            D_value.append(sub)
        # 计算阈值
        F = r*np.std(x, ddof=1)
        # 判断D_value中的每一行中的值比阈值大的个数除以len(x)-m+1的比例
        num = np.sum(D_value>F, axis=1)/(len(X)-m+1-temp)
        # 计算num的对数平均值
        Lm = np.average(np.log(num))
        entropy = abs(entropy) - Lm

    return entropy


def MSE(signal , max_scale:int = 20):
    result = []
    length = len(signal)
    # std = np.std(signal)

    for scale in range(1 , max_scale + 1):
        # 确定截取的长度
        length = int(len(signal) / scale) - 1
        # 分段取平均
        scale_i = signal[ : len(signal) : scale][:length]
        for i in range(1,scale):
            scale_i = scale_i + signal[i: len(signal) : scale][:length]
        scale_i = scale_i / scale
        #计算样本熵
        result.append(Sample_Entropy(scale_i, 2 ,r = 0.15))
        print("scale:" , scale, 'SampEn' , result[-1])
    return result


def P(x,y):
    x1 = np.sum(x) / len(x)
    y1 = np.sum(y) / len(y)
    p = np.sum((x-x1)*(y-y1)) / (np.sqrt(np.sum((x-x1)**2))*np.sqrt(np.sum((y-y1)**2)))
    return p

k1 = []
k2 = []
k3 = []
k4 = []
k5 = []
k6 = []
k7 = []
k8 = []
k9 = []
k10 = []
for i in range(400):
    data = data1[i]
    ceemdan = CEEMDAN()
    eIMFs = ceemdan.ceemdan(data)
    nIMFs = eIMFs.shape[0]
    Y = []
    for j in range(nIMFs):
        p = P(data,eIMFs[j])
        if p > 0.1:
           Y.append(eIMFs[j])
    y_c = np.zeros(2048)
    for k in range(len(Y)):
        y_c = y_c+Y[k]
    print(y_c)
    mse = MSE(y_c,10)
    k1.append(mse[0])
    k2.append(mse[1])
    k3.append(mse[2])
    k4.append(mse[3])
    k5.append(mse[4])
    k6.append(mse[5])
    k7.append(mse[6])
    k8.append(mse[7])
    k9.append(mse[8])
    k10.append(mse[9])

np.savetxt('E:/CEEMDAN-SSA-SVM/MSE1.txt',k1)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE2.txt',k2)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE3.txt',k3)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE4.txt',k4)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE5.txt',k5)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE6.txt',k6)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE7.txt',k7)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE8.txt',k8)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE9.txt',k9)
np.savetxt('E:/CEEMDAN-SSA-SVM/MSE10.txt',k10)

import pandas as pd
from sklearn.multiclass import OneVsRestClassifier
from sklearn.svm import SVC
import numpy as np
import pandas
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, recall_score
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import cross_val_score
import matplotlib.pyplot as plt
import matplotlib as mpl
import os
import glob
mpl.rcParams['font.sans-serif'] = ['KaiTi']   # 保证正常显示中文
mpl.rcParams['font.serif'] = ['KaiTi']        # 保证正常显示中文
mpl.rcParams['axes.unicode_minus'] = False    # 保证负号正常显示
import warnings
warnings.filterwarnings("ignore")

x = []
y = []
data = open('E:/CEEMDAN-SSA-SVM/tezheng.csv',encoding='UTF-8-sig')
for line in data.readlines():
    lineArr = line.strip().split(',')
    # data_train.append(lineArr)
    x.append(lineArr[1:11])
    y.append(lineArr[0])
x_train,x_test,y_train,y_test = train_test_split(x,y,test_size=0.3)
print(y_test)


model2 = SVC(kernel='rbf')
clf2 = model2.fit(x_train, y_train)

print("rbf径向基核函数-训练集：", model2.score(x_train, y_train))
print("rbf径向基核函数-测试集：", model2.score(x_test, y_test))

print(model2.decision_function(x_train))
# 使用predict()可以查看预测结果
print(model2.predict(x_train))


# 输出混淆矩阵
predictions = model2.predict(x_test)
cm = confusion_matrix(y_test, predictions)
print(cm)

# 5折交叉验证
score_SVM = cross_val_score(model2, x_test, y_test, cv=5)  # 五折交叉验证
print(score_SVM)
print('score均值', np.mean(score_SVM))  # 五折交叉验证的均值


# 绘制混淆矩阵
print('预测标签', predictions)
print('真实标签', y_test)
classes = list(map(int, y_test))
# 排序，准确对上分类结果
classes.sort()
# 对比，得到混淆矩阵
confusion = confusion_matrix(y_test, predictions)
# 热度图，后面是指定的颜色块，gray也可以，gray_x反色也可以
plt.imshow(confusion, cmap=plt.cm.Blues)
# 这个东西就要注意了
# ticks 这个是坐标轴上的坐标点
# label 这个是坐标轴的注释说明
indices = range(len(confusion))
# 坐标位置放入
# 第一个是迭代对象，表示坐标的顺序
# 第二个是坐标显示的数值的数组，第一个表示的其实就是坐标显示数字数组的index，但是记住必须是迭代对象
plt.xticks(indices, ['ball', 'inner', 'normal', 'outer'])
plt.yticks(indices, ['ball', 'inner', 'normal', 'outer'])
# 热度显示仪？就是旁边的那个验孕棒啦
plt.colorbar()
# 就是坐标轴含义说明了
plt.xlabel('y_test')
plt.ylabel('predictions')
plt.title('混淆矩阵')
# 显示数据，直观些
for first_index in range(len(confusion)):
    for second_index in range(len(confusion[first_index])):
        plt.text(first_index, second_index, confusion[first_index][second_index])

# 显示
plt.show()

CEEMDAN-多尺度熵（MSE）-SVM滚动轴承故障诊断

1、完全集合经验模态分解CEEMDAN本文就不做多介绍了，此算法是EEMD和EMD改进后得到的，一定程度上有效的抑制了模态混叠等问题。

加入经 EMD 分解后含辅助噪声的 IMF 分量，而不是将高斯白噪声信号直接添加在原始信号中；

2、多尺度熵

多尺度熵MSE的基本原理就是对时间序列进行粗粒化或下采样，即主要是在越来越粗略的时间分辨率下分析时间序列。

再附上样本熵代码

3.利用相关系数对CEEMDAN分解的信号进行重构。

P为计算得到的相关系数，当它大于0.1就保留下来。

4、选取实验数据损伤为0.018、转速1797r/min的西储大学轴承故障诊断数据，外圈故障、内圈故障、滚动体故障和正常状态分别取100个数据，每个数据长度为2048个点。

5、完整代码