#!/usr/bin/env python
# coding: utf-8
# In[1]:
import numpy as np # 导入NumPy库，用于数值计算
import pandas as pd # 导入Pandas库，用于数据处理
import matplotlib.pyplot as plt # 导入Matplotlib库，用于绘图
import pywt # 导入PyWavelets库，用于小波变换
from scipy.signal import welch # 从Scipy的信号模块中导入welch函数，用于计算功率谱密度
from datetime import datetime # 导入datetime模块中的datetime类，用于日期和时间处理
#更多模型咸鱼搜索机器学习之心，支持模型定制
# 从CSV文件读取数据
#data = pd.read_csv("A.csv") # 从CSV文件中读取数据，并将其存储在名为data的DataFrame中
# 读取Excel文件
data = pd.read_excel('A.xlsx')
# 将“time”列转为字符串类型
data["time"] = data["time"].astype(str)#使用csv文件时不需要
#更多模型咸鱼搜索机器学习之心，支持模型定制
t = data["time"].apply(lambda x: (datetime.strptime(x, '%Y-%m-%d %H:%M:%S') - datetime(1970, 1, 1)).total_seconds()).values#时间是用-分开的用这个
#t = data["time"].apply(lambda x: (datetime.strptime(x, '%Y/%m/%d %H:%M') - datetime(1970, 1, 1)).total_seconds()).values#时间格式是/分开的用这个
# 将时间列转换为秒级时间戳，并存储在名为t的NumPy数组中
#更多模型咸鱼搜索机器学习之心，支持模型定制
noisy_signal = data["A"].values # 将"A"列的数据存储在名为noisy_signal的NumPy数组中
wavelet = "db4" # 使用 Daubechies 4 小波
level = 6 # 设置分解层数
coeffs = pywt.wavedec(noisy_signal, wavelet, level=level)
# 对noisy_signal进行小波分解，使用Daubechies 4小波，分解层数为6，将分解系数存储在名为coeffs的列表中
#更多模型咸鱼搜索机器学习之心，支持模型定制
sigma = np.median(np.abs(coeffs[-1])) / 0.6745
# 计算最后一层分解系数的绝对值的中位数
threshold = sigma * np.sqrt(2 * np.log(len(noisy_signal)))
# 根据阈值公式计算阈值
#更多模型咸鱼搜索机器学习之心，支持模型定制
# 应用阈值
coeffs_thresholded = [coeffs[0]] + [pywt.threshold(coeff, threshold) for coeff in coeffs[1:]]
# 重构信号
reconstructed_signal = pywt.waverec(coeffs_thresholded, wavelet)
# 绘制原始信号和重构信号
plt.figure(figsize=(12, 3)) # 创建一个图形窗口，大小为12x3
plt.plot(t, noisy_signal, label="Original Signal") # 绘制原始信号
plt.plot(t, reconstructed_signal, label="Reconstructed Signal (Denoised)") # 绘制去噪后的重构信号
plt.xlabel("Time") # 设置x轴标签为"Time"
plt.ylabel("Amplitude") # 设置y轴标签为"Amplitude"
plt.legend() # 添加图例
#更多模型咸鱼搜索机器学习之心，支持模型定制
# 调整 x 轴刻度密度
x_ticks = np.linspace(t[0], t[-1], 10) # 在时间范围内生成10个刻度点
plt.xticks(x_ticks, [datetime.fromtimestamp(x).strftime('%Y/%m/%d %H:%M') for x in x_ticks], rotation=45) # 设置x轴刻度及标签为对应的时间字符串，并进行45度旋转
plt.show() # 显示图像
# 绘制分解模态
fig, axs = plt.subplots(level + 1, 2, figsize=(12, 2*(level+1)), sharex='col') # 创建子图，行数为level+1，列数为2，图像大小为12x2*(level+1)，共享x轴
for i in range(level + 1):
component = pywt.waverec(coeffs[:i] + [coeffs[i]] + [np.zeros_like(c) for c in coeffs[i+1:]], wavelet) # 重构第i个分解分量
axs[i, 0].plot(t, component) # 绘制分解分量的图像
axs[i, 0].set_title(f"Component {i+1}") # 设置子图标题为"Component i+1"
axs[i, 0].set_ylabel("Amplitude") # 设置y轴标签为"Amplitude"
if i == level:
axs[i, 0].set_xlabel("Time") # 当i为level时，设置x轴标签为"Time"
#更多模型咸鱼搜索机器学习之心，支持模型定制
# 调整 x 轴刻度密度
x_ticks = np.linspace(t[0], t[-1], 10) # 在时间范围内生成10个刻度点
for ax in axs.flat:
ax.set_xticks(x_ticks) # 设置x轴刻度为生成的刻度点
ax.set_xticklabels([datetime.fromtimestamp(x).strftime('%Y/%m/%d %H:%M') for x in x_ticks], rotation=45) # 设置x轴刻度标签为对应的时间字符串，并进行45度旋转
# 计算每个分量的功率谱密度并绘制频谱图
for i in range(level + 1):
component = pywt.waverec(coeffs[:i] + [coeffs[i]] + [np.zeros_like(c) for c in coeffs[i+1:]], wavelet) # 重构第i个分解分量
f_component, Pxx_component = welch(component, nperseg=200) # 计算分解分量的功率谱密度
axs[i, 1].semilogy(f_component, Pxx_component) # 绘制频谱图，使用对数坐标
axs[i, 1].set_title(f"Frequency Spectrum of Component {i+1}") # 设置子图标题为"Frequency Spectrum of Component i+1"
axs[i, 1].set_ylabel("Power Spectral Density") # 设置y轴标签为"Power Spectral Density"
if i == level:
axs[i, 1].set_xlabel("Frequency [Hz]") # 当i为level时，设置x轴标签为"Frequency [Hz]"
plt.tight_layout() # 调整子图的布局，使其紧凑显示
plt.show() # 显示图像
# In[ ]:
