Auto_Diag/DiagnosisLib.py

"""
@Time    : 2024-03-04
@Author  : Wei
@FileName: DiagnosisLib.py
@Software: 故障诊断函数库

"""
#
import numpy as np
import matplotlib.pyplot as plt
from math import pi
from scipy import signal
from scipy import interpolate
import math

class DiagnosisLib:
    
    def __init__(self):
        pass

###############################################################################
###############################################################################
    # 半带滤波器系数
    def halfbandcoefficient(self):
        halfbandcoefficient = np.array([0.00000063557331607068,0,-0.00000627255234261115,0,0.00002474955105003003,0,\
                                        -0.00007140981582044673,0,0.00017237915446870680,0,-0.00036863887543283979,0,\
                                        0.00072039609678193576,0,-0.00131147951710543810,0,0.00225360507770744260,0,\
                                        -0.00369069002024503300,0,0.00580407535380912750,0,-0.00882079299392407940,0,\
                                        0.01302944110367636400,0,-0.01881320502983020000,0,0.02672090662754737100,0,\
                                        -0.03762627155576255700,0,0.05311294860790129200,0,-0.07653440360979835200,0,\
                                        0.11663131038150232000,0,-0.20562506738991029000,0,0.63439761518650584000,1,\
                                        0.63439761518650584000,0,-0.20562506738991029000,0,0.11663131038150232000,0,\
                                        -0.07653440360979835200,0,0.05311294860790129200,0,-0.03762627155576255700,0,\
                                        0.02672090662754737100,0,-0.01881320502983020000,0,0.01302944110367636400,0,\
                                        -0.00882079299392407940,0,0.00580407535380912750,0,-0.00369069002024503300,0,\
                                        0.00225360507770744260,0,-0.00131147951710543810,0,0.00072039609678193576,0,\
                                        -0.00036863887543283979,0,0.00017237915446870680,0,-0.00007140981582044673,0,\
                                        0.00002474955105003003,0,-0.00000627255234261115,0,0.00000063557331607068])
        return halfbandcoefficient

###############################################################################
###############################################################################
    # 计算时域波形均值
    def Average(self,xTimeWave):
        return sum(xTimeWave) / len(xTimeWave)

###############################################################################
###############################################################################
    # 计算时域波形xTimeWave的频谱
    # fs是时域波形的采样频率，单位是Hz
    # 调用代码须保证输入时域波形xTimeWave的数据点数为2的N幂
    # 频谱计算默认采用汉宁窗
    # 输出为频谱的x轴freq即频率，单位是Hz， 和频谱的y轴，单位是时域波形的幅值单位
    # x轴和y轴的数据点数为输入时域波形xTimeWave数据点数的1/2.56倍+1，即1024点时域波形
    # 输出频谱的x轴和y轴分别有401点数据，2048点时域波形时，输出频谱的x轴和y轴分别有801点数据。
    
    def Spectrum(self,xTimeWave,fs):
        DataLen = len(xTimeWave)
        xSpec = abs(np.fft.fft((xTimeWave*np.hanning(DataLen))))/DataLen*4
        freq = np.arange(0, fs / 2.56 + fs / DataLen, fs / DataLen)
        
        return [freq, xSpec[0:int(DataLen/2.56)+1]]
    
###############################################################################
###############################################################################
    # 计算时域波形xTimeWave的幅值谱倒频谱
    # fs是时域波形的采样频率，单位是Hz
    # 调用代码须保证输入时域波形xTimeWave的数据点数为2的N幂
    # 输出为倒谱的x轴quefrency,单位是秒，倒谱的幅值单位是dB
    def Cepstrum(self,xTimeWave,fs):
#        cepstrum = np.abs(np.fft.ifft(np.log(np.absolute(np.fft.fft(xTimeWave)))))
        cepstrum = 20*np.log((np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(xTimeWave))))))/0.1)
        quefrency = np.array(range(len(xTimeWave)))/fs
        
        return [quefrency,cepstrum]
    
###############################################################################
###############################################################################
    # 计算频谱xSpec的固定频率特征值,xSpec是频谱的y轴数组，其数据点数为101，201，401，801，1601，3201，6401，12801
    # fMax是频谱的分析频宽，单位为Hz
    # fTarget是特征频率的目标值，单位为Hz。比如轴的转频，齿轮啮合频率，水泵叶片通过频率等等
    # fRange是特征频率的搜索范围，其输入值根据fRangeMode的不同值而不一样，
    # fRangeMode可以赋值0或1，默认值是0。当fRangeMode = 0时，fRange的值是xx%， 比如 2%(实际输入为0.02) ，
    # 此时搜索范围是 fTarge +/- fTarge*fRange/2。 
    # 当fRangeMode = 1时，fRange的值是频率值，比如2Hz，
    # 此时的搜索范围就是 fTarge +/- fRange/2   
    # Detection是特征值的峰值（此时输入0），有效值（输入1），峰峰值（输入2）
    # numHarmonics为谐波数，0代表不需要考虑谐波即只有特征频率的基频，1代表包括基频和2倍频，2代表包括基频和2、3倍频
    def FixedFreqFeature(self,xSpec, fMax, fTarget, fRange, fRangeMode, Detection, numHarmonics):

        if fRangeMode == '百分比':
            fRange = fRange*fTarget
        if fRangeMode == '数值':
            fRange = fRange
            
#        print(fRange)
        SpecResolution = fMax/(len(xSpec)-1)

        index_SearchRangeMax = int(math.ceil((fTarget + fRange/2)/SpecResolution))
        index_SearchRangeMin = int(math.floor((fTarget - fRange/2)/SpecResolution))
        
        if index_SearchRangeMin < 3:
            index_SearchRangeMin = 3
            
        index_SearchRange = index_SearchRangeMax - index_SearchRangeMin
        
        xArray = xSpec[index_SearchRangeMin:(index_SearchRangeMax+1)]

        while((int(np.argmax(xArray))==0) or (int(np.argmax(xArray))==len(xArray)-1)):
            index_SearchRangeMin = index_SearchRangeMin -1
            index_SearchRange = index_SearchRange + 1
            xArray = xSpec[int(index_SearchRangeMin-math.floor(index_SearchRange/2)):int(index_SearchRangeMin+math.ceil(index_SearchRange/2+1))]


        ActualTargetFreq = self.FundamentalFreqActualValue(xArray, SpecResolution,index_SearchRangeMin)
        # print(ActualTargetFreq)
        
        featureValue = self.FreqFeatureValue(xSpec, SpecResolution, ActualTargetFreq, index_SearchRange, Detection, numHarmonics)
             
        return featureValue 
    
###############################################################################
############################################################################### 
    # 该函数只用于内部调用，其调用函数为 FixedFreqFeature      
    def FundamentalFreqActualValue(self,xArray, SpecResolution, index_SearchRangeMin):
        indexPeak = int(self.PeakFinderInArray(xArray))
        indexCenter = int(indexPeak + index_SearchRangeMin)
        ActualFundamentalFreq = self.FrequencyCompensation(indexCenter,xArray[indexPeak],xArray[indexPeak-1],xArray[indexPeak+1],SpecResolution)
        
        return ActualFundamentalFreq

###############################################################################
###############################################################################      
    # 该函数只用于内部调用，其调用函数为 FundamentalFreqActualValue    
    def PeakFinderInArray(self,xArray):
        NoElements = len(xArray)
        index_peak = int(np.argmin(xArray))                                     # initialize the peak index with the minimal value's index
        peak = xArray[index_peak]                                               # initialize the peak value with the minimal value
        
        for i in range(0, NoElements-2):
           if (xArray[i+1] > xArray[i]) and (xArray[i+1] > xArray[i+2]) and (xArray[i+1] > peak):
               index_peak = i+1
               peak = xArray[i+1]
        return index_peak
    
###############################################################################
###############################################################################
    # 该函数只用于内部调用，其调用函数为 FundamentalFreqActualValue    
    def FrequencyCompensation(self,indexCenter, AmpCenter, AmpLeft, AmpRight, xSpecResolution):
        if AmpRight >= AmpLeft:
            delta_k = (2*AmpRight-AmpCenter)/(AmpRight+AmpCenter)
        else:
            delta_k = (AmpCenter-2*AmpLeft)/(AmpLeft+AmpCenter)
        
        ActualFrequency = (indexCenter + delta_k)*xSpecResolution
        return ActualFrequency
    
###############################################################################
###############################################################################
    # 该函数只用于内部调用，其调用函数为 FixedFreqFeature    
    def FreqFeatureValue(self,xSpec, SpecResolution, ActualTargetFreq, index_SearchRange, Detection, numHarmonics = 0):
        sumation = 0
#        print(ActualTargetFreq)
        featureValue = []
        featureValue = [0 for i in range(numHarmonics+1)]
        for i in range(0,numHarmonics+1):
            index_i = int(round((i+1)*ActualTargetFreq/SpecResolution))

            a = xSpec[int(index_i-math.floor(index_SearchRange/2)):int(index_i+math.ceil(index_SearchRange/2+1))]
            while((int(np.argmax(a))==0) or (int(np.argmax(a))==len(a)-1)):
                index_SearchRange = index_SearchRange + 2
                a = xSpec[int(index_i-math.floor(index_SearchRange/2)):int(index_i+math.ceil(index_SearchRange/2+1))]

            indexPeak = int(self.PeakFinderInArray(a) + index_i-math.floor(index_SearchRange/2))
#            print(int(index_i-math.floor(index_SearchRange/2)),indexPeak,int(index_i+math.ceil(index_SearchRange/2)))
            for j in range(indexPeak-2,indexPeak+2+1):
                sumation = sumation + xSpec[j]*xSpec[j]
            
            featureValue[i] = np.sqrt(sumation/1.5)
            sumation = 0        
#            print(featureValue[i])
        
        if Detection == '峰值':
            pass
        if Detection == '有效值':
            featureValue = [x/math.sqrt(2) for x in featureValue]
        if Detection == '峰峰值':
            featureValue = [x*2 for x in featureValue]
            
        return featureValue

###############################################################################
###############################################################################
    # 计算频域积分
    # xSpecOrig需要对其计算积分的频谱数组，其长度为：101，201，401，801，1601，3201，6401，12801
    # fMax为频谱的分析频宽，单位为Hz
    # fCut为计算积分的低频截止频率，单位为Hz    
    def FreqIntegration(self,xSpecOrig, fMax, fCut):
        nLine = len(xSpecOrig) - 1
        lineResolution = fMax / nLine
        iStarting = int(math.floor(fCut / lineResolution))
        if iStarting == 0:
            iStarting = 1
        xSpecOut = np.zeros(len(xSpecOrig))
        for i in range(iStarting, len(xSpecOrig)):
            xSpecOut[i] = xSpecOrig[i]/(2*pi*lineResolution*i)
        return xSpecOut
    
###############################################################################    
###############################################################################
    # 巴特沃斯高通滤波
    # xTimeWave是时域波形数组， fs是采样频率，单位为Hz
    # fco 为高通截止频率，单位为Hz
    # Order为滤波器阶数，通常为偶数，2，4，6，8一般不超过8
    def HighPassButter(self,xTimeWave, fs, fco, Order):
        sos = signal.butter(Order, fco/(fs/2), 'hp', output='sos')
        filtered = signal.sosfilt(sos, xTimeWave)
        return filtered

###############################################################################
###############################################################################
    # 巴特沃斯低通滤波
    # xTimeWave是时域波形数组， fs是采样频率，单位为Hz
    # fco 为低通截止频率，单位为Hz
    # Order为滤波器阶数，通常为偶数，2，4，6，8一般不超过8     
    def LowPassButter(self,xTimeWave, fs, fco, Order):
        sos = signal.butter(Order, fco/(fs/2), 'lp', output='sos')
        filtered = signal.sosfilt(sos, xTimeWave)
        return filtered

