import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from math import pi
import math
from matplotlib.font_manager import FontProperties
ChineseFont1 = FontProperties(fname = 'C:\\Windows\\Fonts\\simsun.ttc')
from DiagnosisLib import DiagnosisLib

myDiagnosisLib = DiagnosisLib()

###############################################################################
# 1、提取原始振动加速度数据并确定参数
'''
# Acc = pd.read_csv('中国风电合肥莲花山风电_14#风机_主轴前轴承水平1H_128k 加速度波形(0.1-10000)_25600HZ_加速度_1100RPM_20181203014901.txt')
with open("中国风电合肥莲花山风电_14#风机_主轴前轴承水平1H_128k 加速度波形(0.1-10000)_25600HZ_加速度_1100RPM_20181203014901.txt", "r") as file:  #使用绝对路径 开文件
  Acc = file.read()
Acc = np.fromstring(Acc, dtype=float, sep=' ')   # 原始振动加速度时域波形
'''
Acc = pd.read_csv('国能山东--长清风电场_A1_齿轮箱中间轴轴向5A_25600_0.0_20240513031500.csv', header=None)
Acc = Acc.values
Acc = Acc.reshape(-1)     # 转化为一维数组
# 原始数据参数
fspd = 10
ShaftSpd = fspd
fs_acc = 25600
fmax_acc = 10000
t_acc = np.arange(0, 1/fs_acc*len(Acc), 1/fs_acc)

###############################################################################
# 2、计算振动加速度、振动速度、振动加速度包络三种频谱
# 2.1、计算振动加速度频谱
[fAcc, AccSpec] = myDiagnosisLib.Spectrum(Acc,fs_acc)

# 2.2、计算振动速度频谱
fCut = 0.1
fMax = 10000
VelSpec = myDiagnosisLib.FreqIntegration(AccSpec, fMax, fCut)

# 2.3、计算加速度包络频谱
# （1）振动加速度波形做带通滤波(50-1000Hz)
lowcut = 500
highcut = 10000
fs = fs_acc
order = 6
Acc_Band = myDiagnosisLib.BandPassButter(Acc, lowcut, highcut, fs, order)

# （2）带通滤波后的波形做包络(带通范围50-1000Hz)
DownSampleRate = 8   # 降采样为1/8
gE_10 = myDiagnosisLib.EnvLopInTime(Acc_Band, DownSampleRate, Method = "SKF")   # 加速度包络时域波形
t_gE_10 =t_acc[::DownSampleRate]

# （3）对包络波形做频谱(分析频率范围0-100Hz)
fmax_gE_10 = 1000   # 包络频谱的最大分析频率
fs_gE_10 = 2.56*fmax_gE_10
[fAcc_gE, gESpec_10] = myDiagnosisLib.Spectrum(gE_10,fs_gE_10)

###############################################################################

fRangeMode = "百分比"
fRange = 0.05

###############################################################################
# 诊断跟设备旋转频率相关的故障，即：平衡，对中，松动
fTarget = ShaftSpd
numHarmonics = 4
fmax_vel = 1000

Detection = "有效值"
Vel_RunningSpd_HM = []
Vel_RunningSpd_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                        fRangeMode, Detection, numHarmonics)
# Vel_RunningSpd_HM，生成转频的1X-5X的倍频值
# VelSpec，振动速度频谱
# fmax_vel，最大分析频率
# fTarget，目标频率为转频
Vel_alert_mmsRMS = 4                                                                   # 速度总值预警阈值 （mm/s 有效值）
Vel_danger_mmsRMS = 7                                                                  # 速度总值报警阈值 （mm/s 有效值）

defectType = 'Unbalance'
[UnbalanceDetection, UnbalanceSev] = myDiagnosisLib.RunningSpd_Fault_Diag(Vel_RunningSpd_HM, \
                                        Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType)
print ('Unbalance', UnbalanceDetection, UnbalanceSev)

defectType = 'Misalignment'
[MisalignmentDetction, MisalignmentSev] = myDiagnosisLib.RunningSpd_Fault_Diag(Vel_RunningSpd_HM, \
                                       Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType)
