import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft, ifft #FFT库
import xlrd #读取excel库
# 打开采样数据文件
file_name = r'C:\Users\Chaver\Desktop\task\sampled.xlsx'
wb = xlrd.open_workbook(file_name)
# 根据 sheet 名称获取内容
sh1 = wb.sheet_by_name('3V-R-65Hz-data1')
# 读取采集A、B、C 线电压值 (A:Uab,B:Ubc,C:Uca )
phase_A = sh1.col_values(1,1,257) # A相电压
phase_B = sh1.col_values(2,1,257) # B相电压
# 上零点函数
def find_up(data):
up = []
for x in range(2,len(data)-1):
if data[x-2] < 50 and data[x-1] < 50 and data[x] < 50 and data[x+1] - data[x] > 50:
up.append(x)
return up
# 下零点函数
def find_dw(data):
dw = []
for x in range(1,len(data)-2):
if data[x+2] < 50 and data[x+1] < 50 and data[x] < 50 and data[x] - data[x-1] < -40:
dw.append(x)
return dw
# 获取相位差
def angle_error(phase_A,phase_B):
# 加窗 Hanning
phase_A_window = np.hanning(len(phase_A))*phase_A
phase_B_window = np.hanning(len(phase_B))*phase_B
# fft计算
F_phase_A = fft(phase_A_window)
F_phase_B = fft(phase_B_window)
# 角度计算
angle_A = np.angle(F_phase_A,deg=True)
angle_B = np.angle(F_phase_B,deg=True)
# 相位差
Delta_BA = np.subtract(angle_B,angle_A)
Delta_BA = [x+360 if x<-180 else 360-x if x>180 else x for x in Delta_BA]
return Delta_BA
# 获取还原前角度差
angle_original = angle_error(phase_A, phase_B)
print("还原前角度差:", angle_original[1:4])
# 获取A,B零点位置
up_A = find_up(phase_A)
dw_A = find_dw(phase_A)
up_B = find_up(phase_B)
dw_B = find_dw(phase_B)
# 还原全波
# 原始数据近零数值置零
def pre_data(data):
for x in range(len(data)):
if data[x] < 50:
data[x] = 0
return data
# 获取完整的第一个半波
def zero(D, up, dw):
for x in range(len(D)):
if x < up or x > dw:
D[x] = 0
return D
# 还原全波主函数
def wave_recovery(data_1,up,dw):
if up[0] < dw[0]:
half_length = dw[0] - up[0]
# 生成新列表data_2用来存放平移反转后的列表data_1
data_2 = [0] * dw[0] + [None] * (len(data_1) - dw[0])
for x in range(dw[0], len(data_1)):
data_2[x] = data_1[x - half_length]
data_2[x] = -data_2[x]
# 上零点左侧曲线复原,data_3
data_3 = [0] * up[0]
for x in range(up[0], 2 * up[0]):
data_3[x - up[0]] = -data_1[x]
# 将元素逆序呈现
data_3 = data_3[::-1]
# 补齐列表长度,添加零元素
data_3 = data_3 + [0] * (len(data_1) - len(data_3))
# 合并data_1与data_2
sum_1 = np.add(data_1, data_2)
phase_A = np.add(sum_1, data_3)
else:
half_length = up[0] - dw[0]
data_2 = [0] * dw[1] + [None] * (len(data_1) - dw[1])
for x in range(dw[1], len(data_1)):
data_2[x] = data_1[x - half_length]
data_2[x] = -data_2[x]
# 合并data_1和data_2,不知为什么一定要放在data_3的for循环前?
sum_data_1_2 = np.add(data_1, data_2)
# 获取第一个完整半波
data_3 = zero(data_1, up[0], dw[1])
# 将第一个半波向左平移half_length反转,得到data_4
data_4 = [0] * dw[0] + [0] * (len(data_3) - dw[0])
for x in range(dw[0], up[0]):
data_4[x] = data_3[x + half_length]
data_4[x] = -data_4[x]
# 合并data_1与data_2,data_4
phase_A = np.add(sum_data_1_2, data_4)
return phase_A
# 获取data1
phase_A_pre = pre_data(phase_A)
phase_B_pre = pre_data(phase_B)
# 获取还原后数据
phase_A_recovery=wave_recovery(phase_A_pre,up_A,dw_A)
phase_B_recovery=wave_recovery(phase_B_pre,up_B,dw_B)
plt.subplot(2,1,1)
plt.plot(phase_A_recovery,marker='o')
plt.subplot(2,1,2)
plt.plot(phase_B_recovery,marker='o')
# 比较还原前与还原后相位差
angle_recovery = angle_error(phase_A_recovery,phase_B_recovery)
print("还原后角度差: ", angle_recovery[1:4])
plt.show()
半波还原整波——封装函数
最新推荐文章于 2024-08-24 11:00:49 发布