###############################################################################
###############################################################################
    # 计算包络
    # xTimeWave是时域波形数组
    # DownSampleRate是降采样率，调用函数需提供此数值，其计算公式为
    #       原始波形xTimeWave的分析频宽/包络后时域波形的分析频宽，比如原始波形分析频宽是10000Hz(采样频率是25600Hz),包络后时域
    #       波形的分析频宽是1000Hz，那 DownSampleRate = 10000/1000 = 10
    # Method是SKF的包络方法，还是PeakHold包络方法      
    def EnvLopInTime(self,xTimeWave, DownSampleRate, Method = "SKF"):
        xOut = np.zeros(math.floor(len(xTimeWave)/DownSampleRate))
        if Method == "SKF":
            scale = 2.2
            xTimeWave = np.absolute(scale*xTimeWave)
            w = 1/(1.28*DownSampleRate)
            Order = 4
            sos = signal.butter(Order, w, 'lp', output='sos')
            xTimeWave = signal.sosfilt(sos, xTimeWave)
            xOut =xTimeWave[::DownSampleRate]
        if Method == "PeakHold":
            for i in range(0, len(xOut)):
                xOut[i] = np.amax(np.absolute(xTimeWave[i*DownSampleRate:((i+1)*DownSampleRate-1)]))
        return xOut - np.mean(xOut)