print ('Misalignment', MisalignmentDetction, MisalignmentSev)

defectType = 'Looseness'
[LoosenessDetection, LoosenessSev] = myDiagnosisLib.RunningSpd_Fault_Diag(Vel_RunningSpd_HM, \
                                      Vel_alert_mmsRMS, Vel_danger_mmsRMS, defectType)
print ('Looseness',LoosenessDetection, LoosenessSev)

###############################################################################
#3、诊断轴承相关的故障，即：外圈，内圈，滚动体，保持架
# 3.1、外圈故障诊断
BPFO = 4.63
fbr_BPFO = ShaftSpd*BPFO
fTarget = fbr_BPFO
numHarmonics = 4

Detection = "有效值"
Vel_Bearing_Defect_HM = []
Vel_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                                        fRangeMode, Detection, numHarmonics)
fmax_gE_10 = 1000

Detection = "峰值"
gE_Bearing_Defect_HM = []
gE_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                                       fRangeMode, Detection, numHarmonics)

Vel_alert_mmsRMS = 4                                                            # 速度总值预警阈值 （mm/s 有效值）
Vel_danger_mmsRMS = 7                                                           # 速度总值报警阈值 （mm/s 有效值）

gE_alert = 3                                                                    # 包络总值预警阈值（gE）
gE_danger = 10                                                                  # 包络总值报警阈值（gE）

defectType = 'BPFO'
#defectType = 'BPFI'
[BrDefDetection, BrDefSev] = 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, defectType)
print ('BPFO', BrDefDetection, BrDefSev)

# 3.2、内圈故障诊断
BPFI = 6.32
fbr_BPFI = ShaftSpd*BPFI
fTarget = fbr_BPFI

Detection = "有效值"
Vel_Bearing_Defect_HM = []
Vel_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                                        fRangeMode, Detection, numHarmonics)
fmax_gE_10 = 1000

Detection = "峰值"
gE_Bearing_Defect_HM = []
gE_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                                       fRangeMode, Detection, numHarmonics)

defectType = 'BPFI'
[BrDefDetection, BrDefSev] = 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, defectType)
print ('BPFI', BrDefDetection, BrDefSev)

# 3.3、滚动体故障诊断
BSF = 1.78
fbr_BSF = ShaftSpd*BSF
fTarget = fbr_BSF

Detection = "有效值"
Vel_Bearing_Defect_HM = []
Vel_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                                        fRangeMode, Detection, numHarmonics)
fmax_gE_10 = 1000

Detection = "峰值"
gE_Bearing_Defect_HM = []
gE_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                                       fRangeMode, Detection, numHarmonics)

defectType = 'BSF'
[BrDefDetection, BrDefSev] = 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, defectType)
print ('BSF', BrDefDetection, BrDefSev)

# 3.4、保持架故障诊断
FTF = 0.38
fbr_FTF = ShaftSpd*FTF
fTarget = fbr_FTF

fmax_gE_10 = 1000

Detection = "峰值"
gE_Bearing_Defect_HM = []
gE_Bearing_Defect_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                                       fRangeMode, Detection, numHarmonics)

defectType = 'FTF'
[BrDefDetection, BrDefSev] = 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, defectType)
print ('FTF', BrDefDetection, BrDefSev)

###############################################################################
# 4、诊断齿轮相关的故障（齿轮磨损和齿轮断齿）
numTeeth = 21                                                                   # 齿轮齿数
GMF = ShaftSpd*numTeeth
fTarget = GMF
numHarmonics = 2

Detection = "有效值"
Vel_GMF_HM = []
Vel_GMF_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                             fRangeMode, Detection, numHarmonics)

Detection = "峰值"
Acc_GMF_HM = []
Acc_GMF_HM = myDiagnosisLib.FixedFreqFeature(AccSpec, fmax_acc, fTarget, fRange, \
                                             fRangeMode, Detection, numHarmonics)

Detection = "峰值"
gE_GMF_HM = []
gE_GMF_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                            fRangeMode, Detection, numHarmonics)

