机械故障诊断及工业工程故障诊断若干例子（第二篇）

　　MATLAB环境下基于平方包络谱的滚动轴承故障诊断
　　算法程序运行环境为MATLABR2018a，执行基于平方包络谱的滚动轴承故障诊断，也可用于金融时间序列，地震信号，机械振动信号，语音信号，声信号等一维时间序列信号。
　　Python环境下基于麦克风信号与随机森林的机器轴承运行状态识别
　　以往的轴承故障诊断都是围绕内外圈、滚动体等进行的，而该算法程序使用麦克风阵列信号和随机森林模型进行机器轴承运行状态识别，轴承运行状态包括正常状态（goodbearing），外圈单点损伤状态（badbearing），低速运行状态（slowspeedbearing），过度润滑状态（overlubrication），轴承润滑不足状态（underlubrication），轴承烧蚀故障状态（corona），轴承轻微变形（上拱，arching）等状态，看一下所使用的麦克风音频文件。
　　本项目要用到librosa模块，librosa是一个用于音乐和音频分析的python包，可以pipinstalllibrosa
　　首先导入相关模块importlibrosaimportlibrosa。displayfromscipy。ioimportwavfileaswavimportmatplotlib。pyplotaspltimportpandasaspdimportnumpyasnpfromsklearn。ensembleimportRandomForestClassifierpipinstallxgboostfromxgboostimportXGBClassifier
　　导入音频样本并部分可视化audiofileoverlubrication。wav轴承过度润滑状态data，sampleratelibrosa。load（audiofile）plt。figure（figsize（20，5））setsizeofvisualizationlibrosa。display。waveshow（data，srsamplerate）plt。xlabel（Time（samples））plt。ylabel（Amplitude）
　　audiofileunderlubrication。wav轴承润滑不足data，sampleratelibrosa。load（audiofile）plt。figure（figsize（20，5））setsizeofvisualizationlibrosa。display。waveshow（data，srsamplerate）plt。xlabel（Time（samples））plt。ylabel（Amplitude）
　　audiofilecorona。wav轴承烧蚀故障data，sampleratelibrosa。load（audiofile）plt。figure（figsize（20，5））setsizeofvisualizationlibrosa。display。waveshow（data，srsamplerate）plt。xlabel（Time（samples））plt。ylabel（Amplitude）
　　看一下某个样本音频的直方图采样频率print（Thesamplingfrequencyis：str（samplerate）Hz）音频样本的直方图plt。figure（figsize（20，5））plt。hist（data，bins100）plt。xlabel（Amplitude）plt。ylabel（Frequency）
　　大多数值接近0。
　　特征提取
　　对一个样本音频文件进行特征提取（梅尔频率倒谱系数MFCC），MFCC一般是用于表示声音信号的特征集mfccslibrosa。feature。mfcc（ydata，srsamplerate，nmfcc40）mfccs。shape有2371组不同的40个值的原因是因为音频样本的长度查看通过mfcc生成的特征mfccs〔0〕
　　对所有音频文件进行特征提取deffeatureextractor（file）：audio，srlibrosa。load（file）mfccfeatlibrosa。feature。mfcc（yaudio，srsr，nmfcc40）mfccscaledfeatnp。mean（mfccfeat。T，axis0）returnmfccscaledfeat输出标签labels〔arching，badbearing，goodbearing，corona，overlubrication，underlubrication，tracking，slowspeedbearing〕extractedfeat〔〕foriinlabels：filepathi。wavdatafeatureextractor（filepath）extractedfeat。append（data）
　　看一下特征长度len（extractedfeat）
　　显示提取特征的样本及其相应的标签print（str（str（labels〔0〕）：str（extractedfeat〔0〕）））
　　创建pandas数据框架，数据框架的的特征为列，标签为行data｛Features：extractedfeat，Class：labels｝dfpd。DataFrame（data）显示数据框架dfdf。convertdtypes（）df
　　初始化RandomForestClassifier分类器rfRandomForestClassifier（）
　　创建输入和输出数组，输入数组是特征，输出数组是标签Xdf〔Features〕。tolist（）ydf〔Class〕。tolist（）
　　训练模型rf。fit（X，y）
　　准确率rf。score（X，y）
　　准确率为100
　　MATLAB环境下基于双树复小波变换的轴承故障诊断
　　算法程序运行环境为MATLABR2021b，执行基于双树复小波变换的轴承故障诊断，也可用于金融时间序列，地震信号，机械振动信号，语音信号，声信号等一维时间序列信号。
　　双树复小波变换DTCWT为两个独立的两通道滤波器组，在实际应用用，不能随意选择两棵树中使用的尺度小波滤波器。第一棵树｛h0，h1｝的低通（尺度）和高通（小波）滤波器生成一个尺度函数和小波，另一棵树是由第一棵树的尺度函数的近似希尔伯特变换以及相应的高通滤波器生成的小波函数组成，记作｛g0，g1｝。因此由两棵树形成的复值尺度函数和小波函数是近似解析的，DTCWT的冗余度明显小于未抽取的DWT的冗余度。
　　双树复小波变换DTCWT基本理论
　　双树复小波变换DTCWT采用二叉树结构的两路滤波器组进行信号的分解和重构，第一棵树生成实部，第二棵生成虚部，合理设计实、虚部树低通滤波器，满足半采样延迟条件，具有近似平移不变性。两树的滤波器采样频率相同，但是它们之间的延迟恰好是一个采样间隔，这样虚部树中第层的二抽取恰好采到实部树中二抽取所丢掉的采样值，在获得了复小波变换的平移不变性的同时避免了大量的计算并且具有容易实现的优势。下图为层双树复小波的分解和重构过程。
　　虚线上方实部树变换的小波系数和尺度系数可由式（2）和（3）计算
　　相关的参考文献如下：
　　〔1〕HuangTongyuan，XuJia，YangYuling，HanBaoru。RobustZeroWatermarkingAlgorithmforMedicalImagesUsingDoubleTreeComplexWaveletTransformandHessenbergDecomposition〔J〕。Mathematics，2022，10（7）。
　　〔2〕ZhouYilu，FuXiaojin。ImageDenoisingBasedonDualtreeComplexWaveletTransformandConvolutionalNeuralNetwork〔J〕。JournalofPhysics：ConferenceSeries，2021，1995（1）。
　　〔3〕LeiWang，ZhiwenLiu，HongruiCao，XinZhang。Subbandaveragingkurtogramwithdualtreecomplexwaveletpackettransformforrotatingmachineryfaultdiagnosis〔J〕。MechanicalSystemsandSignalProcessing，2020，142（C）。
　　下面开始进行轴承故障诊断
　　首先以70Hz轴承外圈故障为例，故障特征频率因子为BPFI6。587，采样频率为10240HZ，首先看一下原始信号的包络谱
　　可见包络谱的故障特征频率虽然出现了，但是幅值相对较低。
　　下面进行DTCWT分解，分解层数为5
　　看一下实部树所对应的波形
　　实部树所对应的包络谱
　　重点看一下Level2的包络谱
　　故障特征频率的相对幅值较高
　　然后看一下虚部树所对应的分解波形
　　虚部树所对应的包络谱
　　然后看一下50Hz的情形，包络谱如下
　　实部树分解所得到的包络谱如下
　　详细看一下Level1的包络谱
　　可见轴承故障特征频率的相对幅值较高，虚部树与此类似
　　面包多代码
　　https：mbd。puboGeBENHAGEN