###############################################################################
###############################################################################
    # 内部调用函数
    def DownSampleBy2(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave, 1, 2)
        return xTimeWave

    def DownSampleBy4(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,4)
        return xTimeWave
    
    def DownSampleBy5(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,5)
        return xTimeWave
    
    def DownSampleBy8(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,8)
        return xTimeWave
    
    def DownSampleBy10(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        return xTimeWave
    
    def DownSampleBy20(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,2)
        return xTimeWave
    
    def DownSampleBy40(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,4)
        return xTimeWave
    
    def DownSampleBy50(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,5)
        return xTimeWave

    def DownSampleBy80(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,8)
        return xTimeWave

    def DownSampleBy100(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        return xTimeWave

    def DownSampleBy200(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,2)
        return xTimeWave
    
    def DownSampleBy400(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,4)
        return xTimeWave
    
    def DownSampleBy500(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,5)
        return xTimeWave

    def DownSampleBy800(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,8)
        return xTimeWave
    
    def DownSampleBy1000(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        return xTimeWave
    
    def DownSampleBy2000(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,2)
        return xTimeWave
    
    def DownSampleBy4000(self,xTimeWave):
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,10)
        xTimeWave = signal.resample_poly(xTimeWave,1,4)
        return xTimeWave
    

        
    def DownSample(self, xTimeWave, fsIn, fsOut, nPointIn, nPointOut):
#输入:	 xTimeWave:	时域波形数据
#	     Fs:	xTimeWave波形数据的采样频率(Hz)
#	     fs:	降采样后时域波形数据的采样频率(Hz)
#	     nPointIn:	输入波形的数据点数
#	     nPointOut:	返回的波形的数据点数
#
#返回：	返回一个长度为nPointOut的时域波形数组

        
        DownSampleRate = int(fsIn/fsOut)
        
        if int(nPointIn/nPointOut) < DownSampleRate:
            print('Input data length is not enough')
            return -1
        
        if DownSampleRate == 2:
            return self.DownSampleBy2(xTimeWave[0:int(nPointOut*2)])
        
        if DownSampleRate == 4:
            return self.DownSampleBy4(xTimeWave[0:int(nPointOut*4)])
        
        if DownSampleRate == 5:
            return self.DownSampleBy5(xTimeWave[0:int(nPointOut*5)])
        
        if DownSampleRate == 8:
            return self.DownSampleBy8(xTimeWave[0:int(nPointOut*8)])
        
        if DownSampleRate == 10:
            return self.DownSampleBy10(xTimeWave[0:int(nPointOut*10)])
        
        if DownSampleRate == 20:
            return self.DownSampleBy20(xTimeWave[0:int(nPointOut*20)])
        
        if DownSampleRate == 40:
            return self.DownSampleBy40(xTimeWave[0:int(nPointOut*40)])
        
        if DownSampleRate == 50:
            return self.DownSampleBy50(xTimeWave[0:int(nPointOut*50)])
        
        if DownSampleRate == 100:
            return self.DownSampleBy100(xTimeWave[0:int(nPointOut*100)])
        
        if DownSampleRate == 200:
            return self.DownSampleBy200(xTimeWave[0:int(nPointOut*200)])
        
        if DownSampleRate == 400:
            return self.DownSampleBy400(xTimeWave[0:int(nPointOut*400)])
        
        if DownSampleRate == 500:
            return self.DownSampleBy500(xTimeWave[0:int(nPointOut*500)])
        
        if DownSampleRate == 1000:
            return self.DownSampleBy1000(xTimeWave[0:int(nPointOut*1000)])
        
        if DownSampleRate == 2000:
            return self.DownSampleBy2000(xTimeWave[0:int(nPointOut*2000)])
        
        if DownSampleRate == 4000:
            return self.DownSampleBy4000(xTimeWave[0:int(nPointOut*4000)])
        
###############################################################################
###############################################################################

    def OrderAnalysis(self, xTimeWave, Fs, xTimeSpeedPulse, nPulsePerRev, nRev, nLine, TransLevel, TransLevelPara):
#输入:	 xTimeWave:	      时域波形数据
#	     Fs:	          xTimeWave波形数据的采样频率(Hz)
#	     xTimeSpeedPulse: 相邻转速脉冲信号的时间间隔（s）
#	     nPulsePerRev:	  每周的脉冲个数
#	     nRev:	          阶次分析整周数（整周数必须是2的N次幂，比如：4、8、16、32、64、128、256、512等等）
#        nLine：          阶次频谱的谱线数（谱线数可选项为：100、200、400、800、1600、3200、6400）
#        TransLevel：     阶次分析轴与转速测量轴的转速等级关系，
#                         0：表示阶次分析轴与转速分析轴是同一根轴
#                         1：表示阶次分析轴与转速分析轴不在同一根轴上，有1次转速比变化
#                         2：表示阶次分析轴与转速分析轴不在同一根轴上，有2次转速比变化
#                         3：表示阶次分析轴与转速分析轴不在同一根轴上，有3次转速比变化， 以此类推
#        TransLevelPara： 齿轮齿数数组，
#                         当  TransLevel == 0 时，这个数组只有1个元素，可以是[0]  
#                         当  TransLevel == 1 时，有一对齿轮啮合，2根轴。这个数组有2个元素，第一个元素是转速测量轴（1轴）上的齿轮齿数，
#                             第二个元素是与转速测量轴啮合（2轴）的齿轮的齿数，
#                         当  TransLevel == 2 时，有二对齿轮啮合，3根轴。这个数组有4个元素，第一个元素是转速测量轴（2轴）上的齿轮齿数，
#                             第二个元素是与转速测量轴啮合（2轴）的齿轮的齿数，第三个元素是2轴与下一级（3轴）啮合的齿轮的齿数，
#                             第四个元素是3轴与2轴啮合的齿轮的齿数。 
#返回：	等角度采样的角度轴坐标，等角度采样的波形，阶次谱的阶次轴坐标，阶次谱的幅值
        nRevSpeedSensorShaft = int(len(xTimeSpeedPulse)/nPulsePerRev)                                                           # 计算转速测量轴的转速测量信号有多少个整周数
        xTimeSpeedPulse = xTimeSpeedPulse[0:0+int(nRevSpeedSensorShaft*nPulsePerRev)]                                           # 整周数以外多余的脉冲信号去掉
        xTimeSpeedPerRevSpeedSensorShaft = np.sum(np.reshape(xTimeSpeedPulse,(nRevSpeedSensorShaft,nPulsePerRev)),axis = 1)     # 计算转速测量轴每转一周的时间
        
        if TransLevel == 0:
            xTimeSpeedPerRevCurrentShaft = xTimeSpeedPerRevSpeedSensorShaft  
            xTimeSpeedPerRevCurrentShaft_ave = np.average(xTimeSpeedPerRevCurrentShaft)                                                   # 当转速测量轴就是阶次分析的轴时，转速脉冲可以直接使用
            fMax_MAX = 1/xTimeSpeedPerRevCurrentShaft_ave*nLine/nRev
        else:          
            xTemp = np.cumsum(xTimeSpeedPerRevSpeedSensorShaft)                 # 转速脉冲时间累积
            xTemp = np.append(0,xTemp)                                          # 添加起始点                                                         # 当转速测量轴不是阶次分析轴时，传输脉冲信号必须转换到阶次分析轴上
            for i in range (TransLevel):
                up = TransLevelPara[i]                                          # 转速测量轴齿轮齿数
                down = TransLevelPara[i+1]                                      # 啮合轴齿轮齿数
                dataLength = int(len(xTemp)*up)                                 # 转速测量轴增采样后的数据长度
                flipelements = int(len(xTemp)*0.5)                              # 转速信号加长50%防止长度不够
                extention = np.zeros(flipelements)   
                for j in range(flipelements):                                   # 转速信号加长算法，以最后一点为对称中心
                    extention[j] = xTemp[-1]-(xTemp[-1-j-1]-xTemp[-1])  
                xTemp1 = np.append(xTemp,extention)                             # 转速信号加长后的转速脉冲时间累积
                xTemp2 = np.zeros(xTemp1.size*(up),xTemp1.dtype)
                xTemp2[::up] = xTemp1                                           # 脉冲时间累积信号增采样，增采样率为转速测量轴齿轮的齿数
                # plt.figure(30)
                # plt.plot(xTemp2) 
                #==============================================================
                delta1 = 0.00000001
                delta2 = 0.0000001
                freq1 = 0.39/up
                freq2 = 0.61/up
                AA = np.matrix([[-4.278e-1,-4.761e-1,0],[-5.941e-1,7.114e-2,0],[-2.66e-3,5.309e-03,0]])
                d1 = np.log10(delta1)
                d2 = np.log10(delta2)
                D = np.array([1,d1,d1*d1])*AA*np.array([[1],[d2],[d2*d2]])
                bb = np.array([11.01217,0.51244])
                fK = np.dot(np.array([1.0,(d1-d2)]),bb)
                df = np.abs(freq2-freq1)
                N = int(D[0]/df-fK*df+1)
                if np.remainder(N,2) == 0:
                    N = N+1
                upSampleFilterCoefficients = signal.remez(N, [0,0.39,0.61,up/2], [1,0],Hz = up,maxiter = 500)
                # [freq, response] = signal.freqz(upSampleFilterCoefficients)
                # ampl = np.abs(response)
                # plt.figure(31)
                # plt.plot(freq,ampl)
                #==============================================================
                groupDelay = int((len(upSampleFilterCoefficients)+1)/2)
                xTemp3 = signal.convolve(xTemp2,upSampleFilterCoefficients)
                xTemp4 = (xTemp3[groupDelay-1:groupDelay-1+dataLength])*up
                plt.figure(30)
                plt.plot(xTemp4)
                xTemp5 = xTemp4[0::down]
                
                
            xTimeSpeedPerRevCurrentShaft = np.diff(xTemp5) 
            xTimeSpeedPerRevCurrentShaft_ave = np.average(xTimeSpeedPerRevCurrentShaft)                                                   # 当转速测量轴就是阶次分析的轴时，转速脉冲可以直接使用
            fMax_MAX = 1/xTimeSpeedPerRevCurrentShaft_ave*nLine/nRev
        xTimeSpeedPerRevCurrentShaft_min = np.min(xTimeSpeedPerRevCurrentShaft)                                                 # 计算轴转得最快的一周的时间
        fMax_MAX = 1/xTimeSpeedPerRevCurrentShaft_min*nLine/nRev                                                                # 根据最快的一周的时间得出阶次频谱最高分析频宽的频率
        Order = 6
        if fMax_MAX <= Fs/2.56:
            xTimeWave = self.LowPassButter(xTimeWave, Fs, fMax_MAX, Order)      # 低通滤波消除大于fMax_MAX的信号，然后才能做阶次分析
  
        xTimeSpeedCurrentShaftCumulative = np.cumsum(xTimeSpeedPerRevCurrentShaft)                                              # 计算轴每转一周所用时间的累积值
        nVibDataPoint = nLine*2.56                                                                                              # 计算时域波形点数
        up = nVibDataPoint/nRev                                                                                                 # 计算每周等角度采样点数                                                             
#        down = 1
        halfbandWeight = self.halfbandcoefficient() 
        
        groupDelay = int((len(halfbandWeight)+1)/2)
        numIteration = int(np.log10(up)/np.log10(2))
        temp = np.append(0,xTimeSpeedCurrentShaftCumulative)
        flipelements = int(len(temp)*0.5)
        extention = np.zeros(flipelements)         
        for i in range(numIteration):
            dataLength = int(len(temp)*2)
            for j in range(flipelements):
                extention[j] = temp[-1]-(temp[-1-j-1]-temp[-1])  
            temp1 = np.append(temp,extention)
            temp2 = np.insert(temp1,np.arange(1,len(temp1),1),0)
            temp3 = np.append(temp2,0)
            temp4 = signal.convolve(temp3,halfbandWeight)
            temp = temp4[groupDelay-1:groupDelay-1+dataLength]

            
        xTimeNewSamplePoint = temp
   
        flipelements = int(len(xTimeWave)*0.25)
        extention = np.zeros(flipelements)
        for j in range(flipelements):
            extention[j] = xTimeWave[-1]-(xTimeWave[-1-j-1]-xTimeWave[-1])  
        xTimeWave = np.append(xTimeWave,extention)        
        xTimeSamplePoint = np.arange(0,len(xTimeWave)/Fs,1/Fs)                                                                  # 计算输入波形的采样时间点
#        spl = UnivariateSpline(xTimeSamplePoint,xTimeWave) 
#        Waveform = spl(xTimeNewSamplePoint)  
        fx = interpolate.interp1d(xTimeSamplePoint,xTimeWave,kind= 'cubic')
        Waveform = fx(xTimeNewSamplePoint)
        OrderWaveform = Waveform[0:int(nLine*2.56)]

        ct = 360/(nLine*2.56/nRev)
        orderTimeLine = np.arange(0,nLine*2.56*ct,ct)
        fOT = nLine*2.56/nRev
        fMax = fOT/2.56
        fCut = 10*fMax/nLine
        WinType = 0
        AverageType = 0
        OverLapType = 0
        [f, xSpectrum, xSpectrumPhase] = self.Spectrum(OrderWaveform,nLine,fMax,fCut,WinType, AverageType ,OverLapType)                                                     # 计算插值函数
        len(OrderWaveform)
        return [orderTimeLine, OrderWaveform, f, xSpectrum]    

###############################################################################
###############################################################################
    # 内部函数
    def halfBandDesign (self, filterLength, transitionBand ):
    
    	invalidInput = False
    
    	# check if integer
    	if (np.abs(filterLength - int(filterLength)) > 1e-10):
    		print('halfBandDesign.py: filterLength must be an integer')
    		invalidInput = True
    
    	# check if too small
    	if (filterLength < 7):
    		print('halfBandDesign.py: filterLength must be larger than 6')
    		invalidInput = True
    
    	# check if proper length
    	if (np.mod(filterLength+1,4) != 0):
    		print('halfBandDesign.py: filterLength+1 must be divisble by 4')
    		invalidInput = True
    
    	# check range for transition band
    	if (transitionBand <= 0 or transitionBand >= 0.5):
    		print('halfBandDesign.py: transitionBand must be greater than 0 and less than 0.5')
    		invalidInput = True
    
    	if (invalidInput):
    		return []
    
    	else:
    
    		# design a half band filter with remez
    		cutoff = 0.25
    		fPass = cutoff - (transitionBand/2)
    		fStop = cutoff + (transitionBand/2)
    		fVec = [0, fPass, fStop, 0.5]
    		aVec = [1, 0]
    
    		weights = signal.remez(filterLength,fVec,aVec)
    
    		# force zero weights
    		zeroWeightIndicesHalf = np.arange(2,(filterLength-1)/2,2,dtype=int)
    		zeroWeightIndicesNegative = np.concatenate((-zeroWeightIndicesHalf[::-1],zeroWeightIndicesHalf))
    		zeroWeightIndices = zeroWeightIndicesNegative - zeroWeightIndicesNegative[0] + 1
    
    		weights[zeroWeightIndices] = 0
    
    	return weights
###############################################################################
###############################################################################    
        
    # def Pk_Search_BW(self):
    #     BW = 0.03
    #     return BW
    
###############################################################################    
###############################################################################
    # @staticmethod    
    def gE_Bearing_Defect_HM_Threshhold(self,alert, danger, defectType):
        gE_Bearing_Defect_HM_alert = np.zeros(5)
        gE_Bearing_Defect_HM_danger = np.zeros(5)
        if defectType == 'BPFO' :
            gE_Bearing_Defect_HM_alert[0] = alert*0.167
            gE_Bearing_Defect_HM_alert[1] = alert*0.131
            gE_Bearing_Defect_HM_alert[2] = alert*0.098
            gE_Bearing_Defect_HM_alert[3] = alert*0.089
            gE_Bearing_Defect_HM_alert[4] = alert*0.082
            
            gE_Bearing_Defect_HM_danger[0] = danger*0.167
            gE_Bearing_Defect_HM_danger[1] = danger*0.131
            gE_Bearing_Defect_HM_danger[2] = danger*0.098
            gE_Bearing_Defect_HM_danger[3] = danger*0.089
            gE_Bearing_Defect_HM_danger[4] = danger*0.082
            
            return gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger
        
        if defectType == 'BPFI':
            gE_Bearing_Defect_HM_alert[0] = alert*0.088
            gE_Bearing_Defect_HM_alert[1] = alert*0.064
            gE_Bearing_Defect_HM_alert[2] = alert*0.056
            gE_Bearing_Defect_HM_alert[3] = alert*0.035
            gE_Bearing_Defect_HM_alert[4] = alert*0.024
            
            gE_Bearing_Defect_HM_danger[0] = danger*0.088
            gE_Bearing_Defect_HM_danger[1] = danger*0.064
            gE_Bearing_Defect_HM_danger[2] = danger*0.056
            gE_Bearing_Defect_HM_danger[3] = danger*0.035
            gE_Bearing_Defect_HM_danger[4] = danger*0.024
            
            return gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger
            
        if defectType == 'BSF':
            gE_Bearing_Defect_HM_alert[0] = alert*0.088
            gE_Bearing_Defect_HM_alert[1] = alert*0.064
            gE_Bearing_Defect_HM_alert[2] = alert*0.056
            gE_Bearing_Defect_HM_alert[3] = alert*0.035
            gE_Bearing_Defect_HM_alert[4] = alert*0.024
            
            gE_Bearing_Defect_HM_danger[0] = danger*0.088
            gE_Bearing_Defect_HM_danger[1] = danger*0.064
            gE_Bearing_Defect_HM_danger[2] = danger*0.056
            gE_Bearing_Defect_HM_danger[3] = danger*0.035
            gE_Bearing_Defect_HM_danger[4] = danger*0.024
            
        
            return gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger

        if defectType == 'FTF':
            gE_Bearing_Defect_HM_alert[0] = alert*0.088
            gE_Bearing_Defect_HM_alert[1] = alert*0.064
            gE_Bearing_Defect_HM_alert[2] = alert*0.056
            gE_Bearing_Defect_HM_alert[3] = alert*0.035
            gE_Bearing_Defect_HM_alert[4] = alert*0.024
            
            gE_Bearing_Defect_HM_danger[0] = danger*0.088
            gE_Bearing_Defect_HM_danger[1] = danger*0.064
            gE_Bearing_Defect_HM_danger[2] = danger*0.056
            gE_Bearing_Defect_HM_danger[3] = danger*0.035
            gE_Bearing_Defect_HM_danger[4] = danger*0.024
            
            return gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger    

###############################################################################
###############################################################################
    # @staticmethod    
    def Vel_Bearing_Defect_HM_Threshhold(self,alert, danger, defectType):
        Vel_Bearing_Defect_HM_alert = np.zeros(5)
        Vel_Bearing_Defect_HM_danger = np.zeros(5)
        if defectType == 'BPFO':
            Vel_Bearing_Defect_HM_alert[0] = alert*1.0
            Vel_Bearing_Defect_HM_alert[1] = alert*1.0
            Vel_Bearing_Defect_HM_alert[2] = alert*1.0
            Vel_Bearing_Defect_HM_alert[3] = alert*1.0

            
            Vel_Bearing_Defect_HM_danger[0] = danger*1.0
            Vel_Bearing_Defect_HM_danger[1] = danger*1.0
            Vel_Bearing_Defect_HM_danger[2] = danger*1.0
            Vel_Bearing_Defect_HM_danger[3] = danger*1.0
            
            return Vel_Bearing_Defect_HM_alert, Vel_Bearing_Defect_HM_danger
        
        if defectType == 'BPFI':
            Vel_Bearing_Defect_HM_alert[0] = alert*1.0
            Vel_Bearing_Defect_HM_alert[1] = alert*1.0
            Vel_Bearing_Defect_HM_alert[2] = alert*1.0
            Vel_Bearing_Defect_HM_alert[3] = alert*1.0

            
            Vel_Bearing_Defect_HM_danger[0] = danger*1.0
            Vel_Bearing_Defect_HM_danger[1] = danger*1.0
            Vel_Bearing_Defect_HM_danger[2] = danger*1.0
            Vel_Bearing_Defect_HM_danger[3] = danger*1.0
            
            return Vel_Bearing_Defect_HM_alert, Vel_Bearing_Defect_HM_danger
            
        if defectType == 'BSF':
            Vel_Bearing_Defect_HM_alert[0] = alert*1.0
            Vel_Bearing_Defect_HM_alert[1] = alert*0.5
            Vel_Bearing_Defect_HM_alert[2] = alert*0.5
            Vel_Bearing_Defect_HM_alert[3] = alert*0.5

            
            Vel_Bearing_Defect_HM_danger[0] = danger*1.0
            Vel_Bearing_Defect_HM_danger[1] = danger*0.5
            Vel_Bearing_Defect_HM_danger[2] = danger*0.5
            Vel_Bearing_Defect_HM_danger[3] = danger*0.5
            
            return Vel_Bearing_Defect_HM_alert, Vel_Bearing_Defect_HM_danger 
        
        
###############################################################################
###############################################################################
    # @staticmethod    
    def Vel_Bearing_RunningSpd_HM_Threshhold(self,alert, danger, defectType):
        Vel_Bearing_RunningSpd_HM_alert = np.zeros(5)
        Vel_Bearing_RunningSpd_HM_danger = np.zeros(5)

        if defectType == 'BSF':
            Vel_Bearing_RunningSpd_HM_alert[0] = alert*1.0
            Vel_Bearing_RunningSpd_HM_alert[1] = alert*0.5
            Vel_Bearing_RunningSpd_HM_alert[2] = alert*0.5
            Vel_Bearing_RunningSpd_HM_alert[3] = alert*0.5
            Vel_Bearing_RunningSpd_HM_alert[4] = alert*0.5
            
            Vel_Bearing_RunningSpd_HM_danger[0] = danger*1.0
            Vel_Bearing_RunningSpd_HM_danger[1] = danger*0.5
            Vel_Bearing_RunningSpd_HM_danger[2] = danger*0.5
            Vel_Bearing_RunningSpd_HM_danger[3] = danger*0.5
            Vel_Bearing_RunningSpd_HM_danger[4] = danger*0.5

            return Vel_Bearing_RunningSpd_HM_alert, Vel_Bearing_RunningSpd_HM_danger

###############################################################################
###############################################################################
    # @staticmethod                
    def Vel_Unbalance_HM_Threshold(self,alert, danger):
        Vel_Unbalance_HM_alert = np.zeros(5)
        Vel_Unbalance_HM_danger = np.zeros(5)

        Vel_Unbalance_HM_alert[0] = alert*1.0
        Vel_Unbalance_HM_alert[1] = alert*0.5
        Vel_Unbalance_HM_alert[2] = alert*0.5
        Vel_Unbalance_HM_alert[3] = alert*0.5
        Vel_Unbalance_HM_alert[4] = alert*0.5
        
        Vel_Unbalance_HM_danger[0] = danger*1.0
        Vel_Unbalance_HM_danger[1] = danger*0.5
        Vel_Unbalance_HM_danger[2] = danger*0.5
        Vel_Unbalance_HM_danger[3] = danger*0.5
        Vel_Unbalance_HM_danger[4] = danger*0.5
        
        return Vel_Unbalance_HM_alert, Vel_Unbalance_HM_danger

###############################################################################
###############################################################################
    # @staticmethod                
    def Vel_Misalignment_HM_Threshold(self,alert, danger):
        Vel_Misalignment_HM_alert = np.zeros(5)
        Vel_Misalignment_HM_danger = np.zeros(5)
        
        Vel_Misalignment_HM_alert[0] = alert*1.0
        Vel_Misalignment_HM_alert[1] = alert*0.5
        Vel_Misalignment_HM_alert[2] = alert*0.5
        Vel_Misalignment_HM_alert[3] = alert*0.5
        Vel_Misalignment_HM_alert[4] = alert*0.5
        
        Vel_Misalignment_HM_danger[0] = danger*1.0
        Vel_Misalignment_HM_danger[1] = danger*0.5
        Vel_Misalignment_HM_danger[2] = danger*0.5
        Vel_Misalignment_HM_danger[3] = danger*0.5
        Vel_Misalignment_HM_danger[4] = danger*0.5
        
        return Vel_Misalignment_HM_alert, Vel_Misalignment_HM_danger

###############################################################################
###############################################################################
    # @staticmethod                
    def Vel_Looseness_HM_Threshold(self,alert, danger):
        Vel_Looseness_HM_alert = np.zeros(5)
        Vel_Looseness_HM_danger = np.zeros(5)
        
        Vel_Looseness_HM_alert[0] = alert*1.0
        Vel_Looseness_HM_alert[1] = alert*0.5
        Vel_Looseness_HM_alert[2] = alert*0.5
        Vel_Looseness_HM_alert[3] = alert*0.5
        Vel_Looseness_HM_alert[4] = alert*0.5
        
        Vel_Looseness_HM_danger[0] = danger*1.0
        Vel_Looseness_HM_danger[1] = danger*0.5
        Vel_Looseness_HM_danger[2] = danger*0.5
        Vel_Looseness_HM_danger[3] = danger*0.5
        Vel_Looseness_HM_danger[4] = danger*0.5
        
        return Vel_Looseness_HM_alert, Vel_Looseness_HM_danger
    
###############################################################################
###############################################################################
    # @staticmethod    
    def gE_2XLF_HM_Threshhold(self,alert, danger):
        gE_2XLF_HM_alert = np.zeros(5)
        gE_2XLF_HM_danger = np.zeros(5)

        gE_2XLF_HM_alert[0] = alert*0.167
        gE_2XLF_HM_alert[1] = alert*0.131
        gE_2XLF_HM_alert[2] = alert*0.098
        gE_2XLF_HM_alert[3] = alert*0.089
        gE_2XLF_HM_alert[4] = alert*0.082
            
        gE_2XLF_HM_danger[0] = danger*0.167
        gE_2XLF_HM_danger[1] = danger*0.131
        gE_2XLF_HM_danger[2] = danger*0.098
        gE_2XLF_HM_danger[3] = danger*0.089
        gE_2XLF_HM_danger[4] = danger*0.082

        return gE_2XLF_HM_alert, gE_2XLF_HM_danger

###############################################################################
###############################################################################
    # @staticmethod    
    def Vel_2XLF_HM_Threshhold(self,alert, danger):
        Vel_2XLF_HM_alert = np.zeros(1)
        Vel_2XLF_HM_danger = np.zeros(1)

        Vel_2XLF_HM_alert[0] = alert*1.0          
        Vel_2XLF_HM_danger[0] = danger*1.0

        return Vel_2XLF_HM_alert, Vel_2XLF_HM_danger

###############################################################################
###############################################################################
    # @staticmethod            
    def Acc_GMF_HM_Threshhold(self,alert, danger):
        Acc_GMF_HM_alert = np.zeros(3)
        Acc_GMF_HM_danger = np.zeros(3)
        
        Acc_GMF_HM_alert[0] = alert*1.0
        Acc_GMF_HM_alert[1] = alert*1.0
        Acc_GMF_HM_alert[2] = alert*1.0
        
        Acc_GMF_HM_danger[0] = danger*1.0
        Acc_GMF_HM_danger[1] = danger*1.0
        Acc_GMF_HM_danger[2] = danger*1.0
        
        return Acc_GMF_HM_alert, Acc_GMF_HM_danger
        
###############################################################################
###############################################################################
    # @staticmethod            
    def Vel_GMF_HM_Threshhold(self,alert, danger):
        Vel_GMF_HM_alert = np.zeros(3)
        Vel_GMF_HM_danger = np.zeros(3)
        
        Vel_GMF_HM_alert[0] = alert*1.0
        Vel_GMF_HM_alert[1] = alert*1.0
        Vel_GMF_HM_alert[2] = alert*1.0
        
        Vel_GMF_HM_danger[0] = danger*1.0
        Vel_GMF_HM_danger[1] = danger*1.0
        Vel_GMF_HM_danger[2] = danger*1.0
        
        return Vel_GMF_HM_alert, Vel_GMF_HM_danger
    
###############################################################################
###############################################################################
    # @staticmethod            
    def gE_GMF_HM_Threshhold(self,alert, danger):
        gE_GMF_HM_alert = np.zeros(3)
        gE_GMF_HM_danger = np.zeros(3)
        
        gE_GMF_HM_alert[0] = alert*0.167
        gE_GMF_HM_alert[1] = alert*0.131
        gE_GMF_HM_alert[2] = alert*0.098
        
        gE_GMF_HM_danger[0] = danger*0.167
        gE_GMF_HM_danger[1] = danger*0.131
        gE_GMF_HM_danger[2] = danger*0.098
        
        return gE_GMF_HM_alert, gE_GMF_HM_danger
        
###############################################################################
###############################################################################
    # @staticmethod            
    def gE_GTIF_HM_Threshhold(self,alert, danger):
        gE_GTIF_HM_alert = np.zeros(5)
        gE_GTIF_HM_danger = np.zeros(5)
        
        gE_GTIF_HM_alert[0] = alert*0.167
        gE_GTIF_HM_alert[1] = alert*0.131
        gE_GTIF_HM_alert[2] = alert*0.098
        gE_GTIF_HM_alert[3] = alert*0.089
        gE_GTIF_HM_alert[4] = alert*0.082
        
        gE_GTIF_HM_danger[0] = danger*0.167
        gE_GTIF_HM_danger[1] = danger*0.131
        gE_GTIF_HM_danger[2] = danger*0.098
        gE_GTIF_HM_danger[3] = danger*0.089
        gE_GTIF_HM_danger[4] = danger*0.082
        
        return gE_GTIF_HM_alert, gE_GTIF_HM_danger
    
###############################################################################
###############################################################################   
    # @staticmethod     
    def Vel_VPF_HM_Threshhold(self,alert, danger):
        Vel_VPF_HM_alert = np.zeros(2)
        Vel_VPF_HM_danger = np.zeros(2)
        
        Vel_VPF_HM_alert[0] = alert*1.0
        Vel_VPF_HM_alert[1] = alert*1.0
        
        Vel_VPF_HM_danger[0] = danger*1.0
        Vel_VPF_HM_danger[1] = danger*1.0
        
        return Vel_VPF_HM_alert, Vel_VPF_HM_danger

###############################################################################
###############################################################################
    # @staticmethod            
    def Vel_BTF_HM_Threshhold(self,alert, danger):
        Vel_BTF_HM_alert = np.zeros(2)
        Vel_BTF_HM_danger = np.zeros(2)
        
        Vel_BTF_HM_alert[0] = alert*1.0
        Vel_BTF_HM_alert[1] = alert*1.0
        
        Vel_BTF_HM_danger[0] = danger*1.0
        Vel_BTF_HM_danger[1] = danger*1.0
        
        return Vel_BTF_HM_alert, Vel_BTF_HM_danger
    
###############################################################################
###############################################################################
    # @staticmethod            
    def Vel_SHF_HM_Threshhold(self,alert, danger):
        Vel_SHF_HM_alert = np.zeros(1)
        Vel_SHF_HM_danger = np.zeros(1)
        
        Vel_SHF_HM_alert[0] = alert*1.0
        
        Vel_SHF_HM_danger[0] = danger*1.0
        
        return Vel_SHF_HM_alert, Vel_SHF_HM_danger

###############################################################################
###############################################################################
    # @staticmethod            
    def Vel_RBPF_HM_Threshhold(self,alert, danger):
        Vel_RBPF_HM_alert = np.zeros(1)
        Vel_RBPF_HM_danger = np.zeros(1)
        
        Vel_RBPF_HM_alert[0] = alert*1.0
        
        Vel_RBPF_HM_danger[0] = danger*1.0
        
        return Vel_RBPF_HM_alert, Vel_RBPF_HM_danger
    
###############################################################################
###############################################################################
    # @staticmethod            
    def Acc_RBPF_HM_Threshhold(self,alert, danger):
        Acc_RBPF_HM_alert = np.zeros(1)
        Acc_RBPF_HM_danger = np.zeros(1)
        
        Acc_RBPF_HM_alert[0] = alert*1.0
        
        Acc_RBPF_HM_danger[0] = danger*1.0
        
        return Acc_RBPF_HM_alert, Acc_RBPF_HM_danger

###############################################################################
###############################################################################
    # @staticmethod    
    def Bearing_Fault_Diag(self,gE_Bearing_Defect_HM, Vel_Bearing_Defect_HM, Vel_RunningSpd_HM, gE_alert, gE_danger, Vel_alert, Vel_danger, defectType):
        BrDefDetection = False
        BrDefSev = 0
        gE_Bearing_Defect_HM_alert = np.zeros(5) 
        gE_Bearing_Defect_HM_danger = np.zeros(5)
        Vel_Bearing_Defect_HM_alert = np.zeros(5)
        Vel_Bearing_Defect_HM_danger = np.zeros(5)
        Vel_Bearing_RunningSpd_HM_alert = np.zeros(5)
        Vel_Bearing_RunningSpd_HM_danger = np.zeros(5)

###############################################################################
# Check bearing outer race and inner race defect
        
        if (defectType == 'BPFO') or (defectType == 'BPFI'):
            [gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger] = self.gE_Bearing_Defect_HM_Threshhold(gE_alert, gE_danger, defectType)
#            [Vel_Bearing_RunningSpd_HM_alert, Vel_Bearing_RunningSpd_HM_danger] = self.Vel_Bearing_RunningSpd_HM_Threshhold(Vel_alert, Vel_danger, defectType)
            [Vel_Bearing_Defect_HM_alert, Vel_Bearing_Defect_HM_danger] = self.Vel_Bearing_Defect_HM_Threshhold(Vel_alert, Vel_danger, defectType)

# Bearing outer race frequency, severe, confirmed by velocity spectrum, CF = 8, Severity = 10
# AND 
#    OR 
#       ge_BPFO_1 >= Alert
#       ge_BPFO_2 >= Alert
#    OR 
#       Vel_BPFO_1 >= Alert
#       Vel_BPFO_2 >= Alert
#
# above logics apply to BPFI as well
           
            if (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] or gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_alert[1]) and \
                (Vel_Bearing_Defect_HM[0] >= Vel_Bearing_Defect_HM_alert[0] or Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1]):
                    BrDefDetection = True
                    BrDefSev = 8
                    
                    return BrDefDetection, BrDefSev
                
