1.分析结果
包络图是第三行。第一行水平时域信号，第二行垂直时域信号，第三行对第1行信号在故障频点包络计算后的结果。
上方是最终的视图。右下角我把光标对准了低频位置的那个故障峰，与理论值是一致的.
2.源码
2.1 测试例程源码
import gp_xjtu_data
import gp_calc_envelope
import gp_xjtu_toolsdef test_gp_calc_envelope():(x, sig_h, sig_v, fullpathOfFile) = gp_xjtu_data.get_demo_data();descToDraw = []descToDraw.append(x);descToDraw.append(sig_h);descToDraw.append("Time Domain Signal (XJTU-SY Bearing_Datasets 35Hz12KN 1-1)")descToDraw.append("Time(s) 256000sa/s, 32768pt.")descToDraw.append("g(acc)");descToDraw.append(fullpathOfFile + " horizon");signals = []signals.append(descToDraw.copy());descToDraw = []descToDraw.append(x);descToDraw.append(sig_v);descToDraw.append("Time Domain Signal (XJTU-SY Bearing_Datasets 35Hz12KN 1-1)")descToDraw.append("Time(s) 256000sa/s, 32768pt.")descToDraw.append("g(acc)");descToDraw.append(fullpathOfFile + " horizon");signals.append(descToDraw.copy());(signal_final, signal_demodulated, signal_filtered, signal_ac) = gp_calc_envelope.calc_envelope_of(x, sig_h, 180)descToDraw = []descToDraw.append(x);descToDraw.append(signal_final);descToDraw.append("Time Domain Signal - Envelope Wave")descToDraw.append("Time(s) 256000sa/s, 32768pt.")descToDraw.append("g(acc)");descToDraw.append(fullpathOfFile + " horizon");signals.append(descToDraw.copy());gp_xjtu_tools.DrawSignals(signals,1000, fullpathOfFile)test_gp_calc_envelope()
2.2 包络计算单元
import csv
import numpy as np
import matplotlib.pyplot as plt
import gp_xjtu_toolsdef calc_envelope_with_band(x, signal, freq_low, freq_high):signal_ac = gp_xjtu_tools.dc2ac(signal)signal_filtered = gp_xjtu_tools.band_filter(x, signal_ac, freq_low, freq_high)signal_demodulated = gp_xjtu_tools.envelope(signal_filtered)signal_final = gp_xjtu_tools.dc2ac(signal_demodulated)return (signal_final, signal_demodulated, signal_filtered, signal_ac)def calc_envelope_of(x, signal, centerFreq):return calc_envelope_with_band(x, signal, centerFreq*0.8, centerFreq*5)
2.3 数据获取单元
import csv
import numpy as np
import matplotlib.pyplot as pltdef ReadOneChannelRecordFromXJBD(file_path):horizontal_signals = []vertical_signals = []with open(file_path, newline='') as csvfile:reader = csv.reader(csvfile)next(reader) # 跳过标题行for row in reader:horizontal_signals.append(float(row[0]))vertical_signals.append(float(row[1]))return (horizontal_signals, vertical_signals)#将[+-50g的信号归一化为0点在32768处的u16信号]
def double2u16(data_in_array):data = np.array(data_in_array);normalized_data = np.clip(np.array(((data - (-50.0)) / 100.0 * 65535), dtype=np.uint16), 0, 65535)return normalized_data;def resampleData(data_in_array, target_point):# 确定抽样点位置的索引值idx = np.linspace(0, len(data_in_array)-1, target_point, dtype=np.int32)# 使用interp函数进行抽样resampled_data = np.interp(idx, np.arange(len(data_in_array)), data_in_array)return resampled_data;def getSampledataIn2048U16(file_path, dataType):(hdata, vdata) = ReadOneChannelRecordFromXJBD(file_path);if(dataType == "h"):data = resampleData(hdata, 2048)else:data = resampleData(vdata, 2048)return double2u16(data);def get_demo_data():file_path=".\\data\\1-1\\"file_name="123.csv"full_path=file_path + file_name#full_path = r"D:\git2024\20240304_PushVibrationData\XJTU-SY_Bearing_Datasets\35Hz12kN\Bearing1_1\1.csv"(hSig, vSig) = ReadOneChannelRecordFromXJBD(full_path);#hSig = xjtu_file.double2u16(hSig);#vSig = xjtu_file.double2u16(vSig);x = np.arange(0,len(hSig));saps=25.6e3;x_axis = x*1.0/saps;return (x_axis, hSig, vSig, full_path)
2.4 基础工具集 - 计算 - 绘图
import csv
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt, hilbertdef dc2ac(signal):y = signalfft_signal = np.fft.fft(y)fft_signal[0] = 0; #dc=>acfft_signal_no_phase_correction = np.fft.ifft(fft_signal)return fft_signal_no_phase_correction;def band_filter(x, signal, freq_low, freq_high):saps=1.0/x[1]print('saps=', saps)# 定义一个二阶巴特沃斯滤波器b, a = butter(2, [2*freq_low/saps, 2*freq_high/saps], btype='band')# 应用滤波器filtered_signal = filtfilt(b, a, signal)signal_abs = [np.abs(z.real) for z in filtered_signal]return signal_absdef envelope(signal):complex_signal = hilbert(signal);signal_abs = [np.abs(z.real) for z in complex_signal]return signal_abs;def GetFFTOfSignal(x_axis, sig):#calc_fftfft_signal = np.fft.fft(y)fft_signal[0] = 0; #dc=>acfreq = np.arange(len(fft_signal)//2)freq = freq *(1.0/(len(fft_signal)*x[1]));print(freq[-1], len(freq))fft_signal_abs = fft_signal[0:len(fft_signal)//2]fft_signal_abs = [np.abs(z.real*x[1]) for z in fft_signal_abs]fft_x = freq;fft_y = np.abs(fft_signal_abs);ffty_toshow = fft_yffty_raw_ac = fft_signalreturn (fft_x, ffty_toshow, ffty_raw_ac);def GetFFTOfSignal_withShutdownFreq(x_axis, sig, freq_shutdown):#calc_fftfft_signal = np.fft.fft(y)fft_signal[0] = 0; #dc=>acfreq = np.arange(len(fft_signal)//2)freq = freq *(1.0/(len(fft_signal)*x[1]));print(freq[-1], len(freq))fft_signal_abs = fft_signal[0:len(fft_signal)//2]fft_signal_abs = [np.abs(z.real*x[1]) for z in fft_signal_abs]max_freq_index =next((i for i, num in enumerate(freq) if num > freq_shutdown), len(freq))print("max_freq_index = ", max_freq_index, ",values=", freq[max_freq_index])fft_x = freq[:max_freq_index];fft_y = np.abs(fft_signal_abs[:max_freq_index]);ffty_toshow = fft_yffty_raw_ac = fft_signalreturn (fft_x, ffty_toshow, ffty_raw_ac);def DrawSignals(signals, freq_shutdown, signal_memo):cnt = len(signals)# 绘制频谱图plt.figure(figsize=(12, 6))row = cnt;col = 2;sn = 1;for item in signals:x = item[0];y = item[1];sTitle = item[2];sX = item[3];sY = item[4];# 绘制时间域波形plt.subplot(row, 2, sn)sn = sn +1;plt.plot(x, y, color="green")plt.xlabel(sX)plt.ylabel(sY)plt.title(sTitle)#calc_fftfft_signal = np.fft.fft(y)fft_signal[0] = 0; #dc=>acfreq = np.arange(len(fft_signal)//2)freq = freq *(1.0/(len(fft_signal)*x[1]));print(freq[-1], len(freq))fft_signal_abs = fft_signal[0:len(fft_signal)//2]fft_signal_abs = [np.abs(z.real*x[1]) for z in fft_signal_abs]max_freq_index =next((i for i, num in enumerate(freq) if num > freq_shutdown), len(freq))print("max_freq_index = ", max_freq_index, ",values=", freq[max_freq_index])fft_x = freq[:max_freq_index];fft_y = np.abs(fft_signal_abs[:max_freq_index]);#绘制频谱plt.subplot(row, 2, sn)sn = sn +1plt.plot(fft_x, fft_y, color='blue')str_fft = 'Freq Spectrum of Left Figure(datafile=%s)' %(signal_memo)plt.title(str_fft)plt.xlabel('Frequency[cut to %5.1fHz], totalPt=%d' %(freq_shutdown, len(fft_signal_abs)))plt.ylabel('Amplitude')plt.axis('auto') # 启用自动调整大小和缩放功能plt.grid(True)plt.tight_layout()plt.show()
