温度诊断与自动化诊断合并为故障诊断

.gitignore

@@ -0,0 +1,2 @@
+__pycache__/
+*.pyc

DiagnosisLib.py

@@ -0,0 +1,2006 @@
+
+"""
+@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
+#         if fRangeMode == 'Per':
+#             fRange = fRange*fTarget
+#         if fRangeMode == 'Num':
+#             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 
+    
+
+    def FixedFreqFeature(self, xSpec, fMax, fTarget, fRange, fRangeMode, Detection, numHarmonics):
+        if fRangeMode == 'Per':
+            fRange = fRange * fTarget
+        if fRangeMode == 'Num':
+            fRange = 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
+
+        # 检查索引是否有效
+        if index_SearchRangeMin >= len(xSpec) or index_SearchRangeMax >= len(xSpec):
+            return [0] * (numHarmonics + 1)  # 返回全零数组
+
+        index_SearchRange = index_SearchRangeMax - index_SearchRangeMin
+
+        # 检查搜索范围是否有效
+        if index_SearchRange <= 0:
+            return [0] * (numHarmonics + 1)  # 返回全零数组
+
+        xArray = xSpec[index_SearchRangeMin:(index_SearchRangeMax + 1)]
+
+        # 检查 xArray 是否为空
+        if len(xArray) == 0:
+            return [0] * (numHarmonics + 1)  # 返回全零数组
+
+        # 确保 xArray 不为空
+        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))]
+
+            # 再次检查 xArray 是否为空
+            if len(xArray) == 0:
+                return [0] * (numHarmonics + 1)  # 返回全零数组
+
+        ActualTargetFreq = self.FundamentalFreqActualValue(xArray, SpecResolution, index_SearchRangeMin)
+        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]
+        if Detection == 'Peak':
+            pass
+        if Detection == 'RMS':
+            featureValue = [x/math.sqrt(2) for x in featureValue]
+        if Detection == 'PP':
+            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 , Vel_danger, defectType ):
+        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

Temp_Diag.py

@@ -0,0 +1,184 @@
+import numpy as np
+import pandas as pd
+from sklearn.neighbors import BallTree
+from sqlalchemy import create_engine, text
+import math, joblib, os
+
+class MSET_Temp:
+    """
+    MSET + SPRT 温度分析类：
+    - 离线训练：genDLMatrix → save_model
+    - 在线推理：load_model → predict_SPRT
+    """
+
+    def __init__(self,
+                 windCode: str,
+                 windTurbineNumberList: list[str],
+                 startTime: str,
+                 endTime: str):
+        self.windCode = windCode.strip()
+        self.windTurbineNumberList = windTurbineNumberList or []
+        self.startTime = startTime
+        self.endTime   = endTime
+
+        # 离线训练／加载后赋值
+        self.matrixD = None
+        self.healthyResidual = None
+        self.normalDataBallTree = None
+
+        # SPRT 参数（离线训练时设置）
+        self.feature_weight: np.ndarray | None = None
+        self.alpha: float = 0.1
+        self.beta:  float = 0.1
+
+    def _get_data_by_filter(self) -> pd.DataFrame:
+        """
+        在线推理专用：根据 self.windTurbineNumberList & 时间拉数据；
+        如果列表为空，则拉全场数据。
+        """
+        table = f"{self.windCode}_minute"
+        engine = create_engine(
+            #"mysql+pymysql://root:admin123456@106.120.102.238:10336/energy_data_prod"
+             "mysql+pymysql://root:admin123456@192.168.50.235:30306/energy_data_prod"
+        )
+        if self.windTurbineNumberList:
+            turbines = ",".join(f"'{t}'" for t in self.windTurbineNumberList)
+            cond = f"wind_turbine_number IN ({turbines}) AND time_stamp BETWEEN :start AND :end"
+        else:
+            cond = "time_stamp BETWEEN :start AND :end"
+
+        sql = text(f"""
+            SELECT *
+            FROM {table}
+            WHERE {cond}
+            ORDER BY time_stamp ASC
+        """)
+        return pd.read_sql(sql, engine, params={"start": self.startTime, "end": self.endTime})
+
+    def calcSimilarity(self, x: np.ndarray, y: np.ndarray, m: str = 'euc') -> float:
+        if len(x) != len(y):
+            return 0.0
+        if m == 'cbd':
+            return float(np.mean([1.0/(1.0+abs(p-q)) for p,q in zip(x,y)]))
+        diffsq = np.sum((x-y)**2)
+        return float(1.0/(1.0+math.sqrt(diffsq)))
+
+    def genDLMatrix(self, trainDataset: np.ndarray,
+                    dataSize4D=100, dataSize4L=50) -> int:
+        """
+        离线训练：构造 matrixD/matrixL/healthyResidual/BallTree
+        """
+        m, n = trainDataset.shape
+        if m < dataSize4D + dataSize4L:
+            return -1
+
+        # Step1：每维最小/最大入 D
+        D_idx, D = [], []
+        for i in range(n):
+            col = trainDataset[:, i]
+            for idx in (np.argmin(col), np.argmax(col)):
+                D.append(trainDataset[idx].tolist())
+                D_idx.append(idx)
+        # Step2：挑样本至 dataSize4D
+        while len(D_idx) < dataSize4D:
+            free = list(set(range(m)) - set(D_idx))
+            scores = [(np.mean([1-self.calcSimilarity(trainDataset[i], d) for d in D]), i)
+                      for i in free]
+            _, pick = max(scores)
+            D.append(trainDataset[pick].tolist())
+            D_idx.append(pick)
+        self.matrixD = np.array(D)
+
+        # BallTree + healthyResidual
+        self.normalDataBallTree = BallTree(
+            self.matrixD,
+            leaf_size=4,
+            metric=lambda a,b: 1.0 - self.calcSimilarity(a, b)
+        )
+        # healthyResidual
+        ests = []
+        for x in trainDataset:
+            dist, idxs = self.normalDataBallTree.query([x], k=20, return_distance=True)
+            w = 1.0/(dist[0]+1e-1)
+            w /= w.sum()
+            ests.append(np.sum([wi*self.matrixD[j] for wi,j in zip(w,idxs[0])], axis=0))
+        self.healthyResidual = np.array(ests) - trainDataset
+        return 0
+
+    def calcSPRT(self,
+                 newsStates: np.ndarray,
+                 feature_weight: np.ndarray,
+                 alpha: float = 0.1,
+                 beta: float = 0.1,
+                 decisionGroup: int = 5) -> list[float]:
+        """
+        Wald-SPRT 得分
+        """
+        # 新状态残差
+        ests = []
+        for x in newsStates:
+            dist, idxs = self.normalDataBallTree.query([x], k=20, return_distance=True)
+            w = 1.0/(dist[0]+1e-1); w/=w.sum()
+            ests.append(np.sum([wi*self.matrixD[j] for wi,j in zip(w,idxs[0])], axis=0))
+        resN = np.array(ests) - newsStates
+
+        # 加权
+        wN = [np.dot(r, feature_weight) for r in resN]
+        wH = [np.dot(r, feature_weight) for r in self.healthyResidual]
+        mu0, sigma0 = np.mean(wH), np.std(wH)
+        low = math.log(beta/(1-alpha)); high = math.log((1-beta)/alpha)
+
+        flags = []
+        for i in range(len(wN)-decisionGroup+1):
+            seg = wN[i:i+decisionGroup]; mu1=np.mean(seg)
+            si = (sum(seg)*(mu1-mu0)/sigma0**2
+                  - decisionGroup*((mu1**2-mu0**2)/(2*sigma0**2)))
+            si = max(min(si, high), low)
+            flags.append(si/high if si>0 else si/low)
+        return flags
+
+    def predict_SPRT(self,
+                     newsStates: np.ndarray,
+                     decisionGroup: int = 5) -> list[float]:
+        """
+        在线推理：用离线保存的 matrixD/healthyResidual/feature_weight/alpha/beta
+        """
+        return self.calcSPRT(
+            newsStates,
+            self.feature_weight,
+            alpha=self.alpha,
+            beta=self.beta,
+            decisionGroup=decisionGroup
+        )
+
+    def save_model(self, path: str):
+        """
+        Save matrixD, healthyResidual, feature_weight, alpha, beta
+        """
+        os.makedirs(os.path.dirname(path), exist_ok=True)
+        joblib.dump({
+            'matrixD': self.matrixD,
+            'healthyResidual': self.healthyResidual,
+            'feature_weight': self.feature_weight,
+            'alpha': self.alpha,
+            'beta': self.beta,
+        }, path)
+
+    @classmethod
+    def load_model(cls, path: str) -> 'MSET_Temp':
+        """
+        Load + rebuild BallTree
+        """
+        data = joblib.load(path)
+        inst = cls('', [], '', '')
+        inst.matrixD = data['matrixD']
+        inst.healthyResidual = data['healthyResidual']
+        inst.feature_weight = data['feature_weight']
+        inst.alpha = data['alpha']
+        inst.beta  = data['beta']
+        inst.normalDataBallTree = BallTree(
+            inst.matrixD,
+            leaf_size=4,
+            metric=lambda a,b: 1.0 - inst.calcSimilarity(a, b)
+        )
+        return inst

api_diag.py

@@ -0,0 +1,285 @@
+# main.py
+
+import os, glob
+import pandas as pd
+import json
+import shutil
+import pandas as pd
+import uvicorn
+from fastapi import FastAPI, UploadFile, HTTPException
+from fastapi.responses import JSONResponse
+from pydantic import BaseModel, Field, model_validator
+from typing import List, Optional, Union
+from autodiag_class import Auto_diag
+from Temp_Diag import MSET_Temp
+
+app = FastAPI(root_path="/api/diag",title=" Diagnosis API")
+
+# 全局：{ windCode: { turbine: { channel: model, … }, … }, … }
+MODEL_STORE: dict[str, dict[str, dict[str, MSET_Temp]]] = {}
+
+# 英文→中文映射
+cn_map = {
+    'main_bearing_temperature': '主轴承温度',
+    'gearbox_oil_temperature': '齿轮箱油温',
+    'generatordrive_end_bearing_temperature': '发电机驱动端轴承温度',
+    'generatornon_drive_end_bearing_temperature': '发电机非驱动端轴承温度'
+}
+
+class TemperatureInput(BaseModel):
+    windCode: str
+    windTurbineNumberList: List[str]
+    startTime: str
+    endTime: str
+
+    @model_validator(mode='before')
+    def ensure_list(cls, v):
+        raw = v.get('windTurbineNumberList')
+        if isinstance(raw, str):
+            v['windTurbineNumberList'] = [raw]
+        return v
+
+class TemperatureThresholdInput(TemperatureInput):
+    pageNo:   int
+    pageSize: int
+
+@app.on_event("startup")
+def load_all_models():
+    for f in glob.glob("models/*/*/*.pkl"):
+        _, wc, turbine, fname = f.split(os.sep)
+        ch = os.path.splitext(fname)[0]
+        MODEL_STORE.setdefault(wc, {}).setdefault(turbine, {})[ch] = MSET_Temp.load_model(f)
+    print("模型加载完成:", {k: list(v.keys()) for k,v in MODEL_STORE.items()})
+
+@app.post("/temperature/threshold")
+async def route_threshold(inp: TemperatureThresholdInput):
+    """
+    输入：
+    {
+      "windCode": "WOF091200030",
+      "windTurbineNumberList": ["WOG01355"],
+      "startTime": "2024-06-01 00:00",
+      "endTime":   "2024-06-05 01:00",
+      "pageNo":    1,
+      "pageSize":  10
+    }
+    输出：
+    {
+      "data": {
+        "type": "temperature_threshold",
+        "records": [
+          {
+            "wind_turbine_number": "WOG01355",
+            "time_stamp": "2025-06-01 00:05:00",
+            "temp_channel": "主轴承温度",
+            "SPRT_score": 0.12,
+            "status": "正常"
+          },
+          ...
+        ],
+        "totalSize": 42
+      },
+      "code": 200,
+      "message": "success"
+    }
+    """
+    # 1) 校验模型是否存在
+    if inp.windCode not in MODEL_STORE:
+        raise HTTPException(404, f"无模型：{inp.windCode}")
+    # 2) 为每台待分析风机，拉数据并推理
+    records = []
+    for turbine in inp.windTurbineNumberList:
+        if turbine not in MODEL_STORE[inp.windCode]:
+            continue
+        analyzer = MSET_Temp(inp.windCode, [turbine], inp.startTime, inp.endTime)
+        df = analyzer._get_data_by_filter()
+        if df.empty:
+            continue
+        df['time_stamp'] = pd.to_datetime(df['time_stamp'])
+
+        for eng, cn in cn_map.items():
+            if eng not in df.columns:
+                continue
+            sub = df[['time_stamp', eng]].dropna()
+            arr = sub[eng].values.reshape(-1,1)
+            ts  = sub['time_stamp'].dt.strftime("%Y-%m-%d %H:%M:%S").tolist()
+
+            model = MODEL_STORE[inp.windCode][turbine].get(eng)
+            if not model:
+                continue
+            flags = model.predict_SPRT(arr, decisionGroup=1)
+
+            for i, sc in enumerate(flags):
+                records.append({
+                    "wind_turbine_number": turbine,
+                    "time_stamp": ts[i],
+                    "temp_channel": cn,
+                    "SPRT_score": sc,
+                    "status": "危险" if sc>0.99 else "正常"
+                })
+
+    # 分页返回
+    total = len(records)
+    start = (inp.pageNo-1)*inp.pageSize
+    end   = start+inp.pageSize
+    return {
+        "data": {
+            "type": "temperature_threshold",
+            "records": records[start:end],
+            "totalSize": total
+        },
+        "code": 200,
+        "message": "success"
+    }
+
+@app.post("/SPRT/trend")
+async def route_trend(inp: TemperatureInput):
+    """
+    输入：
+    {
+      "windCode": "WOF091200030",
+      "windTurbineNumberList": ["WOG01355"],
+      "startTime": "2024-06-01 00:00",
+      "endTime":   "2024-06-05 01:00"
+    }
+    输出：
+    {
+      "data": {
+        "type": "SPRT_trend",
+        "main_bearing": {"timestamps": [...], "values": [...]},
+        "gearbox_oil": {"timestamps": [...], "values": [...]},
+        "generator_drive_end": {"timestamps": [...], "values": [...]},
+        "generator_nondrive_end": {"timestamps": [...], "values": [...]}
+      },
+      "code": 200,
+      "message": "success"
+    }
+
+    """
+    if inp.windCode not in MODEL_STORE:
+        raise HTTPException(404, f"无模型：{inp.windCode}")
+
+    turbines_out = []
+    for turbine in inp.windTurbineNumberList:
+        if turbine not in MODEL_STORE[inp.windCode]:
+            continue
+        analyzer = MSET_Temp(inp.windCode, [turbine], inp.startTime, inp.endTime)
+        df = analyzer._get_data_by_filter()
+        if df.empty:
+            continue
+        df['time_stamp'] = pd.to_datetime(df['time_stamp'])
+
+        ch_data = {}
+        for eng, key in {
+            'main_bearing_temperature':'main_bearing',
+            'gearbox_oil_temperature':'gearbox_oil',
+            'generatordrive_end_bearing_temperature':'generator_drive_end',
+            'generatornon_drive_end_bearing_temperature':'generator_nondrive_end'
+        }.items():
+            if eng not in df.columns:
+                ch_data[key] = {"timestamps": [], "values": []}
+                continue
+            sub = df[['time_stamp', eng]].dropna()
+            ts  = sub['time_stamp'].dt.strftime("%Y-%m-%d %H:%M:%S").tolist()
+            arr = sub[eng].values.reshape(-1,1)
+
+            model = MODEL_STORE[inp.windCode][turbine].get(eng)
+            vals = model.predict_SPRT(arr, decisionGroup=1) if model else []
+            ch_data[key] = {"timestamps": ts, "values": vals}
+
+        # turbines_out.append({
+        #   "wind_turbine_number": turbine,
+        #     **ch_data
+        # })
+
+    return {
+        "data": {
+            "type": "SPRT_trend",
+            **ch_data
+        },
+        "code": 200,
+        "message": "success"
+    }
+
+if __name__ == "__main__":
+    import uvicorn
+    uvicorn.run("main:app", host="0.0.0.0", port=8000, reload=True)
+
+
+
+# 请求模型定义
+class AutoDiagInput(BaseModel):
+    ids: List[int]      # 数据索引id
+    windCode: str       # 风场编号
+    engine_code: str    # 风机编号
+    autodiagType: str   # 诊断类型
+
+    @model_validator(mode='before')
+    def convert_ids(cls, values):
+        """将单个id转换为列表形式"""
+        if isinstance(values.get('ids'), int):
+            values['ids'] = [values['ids']]
+        return values
+
+class DiagnosisResult(BaseModel):
+    status_codes: List[int]  # 每个id对应的状态码列表
+    max_status: int          # 所有状态码中的最大值
+    count_0: int             # 状态码0的个数
+    count_1: int             # 状态码1的个数
+    count_2: int             # 状态码2的个数
+
+@app.post("/autodiag/{autodiagType}")
+async def perform_diagnosis(autodiagType: str, input_data: AutoDiagInput):
+    """
+    执行自动诊断分析
+    参数:
+        autodiagType: 诊断类型
+        input_data: 包含ids, windCode, engine_code的输入数据
+    返回:
+        诊断结果，包含状态码列表和统计信息
+    """
+    autodiag_map = {
+        "Unbalance": "Unbalance_diag",    # 不平衡诊断
+        "Misalignment": "Misalignment_diag",  # 不对中诊断
+        "Looseness": "Looseness_diag",    # 松动诊断
+        "Bearing": "Bearing_diag",        # 轴承诊断       
+        "Gear": "Gear_diag"              # 齿轮诊断
+    }
+
+    if autodiagType not in autodiag_map:
+        raise HTTPException(status_code=400, detail="非可用的诊断类型")
+    
+    try:
+        # 初始化诊断类
+        autodiag = Auto_diag(input_data.ids, input_data.windCode, input_data.engine_code)
+        
+        # 获取诊断方法
+        func = getattr(autodiag, autodiag_map[autodiagType])
+        
+        # 执行诊断
+        if callable(func):
+            result = func()  # 直接返回格式化后的结果
+            return JSONResponse(content=result)
+            
+    except ValueError as e:
+        # 专门捕获齿轮诊断的错误
+        if "Can not perform gearbox diagnosis" in str(e):
+            return JSONResponse(
+                status_code=200,
+                content={
+                    "code": 400,
+                    "message": str(e)
+                }
+            )
+        elif "当前采集频率不适合进行诊断分析" in str(e):
+            return JSONResponse(
+                status_code=200,
+                content={
+                    "code": 405,
+                    "message": str(e)
+                }
+            )            
+        # 其他ValueError
+        raise HTTPException(status_code=400, detail=str(e))
+    except Exception as e:
+        raise HTTPException(status_code=500, detail=str(e))

autodiag_class.py

@@ -0,0 +1,1130 @@
+import numpy as np
+import pandas as pd
+from DiagnosisLib import DiagnosisLib
+from typing import List, Dict, Optional
+import json
+from sqlalchemy import create_engine, text
+import math
+
+class Auto_diag:
+    def __init__(self, ids: List[int], wind_code: str, engine_code: str):
+        """
+        初始化诊断类
+        参数:
+            ids: 数据ID列表
+            wind_code: 风场编号
+            engine_code: 风机编号
+        """
+        self.ids = ids
+        self.wind_code = wind_code
+        self.engine_code = engine_code
+        self.myDiagnosisLib = DiagnosisLib()        
+        self.diagnosis_results = []  # 存储所有诊断结果
+        self.bearing_frequencies = {}  # 存储特征频率
+        self.teeth_count={} #存储齿数
+        # 从数据库加载振动数据
+        self.load_vibration_data()
+        # 计算转速传动比
+        self.calculate_speed_ratio()
+        # 提取轴承参数、齿轮参数 计算特征频率
+        self.Characteristic_Frequency()
+        # 计算频谱
+        self.calculate_spectrums()
+    
+
+    def load_vibration_data(self):
+        """从数据库中获取振动加速度数据"""
+        try:
+            engine = create_engine('mysql+pymysql://root:admin123456@192.168.50.235:30306/energy_data_prod')
+            #engine = create_engine('mysql+pymysql://root:admin123456@106.120.102.238:10336/energy_data_prod')         
+            with engine.connect() as conn:
+                table_name = self.wind_code + '_wave'
+                # 将所有id拼接到SQL语句中，一次性查询所有数据
+                ids_str = ','.join(map(str, self.ids))
+                sql = f"SELECT * FROM {table_name} WHERE id IN ({ids_str})"
+                df = pd.read_sql(sql, conn)
+                self.df_data = df
+                
+                # 提取公共参数
+                self.fspd = self.df_data['rotational_speed'].astype(float).values  # 转速
+                self.fs_acc = self.df_data['sampling_frequency'].astype(float).values  # 采样频率
+                self.mesure_point_names = self.df_data['mesure_point_name'].values
+                
+                # 处理加速度数据
+                self.Acc = []
+                self.t_acc = []
+                for i, row in self.df_data.iterrows():
+                    measure_data_str = row['mesure_data']
+                    try:
+                        measure_data_list = json.loads(measure_data_str)
+                        acc_data = np.array(measure_data_list, dtype=np.float32)
+                    except json.JSONDecodeError:
+                        stripped_str = measure_data_str.strip('[]')
+                        measure_data_list = [float(x.strip()) for x in stripped_str.split(',') if x.strip()]
+                        acc_data = np.array(measure_data_list, dtype=np.float32)
+                    
+                    self.Acc.append(acc_data)
+                    # 使用转换后的数值类型计算时间轴
+                    self.t_acc.append(np.arange(0, len(acc_data)/self.fs_acc[i], 1/self.fs_acc[i]))
+
+        except Exception as e:
+            raise Exception(f"Failed to load vibration data: {str(e)}")
+        finally:
+            engine.dispose()
+
+
+    def calculate_spectrums(self):
+        """计算各种频谱"""
+        print("计算频谱")
+        self.fAcc = []
+        self.AccSpec = []
+        self.VelSpec = []
+        self.gE_10 = []
+        self.fAcc_gE = []
+        self.gESpec_10 = []
+        if  self.fs_acc[0] < 12800:
+            print("当前采集频率不适合进行诊断分析")
+            raise ValueError("当前采集频率不适合进行诊断分析")      
+             
+        for i, acc in enumerate(self.Acc):
+            # 获取当前 id 对应的采样频率
+            fs_acc = self.fs_acc[i]
+
+            fmax_acc =fs_acc/2.56
+            print(fs_acc)
+            print(acc)
+            lowcut = 500 #最低值是固定的 /51.2
+            highcut =fs_acc /2.56
+            # 计算振动加速度频谱
+            [f, spec] = self.myDiagnosisLib.Spectrum(acc, fs_acc)
+            self.fAcc.append(f)
+            self.AccSpec.append(spec)
+
+            # 计算振动速度频谱
+            fCut = 0.1
+            vel_spec = self.myDiagnosisLib.FreqIntegration(spec, fmax_acc, fCut)
+            self.VelSpec.append(vel_spec)
+
+            # 计算加速度包络频谱
+            order = 6
+            acc_band = self.myDiagnosisLib.BandPassButter(acc, lowcut, highcut, fs_acc, order)
+
+
+            fmax_gE_10 =1000
+            DownSampleRate =int(fs_acc /2.56/fmax_gE_10 ) #降采样
+            gE = self.myDiagnosisLib.EnvLopInTime(acc_band, DownSampleRate, Method="SKF")
+            self.gE_10.append(gE)
+
+            fs_gE_10 = 2.56 * fmax_gE_10
+            [f, spec] = self.myDiagnosisLib.Spectrum(gE, fs_gE_10)
+            self.fAcc_gE.append(f)
+            self.gESpec_10.append(spec)
+    
+    def Unbalance_diag(self):
+        """不平衡诊断"""
+        print("不平衡诊断")
+        status_codes = []
+        
+        for i in range(len(self.ids)):
+            ShaftSpd = self.fspd[i]*self.speed_ratio
+
+            if ShaftSpd ==0:
+                status_codes.append(-1)    
+                continue    
+            try:       
+                fs_acc = self.fs_acc[i]
+                fmax_vel = 1000
+                # fmax_vel =max(1000,fs_acc/2.56)  
+                print('ShaftSpd',ShaftSpd)
+                fTarget = ShaftSpd /60
+                print('fTarget',fTarget)
+                numHarmonics = 4 
+            
+                fRange = 0.05
+                fRangeMode = "Per"
+                Detection = "RMS"
+                print(self.VelSpec[i])
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                   # print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                     
+                print('Vel_RunningSpd_HM',Vel_RunningSpd_HM)           
+                Vel_alert_mmsRMS = 4  # 速度报警阈值
+                Vel_danger_mmsRMS = 7  # 速度危险阈值
+                defectType = 'Unbalance'
+                try:                
+                    detection, severity = self.myDiagnosisLib.RunningSpd_Fault_Diag(
+                        Vel_RunningSpd_HM, Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType
+                    )
+                except Exception as e:
+                    detection =False
+                    severity =0                   
+                print("severity",severity)  
+                if not detection:
+                    status_code = 0
+                elif severity < 6:
+                    status_code = 1
+                else:
+                    status_code = 2
+            except Exception as e:    
+                status_code =0    
+            status_codes.append(status_code)
+            print("status_code",status_code)       
+        # 返回结果和统计信息
+        return {
+            "type": "Unbalance",
+            "results": status_codes,
+            "statistics": {
+                "max_status": max(status_codes),
+                "count_0": status_codes.count(0),
+                "count_1": status_codes.count(1),
+                "count_2": status_codes.count(2)
+            }
+        }
+    
+    def Misalignment_diag(self):
+        """不对中诊断"""
+        status_codes = []
+        print("不对中诊断")        
+        for i in range(len(self.ids)):
+            ShaftSpd = self.fspd[i]*self.speed_ratio
+            if ShaftSpd ==0:
+                status_codes.append(-1)    
+                continue    
+            try:
+                fs_acc = self.fs_acc[i]
+                fmax_vel = 1000
+                # fmax_vel =max(1000,fs_acc/2.56)  
+                fTarget = ShaftSpd /60
+                numHarmonics = 4 
+                fRange = 0.05
+                fRangeMode = "Per"
+                Detection = "RMS"
+                print('ShaftSpd',ShaftSpd)            
+                print('fTarget',fTarget)
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                   # print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                        
+                print('Vel_RunningSpd_HM',Vel_RunningSpd_HM)               
+                Vel_alert_mmsRMS = 4  # 速度报警阈值
+                Vel_danger_mmsRMS = 7  # 速度危险阈值
+
+                try:
+                    defectType = 'Misalignment'
+                    detection, severity = self.myDiagnosisLib.RunningSpd_Fault_Diag(
+                        Vel_RunningSpd_HM, Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType
+                    )
+                except Exception as e:
+                    detection = False
+                    severity =0                  
+                print('severity',severity)                 
+                if not detection:
+                    status_code = 0
+                elif severity < 6:
+                    status_code = 1
+                else:
+                    status_code = 2
+            except Exception as e:   
+                status_code =0
+            status_codes.append(status_code)
+        
+        return {
+            "type": "Misalignment",
+            "results": status_codes,
+            "statistics": {
+                "max_status": max(status_codes),
+                "count_0": status_codes.count(0),
+                "count_1": status_codes.count(1),
+                "count_2": status_codes.count(2)
+            }
+        }
+    
+    def Looseness_diag(self):
+        """松动诊断"""
+        status_codes = []
+        print("松动诊断")         
+        for i in range(len(self.ids)):
+            ShaftSpd = self.fspd[i]*self.speed_ratio
+            if ShaftSpd ==0:
+                status_codes.append(-1)    
+                continue    
+            try:
+                fs_acc = self.fs_acc[i]
+                fmax_vel = 1000
+                # fmax_vel =max(1000,fs_acc/2.56)  
+                fTarget = ShaftSpd /60
+                numHarmonics = 4 
+                fRange = 0.05
+                fRangeMode = "Per"
+                Detection = "RMS"
+
+                print('ShaftSpd',ShaftSpd)
+                print('fTarget',fTarget)         
+                try:   
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                   # print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                       
+                print('Vel_RunningSpd_HM',Vel_RunningSpd_HM)            
+                Vel_alert_mmsRMS = 4  # 速度报警阈值
+                Vel_danger_mmsRMS = 7  # 速度危险阈值
+
+                try:                   
+                    defectType = 'Looseness'                
+                    detection, severity = self.myDiagnosisLib.RunningSpd_Fault_Diag(
+                        Vel_RunningSpd_HM, Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType
+                    )
+                except Exception as e:
+                    detection =False
+                    severity =0                    
+                print('severity',severity)             
+                if not detection:
+                    status_code = 0
+                elif severity < 6:
+                    status_code = 1
+                else:
+                    status_code = 2
+            except Exception as e:    
+                status_code =0            
+            status_codes.append(status_code)
+        
+        return {
+            "type": "Looseness",
+            "results": status_codes,
+            "statistics": {
+                "max_status": max(status_codes),
+                "count_0": status_codes.count(0),
+                "count_1": status_codes.count(1),
+                "count_2": status_codes.count(2)
+            }
+        }
+ 
+    
+    def Bearing_diag(self):
+        """轴承诊断"""
+        print("轴承诊断")
+        # Vel_RunningSpd_HM，生成转频的1X-5X的倍频值
+        # VelSpec，振动速度频谱
+        # fmax_vel，最大分析频率
+        # fTarget，目标频率为转频
+        Vel_alert_mmsRMS = 4  # 速度总值预警阈值 （mm/s 有效值）
+        Vel_danger_mmsRMS = 7  # 速度总值报警阈值 （mm/s 有效值）
+        gE_alert = 3  # 包络总值预警阈值（gE）
+        gE_danger = 10  # 包络总值报警阈值（gE）
+        results = {
+            "BPFO": {"status_codes": []},
+            "BPFI": {"status_codes": []},
+            "BSF": {"status_codes": []},
+            "FTF": {"status_codes": []}
+        }
+        # 获取第一个id的特征频率(所有id相同)
+        id = self.ids[0]
+        bearing_freq = self.bearing_frequencies.get(id, {})
+        BPFO = bearing_freq.get("BPFO", 0)
+        BPFI = bearing_freq.get("BPFI", 0)
+        BSF = bearing_freq.get("BSF", 0)
+        FTF = bearing_freq.get("FTF", 0)
+
+        for i in range(len(self.ids)):
+            ShaftSpd = self.fspd[i] * self.speed_ratio
+            if ShaftSpd ==0:
+                results["BPFO"]["status_codes"].append(-1)
+                results["BPFI"]["status_codes"].append(-1)
+                results["BSF"]["status_codes"].append(-1)
+                results["FTF"]["status_codes"].append(-1)
+                continue    
+            fs_acc = self.fs_acc[i]
+            # fmax_vel =1000
+            #fmax_vel如果直接取1000 还是会出现目标范围落不到搜寻范围的现象
+            fmax_vel = max(1000,fs_acc/2.56)
+            fRange = 0.05
+            fmax_gE_10 = 1000
+            try:
+                """外圈故障诊断"""
+                print("外圈")
+                fbr_BPFO = BPFO              
+                fTarget = fbr_BPFO #在计算轴承特征频率的函数中已经乘过转速   
+
+                if self.wind_code == 'WOF091200030':
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <1:
+                        fmax_gE_10 =10         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000     
+                else:
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <10:
+                        fmax_gE_10 =100         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000                               
+                numHarmonics = 4
+            
+                fRangeMode = "Per"
+                Detection = "RMS"    
+                print('BPFO',BPFO)    
+                print('ShaftSpd',ShaftSpd)     
+                print('speed_ratio',self.speed_ratio)           
+                print('fTarget',fTarget)            
+                print('VelSpec',self.VelSpec[i])
+                print('gESpec_10',self.gESpec_10[i])    
+
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(self.VelSpec[i], fmax_vel, fTarget, fRange, \
+                                                            fRangeMode, Detection, numHarmonics)    
+                except Exception as e:
+                   # print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                                  
+                print('Vel_RunningSpd_HM\n',Vel_RunningSpd_HM)       
+
+                # 速度频谱特征
+                try:
+                    Vel_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 Vel_Bearing_Defect_HM 时出错: {e}")
+                    Vel_Bearing_Defect_HM = [0] * (numHarmonics + 1)                        
+                print('Vel_Bearing_Defect_HM\n',Vel_Bearing_Defect_HM)
+
+                try:
+                    # 包络频谱特征
+                    Detection = "Peak" 
+                    gE_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 gE_Bearing_Defect_HM 时出错: {e}")
+                    gE_Bearing_Defect_HM = [0] * (numHarmonics + 1)                          
+                print('gE_Bearing_Defect_HM\n',gE_Bearing_Defect_HM)
+
+                try:
+                    # 诊断外圈故障
+                    BPFOdetection, BPFOseverity = self.myDiagnosisLib.Bearing_Fault_Diag(
+                        gE_Bearing_Defect_HM, Vel_Bearing_Defect_HM, Vel_RunningSpd_HM, gE_alert, gE_danger,
+                        Vel_alert_mmsRMS, Vel_danger_mmsRMS, 'BPFO'
+                    )
+                except Exception as e:
+                    BPFOdetection = False
+                    BPFOseverity = 0                                    
+                print("BPFOseverity",BPFOseverity)
+                if not BPFOdetection:
+                    BPFOStatusCode = 0
+                elif BPFOseverity < 6:
+                    BPFOStatusCode = 1
+                else:
+                    BPFOStatusCode = 2
+
+            except Exception as e:
+                BPFOStatusCode = 0
+            results["BPFO"]["status_codes"].append(BPFOStatusCode)
+
+            try:
+                """内圈故障诊断"""
+                print("内圈")       
+                fbr_BPFI = BPFI     
+                #fbr_BPFI = ShaftSpd * BPFI
+                fTarget = fbr_BPFI
+                if self.wind_code == 'WOF091200030':
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <1:
+                        fmax_gE_10 =10         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000     
+                else:
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <10:
+                        fmax_gE_10 =100         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000    
+                
+                numHarmonics = 4
+                fRangeMode = "Per"
+                Detection = "RMS"
+                print('BPFI',BPFI)    
+                print('ShaftSpd',ShaftSpd) 
+                print('fTarget',fTarget)  
+
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(self.VelSpec[i], fmax_vel, fTarget, fRange, \
+                                                        fRangeMode, Detection, numHarmonics)            
+                except Exception as e:
+                    #print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)   
+                print('Vel_RunningSpd_HM\n',Vel_RunningSpd_HM)   
+                
+                # 速度频谱特征
+                try:           
+                    Vel_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                   # print(f"计算 Vel_Bearing_Defect_HM 时出错: {e}")
+                    Vel_Bearing_Defect_HM = [0] * (numHarmonics + 1)  
+                print('Vel_Bearing_Defect_HM\n',Vel_Bearing_Defect_HM)   
+
+                # 包络频谱特征
+                try:
+                    Detection = "Peak"
+                    gE_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 gE_Bearing_Defect_HM 时出错: {e}")
+                    gE_Bearing_Defect_HM = [0] * (numHarmonics + 1)                 
+                print('gE_Bearing_Defect_HM\n',gE_Bearing_Defect_HM)     
+
+                try:
+                    # 诊断内圈故障
+                    BPFIdetection, BPFIseverity = self.myDiagnosisLib.Bearing_Fault_Diag(
+                        gE_Bearing_Defect_HM, Vel_Bearing_Defect_HM, Vel_RunningSpd_HM, gE_alert, gE_danger,
+                        Vel_alert_mmsRMS, Vel_danger_mmsRMS, 'BPFI'
+                    )
+                except Exception as e:
+                    BPFIdetection =False
+                    BPFIseverity =0
+                if not BPFIdetection:
+                    BPFIStatusCode = 0
+                elif BPFIseverity < 6:
+                    BPFIStatusCode = 1
+                else:
+                    BPFIStatusCode = 2
+            except Exception as e:  
+                #print(f"处理内圈故障时发生未捕获的异常: {e}")
+                BPFIStatusCode = 0
+            results["BPFI"]["status_codes"].append(BPFIStatusCode)
+
+            try:
+                """滚动体故障诊断"""
+                print("滚动体")      
+                fbr_BSF = BSF          
+                    
+                if self.wind_code == 'WOF091200030':
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <1:
+                        fmax_gE_10 =10         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000     
+                else:
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <10:
+                        fmax_gE_10 =100       
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000    
+            
+                numHarmonics = 4
+                fTarget = fbr_BSF
+                fRangeMode = "Per"
+                Detection = "RMS"
+                print('BSF',BSF)    
+                print('ShaftSpd',ShaftSpd) 
+                print('fTarget',fTarget)   
+
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(self.VelSpec[i], fmax_vel, fTarget, fRange, \
+                                                            fRangeMode, Detection, numHarmonics) 
+                except Exception as e:  
+                   # print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                      
+                print('Vel_RunningSpd_HM\n',Vel_RunningSpd_HM)    
+
+                try:        
+                    # 速度频谱特征
+                    Vel_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:  
+                    #print(f"计算 Vel_Bearing_Defect_HM 时出错: {e}")
+                    Vel_Bearing_Defect_HM = [0] * (numHarmonics + 1)         
+                print('Vel_Bearing_Defect_HM\n',Vel_Bearing_Defect_HM)  
+
+                try:
+                    # 包络频谱特征S
+                    Detection = "Peak"
+                    gE_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:  
+                    #print(f"计算 gE_Bearing_Defect_HM 时出错: {e}")
+                    gE_Bearing_Defect_HM = [0] * (numHarmonics + 1)   
+                print('gE_Bearing_Defect_HM\n',gE_Bearing_Defect_HM)  
+
+                try:
+                    # 诊断滚动体故障
+                    BSFdetection, BSFseverity = self.myDiagnosisLib.Bearing_Fault_Diag(
+                        gE_Bearing_Defect_HM, Vel_Bearing_Defect_HM, Vel_RunningSpd_HM, gE_alert, gE_danger,
+                        Vel_alert_mmsRMS, Vel_danger_mmsRMS, 'BSF'
+                    )
+                except Exception as e:  
+                    BSFdetection=False
+                    BSFseverity=0
+
+                if not BSFdetection:
+                    BSFStatusCode = 0
+                elif BSFseverity < 6:
+                    BSFStatusCode = 1
+                else:
+                    BSFStatusCode = 2
+            except Exception as e:  
+                #print(f"处理滚动体故障时发生未捕获的异常: {e}")
+                BSFStatusCode = 0
+            results["BSF"]["status_codes"].append(BSFStatusCode)
+
+            try:
+                """保持架故障诊断"""
+                print("保持架")
+                fbr_FTF =  FTF           
+                fTarget = fbr_FTF
+
+                if self.wind_code == 'WOF091200030':
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    elif fTarget <1:
+                        fmax_gE_10 =10         
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000     
+                else:
+                    if fTarget >100:
+                        fmax_gE_10 =10000  
+                    # elif fTarget <1:
+                    #     fmax_gE_10 =0.001       
+                    #     fmax_vel=min(1000,fs_acc/2.56)                    
+                    elif fTarget <10:
+                        fmax_gE_10 =10        
+                        fmax_vel=min(1000,fs_acc/2.56)
+                    else:
+                        fmax_gE_10 =1000    
+            
+                numHarmonics = 4 
+                fRangeMode = "Per"
+                Detection = "RMS"
+                print('FTF',FTF)    
+                print('ShaftSpd',ShaftSpd) 
+                print('fTarget',fTarget)   
+
+                try:
+                    Vel_RunningSpd_HM = self.myDiagnosisLib.FixedFreqFeature(self.VelSpec[i], fmax_vel, fTarget, fRange, \
+                                                            fRangeMode, Detection, numHarmonics)
+                except Exception as e:  
+                    #print(f"计算 Vel_RunningSpd_HM 时出错: {e}")
+                    Vel_RunningSpd_HM = [0] * (numHarmonics + 1)                  
+                print('Vel_RunningSpd_HM\n',Vel_RunningSpd_HM)          
+                
+                try:
+                    # 速度频谱特征
+                    Vel_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:  
+                    #print(f"计算 Vel_Bearing_Defect_HM 时出错: {e}")
+                    Vel_Bearing_Defect_HM = [0] * (numHarmonics + 1)   
+                print('Vel_Bearing_Defect_HM\n',Vel_Bearing_Defect_HM)   
+
+                try:
+                    # 包络频谱特征
+                    Detection = "Peak"
+                    gE_Bearing_Defect_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, fRangeMode, Detection, numHarmonics
+                    )
+                except Exception as e:  
+                    #print(f"计算 gE_Bearing_Defect_HM 时出错: {e}")
+                    gE_Bearing_Defect_HM = [0] * (numHarmonics + 1)   
+                print('gE_Bearing_Defect_HM\n',gE_Bearing_Defect_HM)   
+
+                try:
+                    # 诊断保持架故障
+                    FTFdetection, FTFseverity = self.myDiagnosisLib.Bearing_Fault_Diag(
+                        gE_Bearing_Defect_HM, Vel_Bearing_Defect_HM, Vel_RunningSpd_HM, gE_alert, gE_danger,
+                        Vel_alert_mmsRMS, Vel_danger_mmsRMS, 'FTF'
+                    )
+                except Exception as e: 
+                    FTFdetection =False
+                    FTFseverity =0
+                if not FTFdetection:
+                    FTFStatusCode = 0
+                elif FTFseverity < 6:
+                    FTFStatusCode = 1
+                else:
+                    FTFStatusCode = 2
+            except Exception as e:  
+                #print(f"处理保持架故障时发生未捕获的异常: {e}")
+                FTFStatusCode = 0
+            results["FTF"]["status_codes"].append(FTFStatusCode)
+
+        # 计算统计信息
+        status_counts = {"count_0": 0, "count_1": 0, "count_2": 0}
+        
+        for fault_type in results:
+            status_codes = results[fault_type]["status_codes"]
+            results[fault_type]["max_status"] = max(status_codes)
+            
+            # 累加状态计数
+            status_counts["count_0"] += status_codes.count(0)
+            status_counts["count_1"] += status_codes.count(1)
+            status_counts["count_2"] += status_codes.count(2)
+        
+        return {
+            "type": "Bearing",
+            "results": results,
+            "statistics": status_counts
+        }
+
+    def Gear_diag(self):
+        """齿轮诊断"""
+        print('齿轮诊断')
+        mesure_point_name = self.mesure_point_names[0].lower()
+        if 'gearbox' not in mesure_point_name:
+            print("Can not perform gearbox diagnosis")
+            raise ValueError("Can not perform gearbox diagnosis")
+
+        results = {
+            "wear": {"status_codes": []},
+            "crack": {"status_codes": []}
+        }
+        id =self.ids[0]
+        numTeeth = int(self.teeth_count.get(id,0))
+        print(numTeeth)
+        for i in range(len(self.ids)):
+            ShaftSpd =self.fspd[i] *self.speed_ratio
+            if ShaftSpd ==0:
+                results["wear"]["status_codes"].append(-1)                
+                results["crack"]["status_codes"].append(-1)
+
+            print('speed_ratio',self.speed_ratio)
+            print('ShaftSpd',ShaftSpd)
+            try:
+                GMF = ShaftSpd * numTeeth / 60
+                fTarget = GMF
+                print('fTarget',fTarget)
+                #numHarmonics = 2  
+                numHarmonics = 2             
+                fRange = 0.05                  
+                fmax_gE_10 =1000
+
+                fs_acc = self.fs_acc[i]
+                fmax_vel = max(1000,fs_acc/2.56)
+                if fTarget >100:
+                    fmax_gE_10 =10000  
+                elif fTarget <0.1:
+                    fmax_gE_10 =10         
+                elif fTarget <1:
+                    fmax_gE_10 =100
+                    fmax_vel=min(1000,fs_acc/2.56)
+                else:
+                    fmax_gE_10 =1000
+                fmax_acc =fs_acc /2.56   
+                try:
+                    Vel_GMF_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.VelSpec[i], fmax_vel, fTarget, fRange, "Per", "RMS", numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 Vel_GMF_HM 时出错: {e}")
+                    Vel_GMF_HM = [0] * (numHarmonics + 1)                       
+                print('Vel_GMF_HM',Vel_GMF_HM)      
+
+                try:      
+                    Acc_GMF_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.AccSpec[i], fmax_acc, fTarget,fRange, "Per", "Peak", numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 Acc_GMF_HM 时出错: {e}")
+                    Acc_GMF_HM = [0] * (numHarmonics + 1)                      
+                print('Acc_GMF_HM',Acc_GMF_HM) 
+
+                try:
+                    gE_GMF_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, "Per", "Peak", numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 gE_GMF_HM 时出错: {e}")
+                    gE_GMF_HM = [0] * (numHarmonics + 1)                      
+                print('gE_GMF_HM',gE_GMF_HM)  
+
+                fTarget =ShaftSpd /60
+                print('fTarget',fTarget)               
+                if fTarget >100:
+                    fmax_gE_10 =10000  
+                elif fTarget <0.1:
+                    fmax_gE_10 =10         
+                elif fTarget <1:
+                    fmax_gE_10 =100
+                    fmax_vel=min(1000,fs_acc/2.56)
+                else:
+                    fmax_gE_10 =1000           
+                numHarmonics = 4
+                try: 
+                    gE_GTIF_HM = self.myDiagnosisLib.FixedFreqFeature(
+                        self.gESpec_10[i], fmax_gE_10, fTarget, fRange, "Per", "Peak", numHarmonics
+                    )
+                except Exception as e:
+                    #print(f"计算 gE_GTIF_HM 时出错: {e}")
+                    gE_GTIF_HM = [0] * (numHarmonics + 1)                      
+                print('gE_GTIF_HM',gE_GTIF_HM)       
+
+                Acc_GMF_alert_gPk = 0.3  # 齿轮缺陷加速度预警阈值 (g 峰值)
+                Acc_GMF_danger_gPk = 0.5  # 齿轮缺陷加速度报警阈值 (g 峰值)
+                Vel_GMF_alert_mmsRMS = 4.0  # 速度总值预警阈值 （mm/s 有效值）
+                Vel_GMF_danger_mmsRMS = 7.0  # 速度总值报警阈值 （mm/s 有效值）
+                gE_GMF_alert_gE = 3.0  # 齿轮磨损缺陷加速度包络预警阈值 (gE 峰值)
+                gE_GMF_danger_gE = 10.0  # 齿轮磨损缺陷加速度包络报警阈值 (gE 峰值)
+                gE_GTIF_alert_gE = 3.0  # 齿轮断齿缺陷加速度包络预警阈值 (gE 峰值)
+                gE_GTIF_danger_gE = 10.0  # 齿轮断齿缺陷加速度包络报警阈值 (gE 峰值) 
+                try:
+                    # 诊断齿轮故障
+                    (GearToothWearDetection, GearToothWearSev, GearToothCrackBrokenDetection, GearToothCrackBrokenSev) = \
+                        self.myDiagnosisLib.Gear_Fault_Diag(Acc_GMF_HM, Vel_GMF_HM, gE_GMF_HM, gE_GTIF_HM, Acc_GMF_alert_gPk, \
+                                                    Acc_GMF_danger_gPk, Vel_GMF_alert_mmsRMS, Vel_GMF_danger_mmsRMS, gE_GMF_alert_gE,
+                                                    gE_GMF_danger_gE, \
+                                                    gE_GTIF_alert_gE, gE_GTIF_danger_gE)     
+                except Exception as e:
+                    GearToothWearDetection =False
+                    GearToothCrackBrokenDetection =False
+                    GearToothWearSev =0
+                    GearToothCrackBrokenSev =0               
+                print('齿轮磨损', GearToothWearSev, '齿轮断齿', GearToothCrackBrokenSev)            
+                # 处理磨损状态码
+                if not GearToothWearDetection:
+                    wear_status = 0
+                elif GearToothWearSev < 6:
+                    wear_status = 1
+                else:
+                    wear_status = 2
+
+                
+                # 处理裂纹状态码
+                if not GearToothCrackBrokenDetection:
+                    crack_status = 0
+                elif GearToothCrackBrokenSev < 6:
+                    crack_status = 1
+                else:
+                    crack_status = 2
+            except Exception as e:  
+                wear_status=0
+                crack_status=0    
+            results["wear"]["status_codes"].append(wear_status)                   
+            results["crack"]["status_codes"].append(crack_status)
+        
+        # 计算统计信息
+        status_counts = {"count_0": 0, "count_1": 0, "count_2": 0}
+        
+        for fault_type in results:
+            status_codes = results[fault_type]["status_codes"]
+            results[fault_type]["max_status"] = max(status_codes)
+            
+            # 累加状态计数
+            status_counts["count_0"] += status_codes.count(0)
+            status_counts["count_1"] += status_codes.count(1)
+            status_counts["count_2"] += status_codes.count(2)
+        
+        return {
+            "type": "Gear",
+            "results": results,
+            "statistics": status_counts
+        }
+    
+
+    def Characteristic_Frequency(self):
+        """提取轴承、齿轮等参数"""
+        #前端给的ids所对应的轴承参数、齿轮参数是一致的 计算一次即可
+        if not self.bearing_frequencies:  
+            id = self.ids[0]  # 取第一个id计算
+            mesure_point_name = self.mesure_point_names[0]
+            rpm = self.fspd[0] * self.speed_ratio  # 应用转速传动比
+	    
+            # 1、从测点名称中提取部件名称（计算特征频率的部件）
+            str1 = mesure_point_name
+            print(str1)
+
+            # 2、连接233的数据库'energy_show'，从表'wind_engine_group'查找风机编号'engine_code'对应的机型编号'mill_type_code'
+            engine_code = self.engine_code
+            print(engine_code)
+            Engine = create_engine('mysql+pymysql://admin:admin123456@192.168.50.233:3306/energy_show')
+            #Engine = create_engine('mysql+pymysql://admin:admin123456@106.120.102.238:16306/energy_show')        
+            df_sql2 = f"SELECT * FROM wind_engine_group WHERE engine_code = '{engine_code}'"
+            df2 = pd.read_sql(df_sql2, Engine)
+            mill_type_code = df2['mill_type_code'].iloc[0]
+            print(mill_type_code)
+
+            # 3、从相关的表中通过机型编号'mill_type_code'或者齿轮箱编号gearbox_code查找部件'brand'、'model'的参数信息
+            #unit_bearings主轴承参数表 关键词"main_bearing"
+            if 'main_bearing' in str1:
+                print("main_bearing")
+                # df_sql3 = f"SELECT * FROM {'unit_bearings'} where mill_type_code = {'mill_type_code'} "
+                df_sql3 = f"SELECT * FROM unit_bearings WHERE mill_type_code = '{mill_type_code}' "            
+                df3 = pd.read_sql(df_sql3, Engine)
+                if df3.empty:
+                    print("警告: 没有找到有效的机型信息")   
+                if 'front' in str1:
+                    brand = 'front_bearing' + '_brand'  
+                    model = 'front_bearing' + '_model'  
+                    front_has_value = not pd.isna(df3[brand].iloc[0]) and not pd.isna(df3[model].iloc[0]) 
+                    if not  front_has_value:
+                        print("警告: 没有找到有效的品牌信息")
+                elif 'rear' in str1:
+                    brand = 'rear_bearing' + '_brand'  
+                    model = 'rear_bearing' + '_model'  
+                    end_has_value = not pd.isna(df3[brand].iloc[0]) and not pd.isna(df3[model].iloc[0])       
+                    if not  end_has_value:
+                        print("警告: 没有找到有效的品牌信息")                
+                else:
+                    # 当没有指定 front 或 end 时，自动选择有值的轴承信息
+                    front_brand_col = 'front_bearing_brand'
+                    front_model_col = 'front_bearing_model'
+                    rear_brand_col = 'rear_bearing_brand'
+                    rear_model_col = 'rear_bearing_model'
+                    # 检查 front_bearing 是否有值
+                    front_has_value = not pd.isna(df3[front_brand_col].iloc[0]) and not pd.isna(df3[front_model_col].iloc[0])                
+                    # 检查 end_bearing 是否有值
+                    end_has_value = not pd.isna(df3[rear_brand_col].iloc[0]) and not pd.isna(df3[rear_model_col].iloc[0])               
+                    # 根据检查结果选择合适的列
+                    if front_has_value and end_has_value:
+                        # 如果两者都有值，默认选择 front
+                        brand = front_brand_col
+                        model = front_model_col
+                    elif front_has_value:
+                        brand = front_brand_col
+                        model = front_model_col
+                    elif end_has_value:
+                        brand = rear_brand_col
+                        model = rear_model_col
+                    else:
+                        # 如果两者都没有有效值，设置默认值或抛出异常
+                        print("警告: 没有找到有效的轴承信息")
+                        brand = front_brand_col  # 默认使用 front
+                        model = front_model_col  # 默认使用 front
+                print(brand)
+                _brand = df3[brand].iloc[0]
+                _model = df3[model].iloc[0]
+                print(_brand)
+                print(_model)     
+            #unit_dynamo 发电机参数表 关键词generator stator
+            elif 'generator'in str1 or 'stator' in str1:
+                print("generator or 'stator'")
+                # df_sql3 = f"SELECT * FROM {'unit_dynamo'} where mill_type_code = {'mill_type_code'} "
+                df_sql3 = f"SELECT * FROM unit_dynamo WHERE mill_type_code = '{mill_type_code}' "
+                df3 = pd.read_sql(df_sql3, Engine)
+                if 'non' in str1:
+                    brand = 'non_drive_end_bearing' + '_brand'  
+                    model = 'non_drive_end_bearing' + '_model'  
+                else:
+                    brand = 'drive_end_bearing' + '_brand'  
+                    model = 'drive_end_bearing' + '_model'              
+                print(brand)
+                _brand = df3[brand].iloc[0]
+                _model = df3[model].iloc[0]
+                print(_brand)
+                print(_model)       
+
+            #齿轮箱区分行星轮/平行轮 和 轴承两个表
+            elif 'gearbox' in str1:
+                print("gearbox")
+                #根据mill_type_code从unit_gearbox表中获得gearbox_code
+                df_sql3 = f"SELECT * FROM unit_gearbox WHERE mill_type_code = '{mill_type_code}' "
+                df3 = pd.read_sql(df_sql3, Engine)
+                gearbox_code =df3['code'].iloc[0]
+                print(gearbox_code)
+                #如果是行星轮/平行轮 则从unit_gearbox_structure 表中取数据
+                if 'planet'in str1 or 'sun' in str1:
+                    print("'planet' or 'sun' ")
+                    gearbox_structure =1 if 'planet'in str1 else 2
+                    planetary_gear_grade =1
+                    if 'first' in str1:
+                        planetary_gear_grade =1
+                    elif 'second'in str1:
+                        planetary_gear_grade =2
+                    elif 'third'in str1:
+                        planetary_gear_grade =3
+                    # df_sql33 = f"SELECT * FROM unit_gearbox_structure WHERE gearbox_code = '{gearbox_code}' "
+                    df_sql33 = f"""
+                    SELECT bearing_brand, bearing_model,gear_ring_teeth_count
+                    FROM unit_gearbox_structure 
+                        WHERE gearbox_code = '{gearbox_code}' 
+                        AND gearbox_structure = '{gearbox_structure}'
+                        AND planetary_gear_grade = '{planetary_gear_grade}'
+                        """
+                    df33 = pd.read_sql(df_sql33, Engine)
+                    if df33.empty:
+                        print("unit_gearbox_structure没有该测点的参数")
+                    else:
+                        brand = 'bearing' + '_brand'  
+                        model = 'bearing' + '_model'  
+                        teeth_count= df33['gear_ring_teeth_count'].iloc[0]
+                        print('teeth_count',teeth_count)                           
+                        print(brand)
+                        _brand = df33[brand].iloc[0]
+                        _model = df33[model].iloc[0]
+                        has_value = not pd.isna(df33[brand].iloc[0]) and not pd.isna(df33[model].iloc[0])
+                        if has_value:   
+                            print(_brand)
+                            print(_model)   
+                        else:
+                            print("警告: 没有找到有效的轴承信息")
+                #如果是齿轮箱轴承 则从unit_gearbox_bearings 表中取数据
+                elif 'shaft' in str1 or'input' in str1:
+                    print("'shaft'or'input'")
+                    #高速轴 低速中间轴 取bearing_rs/gs均可
+                    parallel_wheel_grade=1
+                    if 'low_speed' in str1:
+                        parallel_wheel_grade =1                
+                    elif 'low_speed_intermediate' in str1:
+                        parallel_wheel_grade =2
+                    elif 'high_speed' in str1:
+                        parallel_wheel_grade =3   
+                    # df_sql33 = f"SELECT * FROM unit_gearbox_bearings WHERE gearbox_code = '{gearbox_code}' "
+                    df_sql33 = f"""
+                    SELECT bearing_rs_brand, bearing_rs_model, bearing_gs_brand, bearing_gs_model,gear_ring_teeth_count
+                    FROM unit_gearbox_bearings 
+                    WHERE gearbox_code = '{gearbox_code}'
+                    AND parallel_wheel_grade = '{parallel_wheel_grade}'
+                    """
+                    df33 = pd.read_sql(df_sql33, Engine)
+                    if not df33.empty:
+                        if 'high_speed' in str1 or 'low_speed_intermediate' in str1:
+                            rs_brand = 'bearing_rs' + '_brand'  
+                            rs_model = 'bearing_rs' + '_model'  
+                            gs_brand = 'bearing_gs' + '_brand'  
+                            gs_model = 'bearing_gs' + '_model'  
+                            rs_has_value = not pd.isna(df33[rs_brand].iloc[0]) and not pd.isna(df33[rs_model].iloc[0])     
+                            gs_has_value = not pd.isna(df33[gs_brand].iloc[0]) and not pd.isna(df33[gs_model].iloc[0])    
+                            if rs_has_value and gs_has_value:
+                                brand = rs_brand
+                                model = rs_model
+                            elif rs_has_value:
+                                brand = rs_brand
+                                model = rs_model
+                            elif gs_has_value:
+                                brand = gs_brand
+                                model = gs_model
+                            else:
+                                print("警告: 没有找到有效的品牌信息")
+                                brand = rs_brand 
+                                model = rs_model 
+                        #低速轴 取bearing_model          
+                        elif 'low_speed'in str1:
+                            brand = 'bearing' + '_brand'  
+                            model = 'bearing' + '_model'
+                    else:
+                        print("警告: 没有找到有效的轴承信息")     
+                    if  not df33.empty:
+                        if 'high_speed' in str1:
+                            teeth_count= df33['gear_ring_teeth_count'].iloc[0]
+                            print('teeth_count',teeth_count)                     
+                    print(brand)
+                    _brand = df33[brand].iloc[0]
+                    _model = df33[model].iloc[0]
+                    print(_brand)
+                    print(_model)   
+            # 4、从表'unit_dict_brand_model'中通过'_brand'、'_model'查找部件的参数信息
+            df_sql4 = f"SELECT * FROM unit_dict_brand_model where manufacture = %s AND model_number = %s"
+            params = [(_brand, _model)]
+            df4 = pd.read_sql(df_sql4, Engine, params=params)
+            n_rolls = df4['rolls_number'].iloc[0]
+            d_rolls = df4['rolls_diameter'].iloc[0]
+            D_diameter = df4['circle_diameter'].iloc[0]
+            theta_deg = df4['theta_deg'].iloc[0]
+
+            # 计算特征频率
+            bearing_freq = self.calculate_bearing_frequencies(n_rolls, d_rolls, D_diameter, theta_deg, rpm)
+            print(bearing_freq)
+            # 将结果赋给所有id
+            
+            for id in self.ids:
+                self.bearing_frequencies[id] = bearing_freq
+                print('aaaa')
+                if 'gearbox' in str1:
+                    self.teeth_count[id] =teeth_count
+    
+
+
+    def calculate_bearing_frequencies(self, n, d, D, theta_deg, rpm):
+        """
+        计算轴承各部件特征频率
+
+        参数：
+        n (int): 滚动体数量
+        d (float): 滚动体直径（单位：mm）
+        D (float): 轴承节圆直径（滚动体中心圆直径，单位：mm）
+        theta_deg (float): 接触角（单位：度）
+        rpm (float): 转速（转/分钟）
+
+        返回：
+        dict: 包含各特征频率的字典（单位：Hz）
+        """
+        # 转换角度为弧度
+        theta = math.radians(theta_deg)
+
+        # 转换直径单位为米（保持单位一致性，实际计算中比值抵消单位影响）
+        # 注意：由于公式中使用的是比值，单位可以保持mm不需要转换
+        ratio = d / D
+
+        # 基础频率计算（转/秒）
+        # f_r = rpm 
+        f_r = rpm / 60.0
+
+        # 计算各特征频率
+        BPFI = (n / 2) * (1 + ratio * math.cos(theta)) * f_r  # 内圈故障频率
+        BPFO = (n / 2) * (1 - ratio * math.cos(theta)) * f_r  # 外圈故障频率
+        BSF = (D / (2 * d)) * (1 - (ratio ** 2) * (math.cos(theta) ** 2)) * f_r  # 滚动体故障频率
+        FTF = 0.5 * (1 - ratio * math.cos(theta)) * f_r  # 保持架故障频率
+
+        return {
+            "BPFI": round(BPFI, 2),
+            "BPFO": round(BPFO, 2),
+            "BSF": round(BSF, 2),
+            "FTF": round(FTF, 2),
+
+        }
+	
+    def calculate_speed_ratio(self):
+        """
+        根据风场编号和测点名称计算转速传动比
+        """
+        # 默认传动比为1
+        self.speed_ratio = 1.0
+        
+        # WOF046400029风场 七台河需要特殊处理
+        if self.wind_code == "WOF046400029":
+            # 获取第一个测点名称
+            mesure_point_name = self.mesure_point_names[0]
+
+            if "gearbox" in mesure_point_name:
+                if "high_speed" in mesure_point_name:
+                    # 高速轴输出端
+                    self.speed_ratio = 1 
+                elif "first" in mesure_point_name:
+                    # 一级行星级
+                    self.speed_ratio = 1 /(5.3 * 5.5 * 3.8)	
+                elif "second" in mesure_point_name:
+                    # 二级行星级
+                    self.speed_ratio = 1 /  (5.5 * 3.8)
+            if "generator" in mesure_point_name:
+                self.speed_ratio =1
+            else:
+                # 非齿轮箱测点
+                self.speed_ratio = 1 / (5.3 * 5.5 * 3.8)	
+	
+	
+	
+	
+