# Bearing outer race frequency, severe, no confirming by velocity spectrum, CF = 6, Severity = 10
# AND 
#    OR 
#       AND 
#           ge_BPFO_1 >= Danger
#           ge_BPFO_2 >= Alert 
#       AND 
#           gE_BPFO_1 >= Danger
#           gE_BPFO_3 >= Alert
#       AND 
#           gE_BPFO_1 >= Danger
#           gE_BPFO_4 >= Alert
#   NOR
#       Vel_BPFO_1 >= Alert
#       Vel_BPFO_2 >= Alert
#
# above logics apply to BPFI as well
                       
            else:
                if ((gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] and gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_alert[1]) or \
                (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] and gE_Bearing_Defect_HM[2] >= gE_Bearing_Defect_HM_alert[2]) or \
                (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] and gE_Bearing_Defect_HM[3] >= gE_Bearing_Defect_HM_alert[3])) and \
                (not(Vel_Bearing_Defect_HM[0] >= Vel_Bearing_Defect_HM_alert[0] or Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1])):
                    BrDefDetection = True
                    BrDefSev = 6                
                    return BrDefDetection, BrDefSev
                
# Bearing outer race frequency, severe, not confirmed by gE spectrum, CF = 6, Severity = 10 
# OR 
#    Vel_BPFO_1 >= Alert
#    Vel_BPFO_2 >= Alert
#
# above logics apply to BPFI as well
                    
                else:
                    if (Vel_Bearing_Defect_HM[0] >= Vel_Bearing_Defect_HM_alert[0] or Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1]):
                        BrDefDetection = True
                        BrDefSev = 6   
                        
                        return BrDefDetection, BrDefSev  

# Bearing outer race frequency, moderate, several defect peaks, no confirming by velocity spectrum, CF = 5, Severity = 5
# AND 
#    OR 
#       AND 
#           ge_BPFO_1 >= Alert
#           ge_BPFO_2 >= Alert 
#       AND 
#           gE_BPFO_1 >= Alert
#           gE_BPFO_3 >= Alert
#       AND 
#           gE_BPFO_1 >= Alert
#           gE_BPFO_4 >= Alert
#   NOR
#       Vel_BPFO_1 >= Alert
#       Vel_BPFO_2 >= Alert
#
# above logics apply to BPFI as well
                        
                    else:
                        if ((gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] and gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_alert[1]) or \
                            (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] and gE_Bearing_Defect_HM[2] >= gE_Bearing_Defect_HM_alert[2]) or \
                            (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] and gE_Bearing_Defect_HM[3] >= gE_Bearing_Defect_HM_alert[3])) and \
                            (not(Vel_Bearing_Defect_HM[0] >= Vel_Bearing_Defect_HM_alert[0] or Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1])):
                            BrDefDetection = True
                            BrDefSev = 3                
                            return BrDefDetection, BrDefSev
                
# Bearing outer race frequency, suspect, only 1x defect peak evaluated, CF = 3, Severity = 5
# ge_BPFO_1 >= Alert
#
# above logics apply to BPFI as well
 
                        else:
                            if (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0]):
                                BrDefDetection = True
                                BrDefSev = 1                
                                return BrDefDetection, BrDefSev
                
            return BrDefDetection, BrDefSev
###############################################################################                
# Check bearing rolling element defect
                    

        if defectType == 'BSF':
            [gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger] = self.gE_Bearing_Defect_HM_Threshhold(gE_alert, gE_danger, defectType)
            [Vel_Bearing_RunningSpd_HM_alert, Vel_Bearing_RunningSpd_HM_danger] = self.Vel_Bearing_RunningSpd_HM_Threshhold(Vel_alert, Vel_danger, defectType)
            [Vel_Bearing_Defect_HM_alert, Vel_Bearing_Defect_HM_danger] = self.Vel_Bearing_Defect_HM_Threshhold(Vel_alert, Vel_danger, defectType)     
#            print(Vel_Bearing_RunningSpd_HM_alert, Vel_Bearing_RunningSpd_HM_danger)
# Bearing rolling element frequency, severe, CF = 9, Severity = 10
# AND
#    AND
#        Vel_BSF_2 >= Alert
#        Vel_BSF_4 >= Alert
#    NOR
#        Vel_RS_2 >= Alert
#        Vel_RS_3 >= Alert
#        Vel_RS_4 >= Alert
#        Vel_RS_5 >= Alert
            
            if (Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1] and Vel_Bearing_Defect_HM[3] >= Vel_Bearing_Defect_HM_alert[3]) and \
                (not(Vel_RunningSpd_HM[1] >= Vel_Bearing_RunningSpd_HM_alert[1] or Vel_RunningSpd_HM[2] >= Vel_Bearing_RunningSpd_HM_alert[2] or \
                 Vel_RunningSpd_HM[3] >= Vel_Bearing_RunningSpd_HM_alert[3] or Vel_RunningSpd_HM[4] >= Vel_Bearing_RunningSpd_HM_alert[4])):
                    BrDefDetection = True
                    BrDefSev = 9                
                    return BrDefDetection, BrDefSev
 
# Bearing rolling element frequency, moderate, CF = 9, Severity = 5
# AND
#    AND
#        gE_BSF_1 >= Alert
#    NOR
#        Vel_BSF_2 >= Alert
#        Vel_BSF_4 >= Alert
#        Vel_RS_2 >= Alert
#        Vel_RS_3 >= Alert
#        Vel_RS_4 >= Alert
#        Vel_RS_5 >= Alert           
                
            else:
                if gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0]  and \
                    (not(Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1] or Vel_Bearing_Defect_HM[3] >= Vel_Bearing_Defect_HM_alert[3] or \
                    Vel_RunningSpd_HM[1] >= Vel_Bearing_RunningSpd_HM_alert[1] or Vel_RunningSpd_HM[2] >= Vel_Bearing_RunningSpd_HM_alert[2] or \
                    Vel_RunningSpd_HM[3] >= Vel_Bearing_RunningSpd_HM_alert[3] or Vel_RunningSpd_HM[4] >= Vel_Bearing_RunningSpd_HM_alert[4])):
                    BrDefDetection = True
                    BrDefSev = 5                
                    return BrDefDetection, BrDefSev

# Bearing rolling element frequency, suspect, CF = 3, Severity = 5
# AND
#    OR
#        gE_BSF_1 = Alert
#        gE_BSF_2 = Alert
#        gE_BSF_4 = Alert
#    NOR
#        gE_BSF_1 = Danger
#        gE_BSF_2 = Danger
#        gE_BSF_4 = Danger
#        Vel_BSF_2 >= Alert
#        Vel_BSF_4 >= Alert
#        Vel_RS_2 >= Alert
#        Vel_RS_3 >= Alert
#        Vel_RS_4 >= Alert
#        Vel_RS_5 >= Alert 
                
                else:
                    if ((gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] and gE_Bearing_Defect_HM[0] < gE_Bearing_Defect_HM_danger[0]) or \
                        (gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_alert[1] and gE_Bearing_Defect_HM[1] < gE_Bearing_Defect_HM_danger[1]) or \
                        (gE_Bearing_Defect_HM[3] >= gE_Bearing_Defect_HM_alert[3] and gE_Bearing_Defect_HM[3] < gE_Bearing_Defect_HM_danger[3] )) and \
                        (not(gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] or gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_danger[1] or gE_Bearing_Defect_HM[3] >= gE_Bearing_Defect_HM_danger[3] or \
                        Vel_Bearing_Defect_HM[0] >= Vel_Bearing_Defect_HM_alert[0] or Vel_Bearing_Defect_HM[1] >= Vel_Bearing_Defect_HM_alert[1] or \
                        Vel_RunningSpd_HM[1] >= Vel_Bearing_RunningSpd_HM_alert[1] or Vel_RunningSpd_HM[2] >= Vel_Bearing_RunningSpd_HM_alert[2] or \
                        Vel_RunningSpd_HM[3] >= Vel_Bearing_RunningSpd_HM_alert[3] or Vel_RunningSpd_HM[4] >= Vel_Bearing_RunningSpd_HM_alert[4])):
                        BrDefDetection = True
                        BrDefSev = 1                
                    
                        return BrDefDetection, BrDefSev
                
            return BrDefDetection, BrDefSev
###############################################################################                
# Check bearing cage defect

        if defectType == 'FTF':
            [gE_Bearing_Defect_HM_alert, gE_Bearing_Defect_HM_danger] = self.gE_Bearing_Defect_HM_Threshhold(gE_alert, gE_danger, defectType)

# Bearing cage frequency, severe, CF = 6, Severity = 10
# OR
#    gE_CG_1 >= Danger
#    gE_CG_2 >= Danger     
            
            if gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] or gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_danger[1]:
                    BrDefDetection = True
                    BrDefSev = 6
                    
                    return BrDefDetection, BrDefSev

# Bearing cage frequency, moderate, CF = 5, Severity = 5
# AND
#    OR
#       gE_CG_1 >= Alert
#       gE_CG_2 >= Alert     
#    NOR
#       gE_CG_1 = Danger
#       gE_CG_2 = Danger           
                    
            else:
                if (gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_alert[0] or gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_alert[1]) and \
                    (not(gE_Bearing_Defect_HM[0] >= gE_Bearing_Defect_HM_danger[0] or gE_Bearing_Defect_HM[1] >= gE_Bearing_Defect_HM_danger[1])):
                    BrDefDetection = True
                    BrDefSev = 3
                    
                    return BrDefDetection, BrDefSev
            
            return BrDefDetection, BrDefSev
        
###############################################################################                
###############################################################################
    # @staticmethod    
    def RunningSpd_Symptom_Rule(self,Vel_RunningSpd_HM, Vel_RunningSpd_HM_alert, Vel_RunningSpd_HM_danger):
        RPM1_MOD = False
        RPM1_SEV = False
        RPM2_MOD = False
        RPM2_SEV = False
        RPM345_MOD = False
        RPM345_SEV = False

# 1x RPM, Moderate, CF=8, Severity = 5
# AND
#     Vel_RS_1 = Alert  
            
        if Vel_RunningSpd_HM[0] >= Vel_RunningSpd_HM_alert[0] and Vel_RunningSpd_HM[0] < Vel_RunningSpd_HM_danger[0]:
            RPM1_MOD = True

# 1x RPM, Severe, CF=8, Severity = 10
# AND
#     Vel_RS_1 = Danger 
                
        if Vel_RunningSpd_HM[0] >= Vel_RunningSpd_HM_danger[0]:
            RPM1_SEV = True
                