fTarget = ShaftSpd
numHarmonics = 4
Detection = "峰值"
gE_GTIF_HM = []
gE_GTIF_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                             fRangeMode, Detection, numHarmonics)

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 峰值)

[GearToothWearDetection, GearToothWearSev, GearToothCrackBrokenDetection, GearToothCrackBrokenSev] = \
    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)
print('齿轮磨损', GearToothWearDetection, GearToothWearSev, '齿轮断齿', GearToothCrackBrokenDetection, \
      GearToothCrackBrokenSev)

###############################################################################
# 5、诊断电机转子条故障
numRotorBar = 44                                                                # 电机转子条数
RBPF = ShaftSpd*numRotorBar                                                     # 电机转子条通过频率
fTarget = RBPF
numHarmonics = 0

Detection = "有效值"
Vel_RBPF_HM = []
Vel_RBPF_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                              fRangeMode, Detection, numHarmonics)

Detection = "峰值"
Acc_RBPF_HM = []
Acc_RBPF_HM = myDiagnosisLib.FixedFreqFeature(AccSpec, fmax_acc, fTarget, fRange, \
                                              fRangeMode, Detection, numHarmonics)

Vel_RBPF_alert_mmsRMS = 4.0
Vel_RBPF_danger_mmsRMS = 7.0
Acc_RBPF_alert_gPk = 2.5
Acc_RBPF_danger_gPk = 3.75

[CrackedRotorbarDetection, CrackedRotorbarSev] = \
    myDiagnosisLib.Cracked_Rotorbar_Diag(Vel_RBPF_HM, Acc_RBPF_HM, Vel_RBPF_alert_mmsRMS, \
                          Vel_RBPF_danger_mmsRMS, Acc_RBPF_alert_gPk, Acc_RBPF_danger_gPk)

print('电机转子条故障',CrackedRotorbarDetection, CrackedRotorbarSev)

###############################################################################
# 6、诊断电机气隙不均匀故障
LF = 50
fTarget = 2*LF

numHarmonics = 0
Detection = "有效值"
Vel_2XLF_HM = []
Vel_2XLF_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                              fRangeMode, Detection, numHarmonics)

numHarmonics = 4
Detection = "峰值"
gE_2XLF_HM = []
gE_2XLF_HM = myDiagnosisLib.FixedFreqFeature(gESpec_10, fmax_gE_10, fTarget, fRange, \
                                             fRangeMode, Detection, numHarmonics)

Vel_2XLF_alert_mmsRMS = 4.0
Vel_2XLF_danger_mmsRMS = 7.0
gE_2XLF_alert_gE = 3
gE_2XLF_danger_gE = 10
[VariableAirgapDetection, VariableAirgapSev] = \
    myDiagnosisLib.Variable_AirGap_Diag(Vel_2XLF_HM, gE_2XLF_HM, \
            Vel_2XLF_alert_mmsRMS, Vel_2XLF_danger_mmsRMS, gE_2XLF_alert_gE, gE_2XLF_danger_gE)
print('电机气隙不均匀故障',VariableAirgapDetection, VariableAirgapSev)

###############################################################################
# 7、诊断离心泵叶片故障
numVane = 13
VPF = ShaftSpd*numVane
fTarget = VPF

numHarmonics = 1
Detection = "有效值"
Vel_VPF_HM = []
Vel_VPF_HM = myDiagnosisLib.FixedFreqFeature(VelSpec, fmax_vel, fTarget, fRange, \
                                              fRangeMode, Detection, numHarmonics)

Vel_VPF_alert_mmsRMS = 4.0
Vel_VPF_danger_mmsRMS = 7.0

[VaneWearDamageDetection, VaneWearDamageSev] = \
    myDiagnosisLib.Vane_WearDamage_Diag(Vel_VPF_HM, Vel_VPF_alert_mmsRMS, \
                                         Vel_VPF_danger_mmsRMS)
print('离心泵叶片故障',VaneWearDamageDetection, VaneWearDamageSev)