# 2x RPM, Moderate, CF=8, Severity = 5
# AND
#     Vel_RS_2 = Alert 
            
        if Vel_RunningSpd_HM[1] >= Vel_RunningSpd_HM_alert[1] and Vel_RunningSpd_HM[1] < Vel_RunningSpd_HM_danger[1]:
            RPM2_MOD = True           

# 2x RPM, Severe, CF=8, Severity = 10
# AND
#     Vel_RS_2 = Danger 
                
        if Vel_RunningSpd_HM[1] >= Vel_RunningSpd_HM_danger[1]:
            RPM2_SEV = True

# 3x, 4x, 5x RPM, moderate, CF=8, Severity = 5
# AND
#    OR
#       Vel_RS_3 = Alert
#       Vel_RS_4 = Alert
#       Vel_RS_5 = Alert 
#    NOR
#       Vel_RS_3 = Danger
#       Vel_RS_4 = Danger
#       Vel_RS_5 = Danger 
        
        if (Vel_RunningSpd_HM[2] >= Vel_RunningSpd_HM_alert[2] and Vel_RunningSpd_HM[2] < Vel_RunningSpd_HM_danger[2] or \
            Vel_RunningSpd_HM[3] >= Vel_RunningSpd_HM_alert[3] and Vel_RunningSpd_HM[3] < Vel_RunningSpd_HM_danger[3] or \
            Vel_RunningSpd_HM[4] >= Vel_RunningSpd_HM_alert[4] and Vel_RunningSpd_HM[4] < Vel_RunningSpd_HM_danger[4]) and \
            (not(Vel_RunningSpd_HM[2] >= Vel_RunningSpd_HM_danger[2] or Vel_RunningSpd_HM[3] >= Vel_RunningSpd_HM_danger[3] or \
             Vel_RunningSpd_HM[4] >= Vel_RunningSpd_HM_danger[4])):
            RPM345_MOD = True
        
# 3x, 4x, 5x RPM, Sever, CF=8, Severity = 10
# OR
#     Vel_RS_3 = Danger
#     Vel_RS_4 = Danger
#     Vel_RS_5 = Danger   

        if Vel_RunningSpd_HM[2] >=  Vel_RunningSpd_HM_danger[2] or Vel_RunningSpd_HM[3] >=  Vel_RunningSpd_HM_danger[3] or \
            Vel_RunningSpd_HM[4] >=  Vel_RunningSpd_HM_danger[4]:
            RPM345_SEV = True
        
        return RPM1_MOD, RPM1_SEV, RPM2_MOD, RPM2_SEV, RPM345_MOD, RPM345_SEV

###############################################################################                
###############################################################################
    # @staticmethod    
# Note: The alert and danger threshold values (in mm/s RMS) are set according to ISO vibration standard for flexible mounting in group 3
                    
    def RunningSpd_Fault_Diag(self,Vel_RunningSpd_HM, Vel_alert = 7.1, Vel_danger = 11, defectType = 'Unbalance'):
        RPM1_MOD = False
        RPM1_SEV = False
        RPM2_MOD = False
        RPM2_SEV = False
        RPM345_MOD = False
        RPM345_SEV = False
        
        UnbalanceDetection = False
        UnbalanceSev = 0
        MisalignmentDetction = False
        MisalignmentSev = 0
        LoosenessDetection = False
        LoosenessSev = 0
        
        Vel_Unbalance_HM_alert = np.zeros(5) 
        Vel_Unbalance_HM_danger = np.zeros(5)
        Vel_Misalignment_HM_alert = np.zeros(5) 
        Vel_Misalignment_HM_danger = np.zeros(5)
        Vel_Looseness_HM_alert = np.zeros(5) 
        Vel_Looseness_HM_danger = np.zeros(5)
        
############################################################################### 
# Check Unbalance             
        if defectType == 'Unbalance':
            [Vel_Unbalance_HM_alert, Vel_Unbalance_HM_danger] = self.Vel_Unbalance_HM_Threshold(Vel_alert, Vel_danger)
            [RPM1_MOD, RPM1_SEV, RPM2_MOD, RPM2_SEV, RPM345_MOD, RPM345_SEV] = self.RunningSpd_Symptom_Rule(Vel_RunningSpd_HM, Vel_Unbalance_HM_alert, Vel_Unbalance_HM_danger)
            # print('RPM1_MOD,RPM1_SEV, RPM2_MOD, RPM2_SEV, RPM345_MOD, RPM345_SEV','   ',RPM1_MOD,'   ', RPM1_SEV,'   ', RPM2_MOD, '   ',RPM2_SEV, '   ',RPM345_MOD, '   ',RPM345_SEV)

# Unbalance, severe, rather confident, 1x RPM peak dominates, CF=8, Severity = 10
# AND
#    RPM1_SEV == True
#    NOR
#       RPM2_MOD == True                
#       RPM2_SEV == True
#       RPM345_MOD == True
#       RPM345_SEV == True                
                
            if RPM1_SEV == True and (not (RPM2_MOD == True or RPM2_SEV == True or RPM345_MOD == True or RPM345_SEV == True )):
                UnbalanceDetection = True
                UnbalanceSev = 8

            else:

# Unbalance, severe, somewhat confident 1x RPM peak dominates but low level of 2xRPM is possible, CF=6, Severity = 10
# AND
#    RPM1_SEV == True
#    NOR
#       RPM2_SEV == True
#       RPM345_MOD == True
#       RPM345_SEV == True
                
                if RPM1_SEV == True and (not(RPM2_SEV == True or RPM345_MOD == True or RPM345_SEV == True )):
                    UnbalanceDetection = True
                    UnbalanceSev = 6                


                else:
# Unbalance, moderate, rather confident since 1x RPM peak dominates, CF=8, Severity = 5
# AND
#    RPM1_MOD == True
#    NOR
#       RPM2_MOD == True
#       RPM2_SEV == True
#       RPM345_MOD == True
#       RPM345_SEV == True
            
                    if RPM1_MOD == True and (not(RPM2_MOD == True or RPM2_SEV == True or RPM345_MOD == True or RPM345_SEV == True )):
                        UnbalanceDetection = True
                        UnbalanceSev = 5

            return UnbalanceDetection, UnbalanceSev
###############################################################################  
# Check Misalignment

        if defectType == 'Misalignment':        
            [Vel_Misalignment_HM_alert, Vel_Misalignment_HM_danger] = self.Vel_Misalignment_HM_Threshold(Vel_alert, Vel_danger)
            [RPM1_MOD, RPM1_SEV, RPM2_MOD, RPM2_SEV, RPM345_MOD, RPM345_SEV] = self.RunningSpd_Symptom_Rule(Vel_RunningSpd_HM, Vel_Misalignment_HM_alert, Vel_Misalignment_HM_danger)

# Misalignment, moderate  CF = 8, Severity = 10
# AND
#    RPM2_SEV == True
#    OR
#      RPM1_MOD == True
#      RPM1_SEV == True
#      RPM345_MOD == True
#   NOR     
#      RPM345_SEV == True             
           
            if RPM2_SEV == True and (RPM1_MOD == True or RPM1_SEV == True or RPM345_MOD == True) and (not(RPM1_SEV == True)):
                MisalignmentDetction = True
                MisalignmentSev = 8

# Misalignment, moderate  CF = 6, Severity 10
# AND
#    RPM2_SEV == True
#    NOR
#       RPM1_MOD == True
#       RPM1_SEV == True
#       RPM345_MOD == True
#       RPM345_SEV == True             

            else: 
                if RPM2_SEV == True and (not(RPM1_MOD == True or RPM1_SEV == True or RPM345_MOD == True or RPM345_SEV == True)):
                    MisalignmentDetction = True
                    MisalignmentSev = 6

# Misalignment, moderate  CF = 8, Severity = 5
# AND
#    RPM2_MOD == True
#    OR
#      RPM1_MOD == True
#      RPM345_MOD == True
#   NOR
#      RPM1_SEV == True
#      RPM345_SEV == True        
                else:
                    if RPM2_MOD == True and (RPM1_MOD == True or RPM345_MOD == True) and (not(RPM1_SEV == True or RPM345_SEV == True)):
                        MisalignmentDetction = True
                        MisalignmentSev = 5

# Misalignment, moderate  CF = 6, Severity = 5
# AND
#    RPM2_MOD == True
#    NOR
#       RPM1_MOD == True
#       RPM1_SEV == True
#       RPM345_MOD == True
#       RPM345_SEV == True     
            
                    else:
                        if RPM2_MOD == True and (not (RPM1_MOD == True or RPM1_SEV == True or RPM345_MOD == True or RPM345_SEV == True)):
                            MisalignmentDetction = True
                            MisalignmentSev = 4
        
            return MisalignmentDetction, MisalignmentSev
###############################################################################
# Check Looseness
        if defectType == 'Looseness':          
            [Vel_Looseness_HM_alert, Vel_Looseness_HM_danger] = self.Vel_Looseness_HM_Threshold(Vel_alert, Vel_danger)
            [RPM1_MOD, RPM1_SEV, RPM2_MOD, RPM2_SEV, RPM345_MOD, RPM345_SEV] = self.RunningSpd_Symptom_Rule(Vel_RunningSpd_HM, Vel_Looseness_HM_alert, Vel_Looseness_HM_danger)        

# Mechanical looseness, sever  CF = 8, Severity = 10
# OR
#   AND
#        RPM1_SEV == True           
#       OR
#          RPM345_MOD == True
#          RPM345_SEV == True
#   AND
#        RPM2_SEV == True           
#       OR
#          RPM345_MOD == True
#          RPM345_SEV == True  
#   AND
#        RPM345_SEV == True           
#       OR
#          RPM1_MOD == True
#          RPM1_SEV == True
#   AND
#        RPM345_SEV == True           
#       OR
#          RPM2_MOD == True
#          RPM2_SEV == True
 
            if (RPM1_SEV == True and (RPM345_MOD == True or RPM345_SEV == True )) or (RPM2_SEV == True and (RPM345_MOD == True or RPM345_SEV == True )) or \
                (RPM345_SEV == True and (RPM1_MOD == True or RPM1_SEV == True )) or (RPM345_SEV == True and (RPM2_MOD == True or RPM2_SEV == True )):
                    LoosenessDetection = True
                    LoosenessSev = 8   

# Mechanical looseness, sever  CF = 6, Severity = 10
# AND
#    RPM345_SEV == True    
#    NOR
#       RPM1_MOD == True
#       RPM1_SEV == True
#       RPM2_MOD == True
#       RPM2_SEV == True
            
            else:
                if RPM345_SEV == True and (not( RPM1_MOD == True or RPM1_SEV == True or RPM2_MOD == True or RPM2_SEV == True )):
                    LoosenessDetection = True
                    LoosenessSev = 6

# Mechanical looseness, moderate  CF = 8, Severity = 5
# AND
#    OR
#       AND 
#           RPM1_MOD == True
#           RPM345_MOD == True
#       AND 
#           RPM2_MOD == True
#           RPM345_MOD == True
#    NOR
#       RPM1_SEV == True
#       RPM2_SEV == True   
#       RPM345_SEV == True     

                else:
                    if ((RPM1_MOD == True and RPM345_MOD == True )or(RPM2_MOD == True and RPM345_MOD == True))and\
                        (not(RPM1_SEV == True or RPM2_SEV == True or RPM345_SEV == True)):
                        LoosenessDetection = True
                        LoosenessSev = 5
  
# Mechanical looseness, moderate  CF = 6, Severity = 5
# AND
#    RPM345_MOD == True    
#    NOR
#       RPM1_MOD == True
#       RPM1_SEV == True
#       RPM2_MOD == True
#       RPM2_SEV == True
            
                    else:
                        if RPM345_MOD == True and (not( RPM1_MOD == True or RPM1_SEV == True or RPM2_MOD == True or \
                                                       RPM2_SEV == True )):
                            LoosenessDetection = True
                            LoosenessSev = 4

    
            return LoosenessDetection, LoosenessSev
        
###############################################################################                
###############################################################################
    # @staticmethod    
    def GMF_GTIF_Symptom_Rule(self,Vel_GMF_HM, Acc_GMF_HM, gE_GMF_HM, gE_GTIF_HM, Vel_GMF_HM_alert, Vel_GMF_HM_danger, \
                              Acc_GMF_HM_alert, Acc_GMF_HM_danger, gE_GMF_HM_alert, gE_GMF_HM_danger, gE_GTIF_HM_alert, \
                                  gE_GTIF_HM_danger):
        GMF_MOD = False
        GMF_SEV = False
        GTIF = False

        
# Gear mesh frequency and harmonics, moderate, CF=8, Severity = 5
# AND
#    OR
#       Vel_GMF_1 = Alert
#       Vel_GMF_2 = Alert
#       Vel_GMF_3 = Alert 
#       Acc_GMF_1 = Alert
#       Acc_GMF_2 = Alert
#       Acc_GMF_3 = Alert 
#       gE_GMF_1 = Alert
#       gE_GMF_2 = Alert
#       gE_GMF_3 = Alert         
#    NOR
#       Vel_GMF_1 = Danger
#       Vel_GMF_2 = Danger
#       Vel_GMF_3 = Danger 
#       Acc_GMF_1 = Danger
#       Acc_GMF_2 = Danger
#       Acc_GMF_3 = Danger
#       gE_GMF_1 = Danger
#       gE_GMF_2 = Danger
#       gE_GMF_3 = Danger
            
        if (Vel_GMF_HM[0] >= Vel_GMF_HM_alert[0] and Vel_GMF_HM[0] < Vel_GMF_HM_danger[0] or \
            Vel_GMF_HM[1] >= Vel_GMF_HM_alert[1] and Vel_GMF_HM[1] < Vel_GMF_HM_danger[1] or \
            Vel_GMF_HM[2] >= Vel_GMF_HM_alert[2] and Vel_GMF_HM[2] < Vel_GMF_HM_danger[2] or \
            Acc_GMF_HM[0] >= Acc_GMF_HM_alert[0] and Acc_GMF_HM[0] < Acc_GMF_HM_danger[0] or \
            Acc_GMF_HM[1] >= Acc_GMF_HM_alert[1] and Acc_GMF_HM[1] < Acc_GMF_HM_danger[1] or \
            Acc_GMF_HM[2] >= Acc_GMF_HM_alert[2] and Acc_GMF_HM[2] < Acc_GMF_HM_danger[2] or \
            gE_GMF_HM[0] >= gE_GMF_HM_alert[0] and gE_GMF_HM[0] < gE_GMF_HM_danger[0] or \
            gE_GMF_HM[1] >= gE_GMF_HM_alert[1] and gE_GMF_HM[1] < gE_GMF_HM_danger[1] or \
            gE_GMF_HM[2] >= gE_GMF_HM_alert[2] and gE_GMF_HM[2] < gE_GMF_HM_danger[2]) \
            and \
            (not(Vel_GMF_HM[0] >= Vel_GMF_HM_danger[0] or \
                  Vel_GMF_HM[1] >= Vel_GMF_HM_danger[1] or \
                  Vel_GMF_HM[2] >= Vel_GMF_HM_danger[2] or \
                  Acc_GMF_HM[0] >= Acc_GMF_HM_danger[0] or \
                  Acc_GMF_HM[1] >= Acc_GMF_HM_danger[1] or \
                  Acc_GMF_HM[2] >= Acc_GMF_HM_danger[2] or \
                  gE_GMF_HM[0] >= gE_GMF_HM_danger[0] or \
                  gE_GMF_HM[1] >= gE_GMF_HM_danger[1] or \
                  gE_GMF_HM[2] >= gE_GMF_HM_danger[2])):
            
            GMF_MOD = True

# Gear mesh frequency and harmonics, sever, CF=8, Severity = 10
# OR
#     Vel_GMF_1 = Danger
#     Vel_GMF_2 = Danger
#     Vel_GMF_3 = Danger
#     Acc_GMF_1 = Danger
#     Acc_GMF_2 = Danger
#     Acc_GMF_3 = Danger 
#     gE_GMF_1 = Danger
#     gE_GMF_2 = Danger
#     gE_GMF_3 = Danger         

            
        if (Vel_GMF_HM[0] >= Vel_GMF_HM_danger[0] or \
            Vel_GMF_HM[1] >= Vel_GMF_HM_danger[1] or \
            Vel_GMF_HM[2] >= Vel_GMF_HM_danger[2] or \
            Acc_GMF_HM[0] >= Acc_GMF_HM_danger[0] or \
            Acc_GMF_HM[1] >= Acc_GMF_HM_danger[1] or \
            Acc_GMF_HM[2] >= Acc_GMF_HM_danger[2] or \
            gE_GMF_HM[0] >= gE_GMF_HM_danger[0] or \
            gE_GMF_HM[1] >= gE_GMF_HM_danger[1] or \
            gE_GMF_HM[2] >= gE_GMF_HM_danger[2]):
            
            GMF_SEV = True
            
# Gear tooth impact frequency and harmonics, CF=8, Severity = 10
# OR 
#     gE_GTIF_1 >= Alert
#     gE_GTIF_2 >= Alert
#     gE_GTIF_3 >= Alert              
#     gE_GTIF_4 >= Alert
#     gE_GTIF_5 >= Alert                       
                
        if (gE_GTIF_HM[0] >= gE_GTIF_HM_alert[0] or \
            gE_GTIF_HM[1] >= gE_GTIF_HM_alert[1] or \
            gE_GTIF_HM[2] >= gE_GTIF_HM_alert[2] or \
            gE_GTIF_HM[3] >= gE_GTIF_HM_alert[3] or \
            gE_GTIF_HM[4] >= gE_GTIF_HM_alert[4]):
            
            GTIF = True
            
         
        return   GMF_MOD, GMF_SEV, GTIF

###############################################################################
###############################################################################  
    # @staticmethod       
    def Gear_Fault_Diag(self,Acc_GMF_HM, Vel_GMF_HM, gE_GMF_HM, gE_GTIF_HM, Acc_GMF_alert_gPk = 0.299, Acc_GMF_danger_gPk = 0.344, \
                                                 Vel_GMF_alert_mmsRMS = 4.0, Vel_GMF_danger_mmsRMS = 8.0, gE_GMF_alert_gE = 3.0, \
                                                     gE_GMF_danger_gE = 6.0, gE_GTIF_alert_gE = 3.0, gE_GTIF_danger_gE = 6.0):
        
        GMF_MOD = False
        GMF_SEV = False
        GTIF = False

        
        Acc_GMF_HM_alert =[]
        Acc_GMF_HM_danger = []
        Vel_GMF_HM_alert =[]
        Vel_GMF_HM_danger = []
        gE_GMF_HM_alert = []
        gE_GMF_HM_danger = []        
        gE_GTIF_HM_alert = []
        gE_GTIF_HM_danger = []
        
        GearToothWearDetection = False
        GearToothWearSev = 0
        
        GearToothCrackBrokenDetection = False
        GearToothCrackBrokenSev = 0
        
        [Acc_GMF_HM_alert, Acc_GMF_HM_danger] = self.Acc_GMF_HM_Threshhold(Acc_GMF_alert_gPk, Acc_GMF_danger_gPk)
        [Vel_GMF_HM_alert, Vel_GMF_HM_danger] = self.Vel_GMF_HM_Threshhold(Vel_GMF_alert_mmsRMS, Vel_GMF_danger_mmsRMS)
        [gE_GMF_HM_alert, gE_GMF_HM_danger] = self.gE_GMF_HM_Threshhold(gE_GMF_alert_gE, gE_GMF_danger_gE)
        [gE_GTIF_HM_alert, gE_GTIF_HM_danger] = self.gE_GTIF_HM_Threshhold(gE_GTIF_alert_gE, gE_GTIF_danger_gE)
        
        [GMF_MOD, GMF_SEV, GTIF] = self.GMF_GTIF_Symptom_Rule(Vel_GMF_HM, Acc_GMF_HM, gE_GMF_HM, gE_GTIF_HM, Vel_GMF_HM_alert, \
                                                    Vel_GMF_HM_danger, Acc_GMF_HM_alert, Acc_GMF_HM_danger, \
                                                    gE_GMF_HM_alert, gE_GMF_HM_danger, gE_GTIF_HM_alert, gE_GTIF_HM_danger)

# Gear wear, misalignment, or eccentricity, sever, CF = 8, severity = 10
# AND
#   GMF_SEV = True
            
        if GMF_SEV == True:
            GearToothWearDetection = True
            GearToothWearSev = 10
            
# Gear wear, misalignment, or eccentricity, moderate, CF = 8, severity = 5
# AND
#   GMF_MOD = True
        
        else:
            if GMF_MOD == True:
                GearToothWearDetection = True
                GearToothWearSev = 5

# Gear tooth broken or cracked, CF = 8, severity = 5
# AND
#   GTIF = True
#   NOR
#       GMF_MOD = True
#       GMF_SEV = True
            
        if GTIF == True and (not(GMF_MOD == True or GMF_SEV == True)):
            GearToothCrackBrokenDetection = True
            GearToothCrackBrokenSev = 5
            
        return GearToothWearDetection, GearToothWearSev, GearToothCrackBrokenDetection, GearToothCrackBrokenSev
    


###############################################################################                
###############################################################################
    # @staticmethod 
# Note: the default alert and danger threshold values (in mm/s RMS for velocity and g peak for acceleration) are set according a DSS AC motor model, it can be adjusted if necessary       
    def RBPF_Symptom_Rule(self,Vel_RBPF_HM, Acc_RBPF_HM, Vel_RBPF_HM_alert, Vel_RBPF_HM_danger, Acc_RBPF_HM_alert, \
                          Acc_RBPF_HM_danger):
        RBPF_MOD = False
        RBPF_SEV = False
        
# Rotorbar pass frequency, sever, CF=8, severity = 10
# OR
#    Vel_RBPF_1 == Danger
#    Acc_RBPF_1 == Danger
            
        if (Vel_RBPF_HM[0] >= Vel_RBPF_HM_danger[0] or Acc_RBPF_HM[0] >= Acc_RBPF_HM_danger[0]):
            RBPF_SEV = True

# Rotorbar pass frequency, moderate, CF=8, severity = 5
# OR
#    Vel_RBPF_1 == Alert
#    Acc_RBPF_1 == Alert
        
        else:
            if (Vel_RBPF_HM[0] >= Vel_RBPF_HM_alert[0] or Acc_RBPF_HM[0] >= Acc_RBPF_HM_alert[0]) and  Vel_RBPF_HM[0] < Vel_RBPF_HM_danger[0] and Acc_RBPF_HM[0] < Acc_RBPF_HM_danger[0]:
                RBPF_MOD = True
      
        return    RBPF_MOD, RBPF_SEV
    


###############################################################################
###############################################################################
    # @staticmethod 
# Note: the default alert and danger threshold values (in mm/s RMS for velocity and g peak for acceleration) are set according a DSS AC motor model, it can be adjusted if necessary 
    def Cracked_Rotorbar_Diag(self,Vel_RBPF_HM, Acc_RBPF_HM, Vel_RBPF_alert_mmsRMS = 0.064*25.4, \
                              Vel_RBPF_danger_mmsRMS = 0.127*25.4, Acc_RBPF_alert_gPk = 2.5, Acc_RBPF_danger_gPk = 3.75):
        RBPF_MOD = False
        RBPF_SEV = False

        Vel_RBPF_HM_alert = []
        Vel_RBPF_HM_danger = []
        Acc_RBPF_HM_alert = []
        Acc_RBPF_HM_danger = []
        
        CrackedRotorbarDetection = False
        CrackedRotorbarSev = 0

        [Acc_RBPF_HM_alert, Acc_RBPF_HM_danger] = self.Acc_RBPF_HM_Threshhold(Acc_RBPF_alert_gPk, Acc_RBPF_danger_gPk)           
        [Vel_RBPF_HM_alert,Vel_RBPF_HM_danger] = self.Vel_RBPF_HM_Threshhold(Vel_RBPF_alert_mmsRMS, Vel_RBPF_danger_mmsRMS)

                
        [RBPF_MOD, RBPF_SEV] = self.RBPF_Symptom_Rule(Vel_RBPF_HM, Acc_RBPF_HM, Vel_RBPF_HM_alert, Vel_RBPF_HM_danger, Acc_RBPF_HM_alert, Acc_RBPF_HM_danger)

# Cracked rotorbar, sever  CF = 8, severity = 10
# AND
#    RBPF_SEV = True            
        if RBPF_SEV == True:
            CrackedRotorbarDetection = True
            CrackedRotorbarSev = 8 

# Cracked rotorbar, moderate  CF = 8, severity = 5
# AND
#    RBPF_MOD = True
#    NOR
#       RBPF_SEV = True        

        if RBPF_MOD == True and (not(RBPF_SEV == True)):
                CrackedRotorbarDetection = True
                CrackedRotorbarSev = 5

        return CrackedRotorbarDetection, CrackedRotorbarSev
    
###############################################################################                
###############################################################################
    # @staticmethod 
# Note: the default velocity alert and danger threshold (in mm/s RMS for velocity) is set according a DSS AC motor model, it can be adjusted if necessary       
    def LF2X_Symptom_Rule(self,Vel_2XLF_HM, gE_2XLF_HM, Vel_2XLF_HM_alert, Vel_2XLF_HM_danger, gE_2XLF_HM_alert, gE_2XLF_HM_danger):
        LF2X_MOD = False
        LF2X_SEV = False        

# 2X line frequency, moderate, CF=8, severity = 5
# AND
#    Vel_2XLF_1 = alert
        
        if Vel_2XLF_HM[0] >= Vel_2XLF_HM_alert[0] and Vel_2XLF_HM[0] < Vel_2XLF_HM_danger[0]:
            LF2X_MOD = True

# 2X line frequency, moderate, several defect peaks, no confirming velocity, CF=8, severity = 5
# AND
#   OR
#     AND
#       gE_2XLF_1 = alert
#       gE_2XLF_2 = alert
#     AND
#       gE_2XLF_1 = alert
#       gE_2XLF_3 = alert      
#     AND
#       gE_2XLF_1 = alert
#       gE_2XLF_4 = alert
#     AND
#       gE_2XLF_1 = alert
#       gE_2XLF_5 = alert
#   NOR
#     Vel_2XLF_1 >= alert

        else:
            if (((gE_2XLF_HM[0] >= gE_2XLF_HM_alert[0] and gE_2XLF_HM[0] < gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[1] >= gE_2XLF_HM_alert[1] and gE_2XLF_HM[1] < gE_2XLF_HM_danger[1])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_alert[0] and gE_2XLF_HM[0] < gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[2] >= gE_2XLF_HM_alert[2] and gE_2XLF_HM[2] < gE_2XLF_HM_danger[2])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_alert[0] and gE_2XLF_HM[0] < gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[3] >= gE_2XLF_HM_alert[3] and gE_2XLF_HM[3] < gE_2XLF_HM_danger[3])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_alert[0] and gE_2XLF_HM[0] < gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[4] >= gE_2XLF_HM_alert[4] and gE_2XLF_HM[4] < gE_2XLF_HM_danger[4]))) and \
                (not(Vel_2XLF_HM[0] >= Vel_2XLF_HM_alert[0])):
                LF2X_MOD = True

# 2X line frequency, severe, CF=8, severity = 10
# AND
#    Vel_2XLF_1  = Danger
            
        if Vel_2XLF_HM[0] >= Vel_2XLF_HM_danger[0]:
            LF2X_SEV = True

# 2X line frequency, severe, several defect peaks, no confirming velocity, CF=8, severity = 10
# AND
#   OR
#     AND
#       gE_2XLF_1 >= danger
#       gE_2XLF_2 >= alert
#     AND
#       gE_2XLF_1 >= danger
#       gE_2XLF_3 >= alert      
#     AND
#       gE_2XLF_1 >= danger
#       gE_2XLF_4 >= alert
#     AND
#       gE_2XLF_1 >= danger
#       gE_2XLF_5 >= alert
#   NOR
#     Vel_2XLF_1 >= alert
            
        else:
            if (((gE_2XLF_HM[0] >= gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[1] >= gE_2XLF_HM_alert[1])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[2] >= gE_2XLF_HM_alert[2])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[3] >= gE_2XLF_HM_alert[3])) or \
                ((gE_2XLF_HM[0] >= gE_2XLF_HM_danger[0]) and (gE_2XLF_HM[4] >= gE_2XLF_HM_alert[4]))) and \
                (not(Vel_2XLF_HM[0] < Vel_2XLF_HM_alert[0])):
                    LF2X_SEV = True            

# 2X line frequency, severe, confirmed by velocity spectrum, CF=8, severity = 10
# AND
#   OR
#     gE_2XLF_1 >= Alert
#     gE_2XLF_2 >= Alert
#   OR
#     Vel_2XLF_1 >= Alert
              
              
            else:
                 if (gE_2XLF_HM[0] >= gE_2XLF_HM_alert[0] or gE_2XLF_HM[1] >= gE_2XLF_HM_alert[1]) and \
                     Vel_2XLF_HM[0] >= Vel_2XLF_HM_alert[0]:
                     LF2X_SEV = True        
        
        return LF2X_MOD, LF2X_SEV
    
###############################################################################
###############################################################################
    # @staticmethod 
    def Variable_AirGap_Diag(self,Vel_2XLF_HM, gE_2XLF_HM, Vel_2XLF_alert = 0.064*25.4, Vel_2XLF_danger = 0.127*25.4, gE_2XLF_alert = 3, gE_2XLF_danger = 10):
        LF2X_MOD = False
        LF2X_SEV = False

        gE_2XLF_HM_alert = []
        gE_2XLF_HM_danger = []
        Vel_2XLF_HM_alert = []
        Vel_2XLF_HM_danger = []        
        
        VariableAirgapDetection = False
        VariableAirgapSev = 0
           
        [gE_2XLF_HM_alert,gE_2XLF_HM_danger] = self.gE_2XLF_HM_Threshhold(gE_2XLF_alert, gE_2XLF_danger)
        [Vel_2XLF_HM_alert,Vel_2XLF_HM_danger] = self.Vel_2XLF_HM_Threshhold(Vel_2XLF_alert, Vel_2XLF_danger)
        [LF2X_MOD, LF2X_SEV] = self.LF2X_Symptom_Rule(Vel_2XLF_HM, gE_2XLF_HM, Vel_2XLF_HM_alert, Vel_2XLF_HM_danger,gE_2XLF_HM_alert, gE_2XLF_HM_danger)

# Variable air gap, sever  CF = 8, severity = 10
# AND
#    LF2X_SEV = True           
        if LF2X_SEV == True:
            VariableAirgapDetection = True
            VariableAirgapSev = 8  

# Variable air gap, moderate  CF = 8, severity = 5
# AND
#    LF2X_MOD = True
#    NOR
#       LF2X_SEV = True        
        else:
            if LF2X_MOD == True and LF2X_SEV == False:
                VariableAirgapDetection = True
                VariableAirgapSev = 5
        
        return VariableAirgapDetection, VariableAirgapSev

###############################################################################
###############################################################################   
       
    def VPF_Symptom_Rule(self,Vel_VPF_HM, Vel_VPF_HM_alert, Vel_VPF_HM_danger):
        VAN_VPFH_MOD = False
        VAN_VPFH_SEV = False

# Vane pass frequency and harmonics (VPFH), sever, CF=8, Severity = 10
# OR
#    Vel_VPF_1 = Danger 
#    Vel_VPF_2 = Danger
                
        if Vel_VPF_HM[0] >= Vel_VPF_HM_danger[0] or Vel_VPF_HM[1] >= Vel_VPF_HM_danger[1] :
            VAN_VPFH_SEV = True
            

# Vane pass frequency and harmonics (VPFH), moderate, CF=8, Severity = 5
# AND
#    OR
#      Vel_VPF_1 = Alert 
#      Vel_VPF_2 = Alert
#    NOR
#      Vel_VPF_1 = Danger 
#      Vel_VPF_2 = Danger        
            
        
        if (Vel_VPF_HM[0] >= Vel_VPF_HM_alert[0] and Vel_VPF_HM[0] < Vel_VPF_HM_danger[0] or\
            Vel_VPF_HM[1] >= Vel_VPF_HM_alert[1] and Vel_VPF_HM[1] < Vel_VPF_HM_danger[1]) and\
            (not(Vel_VPF_HM[0] >= Vel_VPF_HM_danger[0] or Vel_VPF_HM[1] >= Vel_VPF_HM_danger[1])):
            VAN_VPFH_MOD = True

        return VAN_VPFH_MOD, VAN_VPFH_SEV
    
###############################################################################
###############################################################################
    # @staticmethod 
    def Vane_WearDamage_Diag(self,Vel_VPF_HM, Vel_VPF_alert_mmsRMS = 0.159*25.4, \
                                             Vel_VPF_danger_mmsRMS = 0.318*25.4):
        VAN_VPFH_MOD = False
        VAN_VPFH_SEV = False

        Vel_VPF_HM_alert = []
        Vel_VPF_HM_danger = []
        
        VaneWearDamageDetection = False
        VaneWearDamageSev = 0
           
        [Vel_VPF_HM_alert,Vel_VPF_HM_danger] = self.Vel_VPF_HM_Threshhold(Vel_VPF_alert_mmsRMS, Vel_VPF_danger_mmsRMS)
        [VAN_VPFH_MOD, VAN_VPFH_SEV] = self.VPF_Symptom_Rule(Vel_VPF_HM, Vel_VPF_HM_alert, Vel_VPF_HM_danger)

# Vane wear or damage, severe CF = 8, severity = 10
# OR
#    VAN_VPFH_SEV= True
       
        if VAN_VPFH_SEV == True:
            VaneWearDamageDetection = True
            VaneWearDamageSev = 10

# Vane wear or damage, moderate  CF = 6, severity = 5
# OR
#    VAN_VPFH_MOD= True
       
        else:
            if VAN_VPFH_MOD == True:
                VaneWearDamageDetection = True
                VaneWearDamageSev = 5
        
        return VaneWearDamageDetection, VaneWearDamageSev

###############################################################################
###############################################################################
    # 巴特沃斯带通滤波
    # xTimeWave是时域波形数组， fs是采样频率，单位为Hz
    # lowcut 为低截止频率，单位为Hz
    # highcut 为高截止频率，单位为Hz
    # Order为滤波器阶数，通常为偶数，2，4，6，8一般不超过8
    def BandPassButter(self,xTimeWave, lowcut, highcut, fs, order):
        low = lowcut * 2 / fs
        high = highcut * 2 / fs
        sos = signal.butter(order, [low, high], 'bandpass', output='sos')
        filtered = signal.sosfilt(sos, xTimeWave)
        return filtered